Combining neutrino experimental light-curves for pointing to the next Galactic Core-Collapse Supernova

The multi-messenger observation of the next galactic core-collapse supernova will shed light on the different physical processes involved in these energetic explosions. Good timing and pointing capabilities of neutrino detectors would help in the search for an electromagnetic or gravitational-wave counterparts. An approach for the determination of the arrival time delay of the neutrino signal at different experiments using a direct detected neutrino light-curve matching is discussed. A simplified supernova model and detector simulation are used for its application. The arrival time delay and its uncertainty between two neutrino detectors are estimated with chi-square and cross-correlation methods. The direct comparison of the detected light-curves offers the advantage to be model-independent. Millisecond time resolution on the arrival time delay at two different detectors is needed. Using the computed time delay between different combinations of currently operational and future detectors, a triangulation method is used to infer the supernova localisation in the sky. The combination of IceCube, Hyper-Kamiokande, JUNO and KM3NeT/ARCA provides a 90% confidence area of about 340 deg$^2$. These low-latency analysis methods can be implemented in the SNEWS alert system.


Introduction
The two dozen of neutrinos observed from the SN1987A explosion at the Large Magellanic Cloud indicate that Core-Collapse Supernovae (CCSN) are sources of 1-100 MeV neutrinos, with a total energy emitted on the order of 3 × 10 53 erg. In such explosive phenomena, 99% of the total gravitational binding energy of the stellar collapse is released through neutrinos [1]. Neutrinos leave the stellar core several hours before photons. Fast and precise direction reconstruction of the neutrino flux is important for a potential multimessenger follow-up.
In this work, several water Cherenkov detectors are considered: the underground detectors, Super-Kamiokande [2] and Hyper-Kamiokande [3], and the highenergy neutrino telescopes, IceCube [4] and KM3NeT [5]. The main interaction channel in water Cherenkov detectors is inverse beta decay (IBD) of electron antineutrinos (ν e ) on proton targets. The positrons angular distribution is slightly forward-peaked and high energy events can be selected in order to exploit this directionality [6]. A good pointing for CCSN is hardly achievable in the water Cherenkov detectors with such technique. The JUNO scintillator detector is sensitive to both IBD and elastic scattering (ES) [7]. Moreover, for the IBD channel JUNO offers the ability to identify positrons and neutrons and reconstruct their positions. The direction along the line connecting the positions of both products can be used to infer the neutrino direction [8]. Detectors sensitive to ES interactions can provide information on the CCSN localisation. Nowadays, Super-Kamiokande is the only running detector with enough sensitivity to the ES channel to be able to point by itself to the source [9].
By exploiting the time delay between the arrival of the signal at different detector sites, it is possible to use a triangulation method to infer the source localisation on the sky [10,11,12,13,14,15]. In the triangulation method proposed in [10], the uncertainty on the arrival time delay between a pair of detectors is inferred from the total number of events and the CCSN duration. A Kolmogorov test for direct light-curve comparison between detectors was later mentioned in [11], with the conclusion that low event statistics detectors available at that time were not useful for triangulation. In [12], the triangulation method is revisited and a rough estimate of the arrival time uncertainty for each detector is computed assuming a generic neutrino light-curve with an exponential rise. In [13], a more detailed lightcurve template is used. The use of the time delay estimate between the first detected events in each experiment is proposed in [15]. The latter method can be implemented for real-time CCSN localisation. In order to reach a good accuracy, an almost backgroundfree experiment is required. This implies that the use of large volume Cherenkov neutrino telescopes, as Ice-Cube and KM3NeT, is more difficult. Timing with the first events is earlier proposed in [14] together with the exponential rise fit.
In this work, a model-independent approach that relies on matching the detected neutrino light-curves is elaborated. Such an approach requires data sharing between the detectors, which can be done within the SNEWS global network [16].
The paper is organized as follows. The simulation framework to model CCSN neutrino detection rate at different detectors is introduced in Section 2. In Section 3, the statistical methods used to estimate the time difference of the neutrino arrival at the different sites are described. In Section 4, the arrival time delay resolution for different detector pairs is discussed. Different combinations of three and four detectors are then used to evaluate the CCSN localisation uncertainty area. Conclusions are provided in Section 5.

Simulation
The CCSN neutrino detection rates at different observatories is estimated with a simple simulation described in this section. It includes a simplified neutrino luminosity curve and the detector parameters describing their detection efficiency and background.

Supernova neutrino fluxes
The all-flavour neutrino total luminosity is set to L 0 ν = 3 × 10 53 erg, as in [7]. It is equally divided among all neutrino flavours and between neutrinos and antineutrinos. This is consistent with the more accurate estimate of theν e fraction of about 0.14 that includes oscillations in the supernova mantle [6].
The time evolution of the neutrino luminosity is taken from [17] and it is described forν e as: where N is used to normalise the function as follows: For these studies the chosen parameters are: t a = 0.035 ms, t c = 0.2 ms, n a = 2, n p = 20 and n c = 1.5. The light-curve simulated using these parameters approximately follows the predicted accretion phaseν e luminosity in [18]. The chosen value of n c = 1.5 guarantees the convergence of the luminosity integral. The luminosity curve with the assumed parameters is shown in Figure 1. The energy distribution can be described by a quasithermal distribution [18]: which depends on the average neutrino energy,Ẽ ν , and the spectral pinching shape parameter, α. Both parameters are generally varying with time. The typical range of the pinching shape parameter is 2 α 5 [18].
Only the variation in time of the detection rates is relevant for this study. Therefore, the simulation follows the luminosity time evolution with a steady energy neutrino spectrum withẼ ν = 14 MeV and α=3. The neutrino luminosity can be converted into a flux dividing it by the average neutrino energy. Assuming an isotropic neutrino emission from a source at a distance d, the differential neutrino flux is:

CCSN neutrino detection rates
As it is the main channel for the considered water detectors, only the inverse beta decay interaction channel is simulated. The instantaneous event rate in the detector is estimated as the product of the differential neutrino flux, dΦ dEν , the IBD cross-section, σ IBD [6], the number of targets per unit volume, 2 ρN A µ H 2 O , the detection efficiency, det , and the detector volume V det : The factor 2 comes from the two hydrogen atoms in a water molecule, µ = 18 g/mol is the molar mass of water, ρ is the density of water and N A is the Avogadro number.
The detector properties are converted into the detector effective mass, M eff , in the following way: The detection efficiency depends on the neutrino energy, or more precisely, on the energy of the detectable interaction products. A constant efficiency is assumed above an energy threshold, E min ν . Since this value is around the Cherenkov threshold, and below the energy range where supernova neutrinos are expected, E min ν = 0 is set for simplicity. This assumption removes the energy dependence of the M eff . The detection rate is calculated as: where the simplified conversion parameter, I, is the same for all the detectors. For a CCSN at 10 kpc, this factor becomes: The rates for future scintillator or other non-water detectors can be calculated using Equation 7, where M eff is estimated in water equivalent units and N target is the total number of targets in the detector:

Simplified detector model
The detector model is described by two parameters: the supernova detection effective mass, M eff , and the background rate, R bg . The signal rates are estimated for each detector from Equation 7. Both the signal and the background rates are translated into an expected number of events per time bin. In order to simulate experimental fluctuations in the detected neutrino lightcurve, the number of events in each time bin is sampled assuming a Poisson distribution. For the Super-Kamiokande and Hyper-Kamiokande detectors, the effective masses are taken from [19,20,21]. For JUNO, the expected number of proton targets is N target = 1.5 × 10 33 [22], similar to Super-Kamiokande. The same effective mass in water equivalent units, M eff = 22.5 kton, is used for both detec-tors. In this study, the background rates are negligible for Super-Kamiokande, Hyper-Kamiokande [23] and JUNO [22].
The assumed IceCube detector effective volume is 3.5 Mton [24]. The detector consists of 5160 optical modules, each containing one 10-inch photomultiplier. The background rate per optical module is about 300 Hz, taking an average value between standard and high efficiency optical modules.
The optical module of KM3NeT consists of 31 3inch photomultipliers. Using nanosecond scale coincidences between the photomultipliers on the same optical module, the variable background from the bioluminescence can be suppressed [25]. The remaining background contribution is mostly coming from 40 K decays in sea water and it amounts to a total rate of R K40 OM ∼500 Hz per optical module [26]. The effective mass of the KM3NeT detectors is estimated as M eff ≈ 180 kton for KM3NeT/ARCA and 90 kton for KM3NeT/ORCA from [25, Fig. 4, left]. The ARCA and ORCA detectors consist of 4140 and 2070 optical modules, respectively. The coincidence selection translates into a lower effective CCSN detection volume for KM3NeT/ARCA compared to IceCube, even if both have an instrumented volume of ∼1 km 3 .
The effective mass and the background rate for the different detectors used in this work are summarised in Table 1. Some examples of the simulated detector light-curves are provided in Figure 2.

Description of the methods
Two methods to estimate the time difference of the neutrino arrival times between two detectors have been investigated: the chi-square (Section 3.1) and the crosscorrelation (Section 3.2). Performance estimation with large number of simulations is described in Section 3.3 using a simplified CCSN model. Bootstrapping procedure using detected light-curve is introduced in Section 3.4. The triangulation method for the source localisation is reviewed in Section 3.5.

Chi-square method
A method based on χ 2 minimisation has been developed and tested to infer the delay between two lightcurves. This technique is one of the traditional ways to verify the compatibility of two distributions and for parameter estimation. In this case, the estimated parameter is the time delay between the two detected  Assuming a normal distribution of the number of events in each time bin the χ 2 expression is defined as follows: where τ is the time shift between the detected neutrino light-curves, n ti−τ is the number of observed events by the first detector in the time bin t i − τ , m ti is the number of observed events by the second detector in the time bin t i , E(n ti−τ − m ti ) and V (n ti−τ − m ti ) are the expectation value and the variance of the difference in the number of events, respectively. For two normal distributions, the variance of the difference between the number of events corresponds to the sum of their squared standard deviations: The value of τ = T match 0 obtained at the χ 2 minimum provides the best estimate of the true time delay between the two detectors, T 0 .
A different expression of the χ 2 was used in [13]. It assumes that the number of events in each time bin follows a Poisson distribution with a mean value equal to the average expected number of events computed from a CCSN neutrino emission model. This Poisson χ 2 expression cannot be used to compare two experimental light-curves since it requires a known expectation value. Since the number of detected events follows a Poisson distribution, the bin size should be optimized in our method to make the Gaussian approximation valid.
If the light-curves are normalised in a way that signal and background expectations are almost identical for both detectors, then the expectation value E(n ti−τ −m ti ) can be considered null for any t i when τ corresponds to the true shift. To achieve such normalisation, the mean detector background is subtracted from the light-curve so that all combined detectors have a common baseline centred at 0. The background expectation value can be taken from a time window in the detected neutrino light-curve where CCSN signal is not present (off-signal zone) or it can be provided by each experiment based on a longer detector monitoring. In this study, an off-signal zone of 1 s duration is chosen for the background estimate. In order to account for the different detector efficiencies, a normalisation to the detector effective mass is applied. Alternatively, this effect can be taken into account by normalising each light-curve after background subtraction to have unit integral.
The signal arrival time is not known a priori in this model independent analysis. Therefore, a reference time is evaluated to define the window for the χ 2 calculation. The time window of [-300, 300] ms centred at the maximum of the detected light-curve from one of the detectors is chosen for this analysis. The interval covers the transition between the background and the CCSN neutrino signal as well as the accretion phase, for which most of the neutrino emission is expected [18]. Using a time window too long may lead to a degradation of the performance since the optimisation will be more sensitive to background fluctuations. The χ 2 calculation window is fixed with respect to the less performing detector to preserve its background statistics during τ scan. The choice of the fixed detector and the detector for the reference time definition is done to minimize the estimated time delay uncertainty. The procedure for the uncertainty estimation is given in Section 3.3.
The detected light-curves are provided as histograms with a fixed bin size. A bin width of 0.1 ms is chosen, and the same value is used for the time shift step τ , in order to reach O(ms) required resolution. The numbers of events, n ti−τ and m ti , are calculated by summing the events of contiguous bins. The resulting effective bin size is optimised for each detector pair to reach the minimum of arrival time delay uncertainty, δt. Steps of t i+1 − t i can be optimised in order to increase the calculation speed.

Normalised cross-correlation
The cross-correlation can be used for matching two sets of time series [27] or for matching a time series with a template model [28]. The discrete cross-correlation is defined as: where n ti and m ti−τ are the number of observed events by the first and the second detector in time bins t i and t i − τ , respectively. N is the number of bins in the search window [t min , t max ], and τ is the time delay between the two light-curves. The function will present a maximum at τ = T match 0 , allowing to estimate the time delay between the two detectors, T 0 .
In order to account for the different effective masses and background rates of each detector, the zero-normalised cross-correlation is used. Each curve is normalised by subtracting its mean value and scaled by its standard deviation [27], which can be computed in the search window: One of the advantages of this method is that fast Fourier transformations can be used to speed up the calculations and the improvement of the response time can be significant for a large number of bins compared to the chi-square method.

Procedure for the performance estimation
Simulations of the different detected neutrino lightcurves are performed following the model described in Section 2. An artificial delay of the neutrino signal observed between two sites, T 0 , is applied to the first light-curve. The result is not expected to change when exchanging the two detectors. The best estimate T match 0 is inferred with the proposed methods (see Sections 3.1 and 3.2). Finally, T match 0 , is compared to the true T 0 value.
The distribution of T match 0 − T 0 is built from a large number of realisations of the simulated light-curves. To confirm the absence of systematic effects, the distribution of T match 0 − T 0 should be compatible with a normal distribution with mean zero and standard deviation independent of T 0 . The width of the distribution provides an estimate of the uncertainty on the T 0 measurement, δt.

Bootstrapping parameters tuning
The proposed methods present some parameters that can be slightly tuned for different CCSN models, distances to the source and detectors. In most cases, it was verified that the degradation of the time precision is not significant and some common set of parameters can be identified for future combined analysis.
In order to reach the best performance once the supernova is detected, the following bootstrapping procedure can be set up. The detected neutrino light-curve from the best performing detector is used as the model. With this model, the detected neutrino light-curves for all the detectors can be simulated and the procedure to estimate the performance described in Section 3.3 can be repeated in order to tune the parameters.

The triangulation method
The time difference between the CCSN neutrino signal arrival at two detectors located at r i and r j can be expressed as: t ij = ( r i − r j ) · n/c, where n is the unit vector that indicates the emission direction. This vector is calculated in the geographic horizontal coordinate system from the CCSN right ascension, α, declination, δ, and the Greenwich mean sidereal time expressed as angle, γ: On March 21 2000 at noon the J2000.0 equatorial coordinate system matches with the geographic one since γ = 0 • . For this time Equation 15 is simplified to the same expression used in [13]. The position of the detector k can be inferred from its latitude (φ k ) and longitude (λ k ) angles, and the Earth radius (R Earth ): The probability that the scanned angles (α, δ) coincide with the equatorial coordinates of the CCSN is given by the following χ 2 function: assuming that there is no systematic shift in the T match 0,ij determination. The minimum of the function gives the best estimate for the angles (α, δ) of the searched CCSN location in the sky. From Equation 17, one can note that the performance depends on the uncertainty of the measured time delay δt ij of each detector pair. Different detector pairs can be combined into a total χ 2 by summing each contribution: The χ 2 (α, δ) function is converted into a p-value, p(α, δ) = p(χ 2 (α, δ) ≤ χ 2 min ), which is the probability of observing a χ 2 smaller or equal to the minimum χ 2 value obtained scanning all possible directions (α, δ). The 90% confidence level (C.L.) error box of the source localization area is determined as a collection of all points on the sky with a p-value p(α, δ) < 0.9.

Results and performance comparisons
In this section, the estimated resolution of the arrival time delay for different detector pairs is shown. Different combinations of three and four detectors are then used to estimate the resulting CCSN localisation uncertainty area.

Time uncertainty results
The results of the chi-square method are given for two different light-curve normalisations in Tables 2 and Table 3. For the first normalisation, the true background value is assumed for the baseline subtraction and the scaling is done according the true effective mass. This corresponds to the ideal case in which the detector efficiency and the background estimations are known a priori. For the second one, the background expectation value is computed as the average rate in a 1 s off-signal region. The light-curve area in the time window of [-300, 300] ms around the light-curve maximum is normalised to one. By comparing Table 2 and Table 3 results, the realistic experimental curve normalisation gives similar performances compared to the ideal case. This justifies the proposed window for background estimation and the window used for light-curve normalisation. It is verified that the mean of the T match 0 − T 0 distribution is compatible with zero within the statistical uncertainties. The δt obtained from the simulations with a fixed and random delay times are compatible. This confirms that for the assumed model and detector response the chi-square method provides an unbiased estimation of the time delay between the signal arrival at the two detector sites.
For KM3NeT, the effective mass corresponding to the ARCA detector is used for the performance estimation. The maximum signal time delay between the ARCA and ORCA sites, (∼3 ms), is on the same order as the estimated uncertainties for any of the combinations involving the ARCA detector. This prevents a simple merging of the two KM3NeT light-curves.
The results of the chi-square method shown in Tables 2 and 3 are computed using the optimized step and effective bin sizes of 50 ms. Using 10 ms for the step and the bin size for ARCA/IceCube and ARCA/Hyper-Kamiokande combinations provides slightly better results, reaching the value of 6.20±0.15 and 6.30±0.15 ms, respectively. Adding up to 3 Hz of background rates does not decrease the performance for Super-Kamiokande and JUNO detectors.
The performance of the cross-correlation method using the same search time window as for the chi-square method is shown in Table 4. The optimal effective bin size is found to be 10 ms for all combinations. The results are compatible with the chi-square method, so cross-correlation represents a viable alternative.

Results of the triangulation: localisation skymaps
Using the uncertainties estimated for each detector combination, the triangulation algorithm is applied to reconstruct the CCSN location on the sky. To estimate the performance of different detector combinations, a CCSN on vernal equinox at noon is assumed for simplicity. The source direction of the Galactic Centre with equatorial coordinates (α 0 , δ 0 ) = (−94.4 • , −28.9 • ) is chosen as an example. The confidence area is computed from the scan of t ij (α, δ) evaluated on a grid of 256 directions on the sky looking for the best match, assuming that the true value is assigned to the fitted value T match 0,ij = t ij (α 0 , δ 0 ). The source confidence areas are shown in Figures 3  and 4 using the combination of four and three detectors, respectively. These results are obtained using the uncertainties from the chi-square method in Table 3. The size of the error box at 68 and 90% C.L. are provided in Table 5. The results in Table 5 indicate that a favourable position of a detector with respect to the source location Table 2: Uncertainty δt in milliseconds obtained with the chi-square method using ideal normalization of the detector neutrino light-curves. The detector pairs are listed in row and column names. The arrow points to the detector name that is used for the light-curve peak definition, which is also shifted during the χ 2 scan.  Table 3: Uncertainty δt in milliseconds obtained with the chi-square method using average background subtraction and unity normalization of the detector neutrino light-curves. The detector pairs are listed in row and column names. The arrows point to the detector name that is used for the light-curve peak definition, which is also shifted during the χ 2 scan.   and other detectors may compensate for worse time resolution, for example when replacing KM3NeT/ARCA with JUNO in the combination of three detectors. The confidence areas may be further reduced considering their intersection with the Galactic Plane.
The results of our work can be compared to the latest triangulation studies. In [13], the estimate of the 68% C.L. area is ∼66 deg 2 . This result is better than the one obtained in our work (220 deg 2 ), but more than four detectors were combined using IBD and ES channels and the uncertainty estimate relies on the matching with a light-curve template known a priori. The result of our work represents a model independent data analysis proposal. In [14] the timing with the first IBD events has the best performance comparing to the exponential rise fit. In the latest results with the first events observation method [15], the time delay uncertainty for Super-Kamiokande and JUNO combination is 5.7 ms, which is larger than 2.8 ms estimated in our work. This method requires an evaluation and correction for several biases due to background rates and the steepness of the luminosity curve rise. Biases in the method proposed here may appear due to the detector efficiency varying with the neutrino energy and in time, or with respect to the event rates. A proper estimation of such effects is possible considering simulations with a detailed emission model and precise detector parameters. This new method relies on the agreement among the different collaborations for making their full lightcurve available to SNEWS while the method in [15] only requires sending the time information of the first events.
The results of this work can also be compared with the expected performance of Super-Kamiokande using the directionality information from the elastic scattering channel. The 68% C.L. area for the combination of four detectors in this work is smaller than the area expected with the actual Super-Kamiokande configuration (∼314 deg 2 [9]). With Gadolinium doping, this area might be reduced down to ∼13 deg 2 [9]. The expected CCSN 68% C.L. area for the JUNO IBD events reconstruction analysis will be better than 254 deg 2 , comparable with our result, based on the matching of light-curves from four detectors. The triangulation method can be proposed as a low latency analysis. The confidence areas from Super-Kamiokande, JUNO and this triangulation analysis are independent and can be combined in order to obtain a joint refined measure-

Conclusion
Detectors in current and future operation will be capable of collecting a significant number of neutrinos from the next galactic CCSN explosion. This will allow for detailed studies of the time profile. The determination of the neutrino arrival time is crucial for the source localisation that may help for a potential multi-messenger follow-up.
The detected neutrino light-curves can be used to deduce the delay in the arrival time of the supernova emission at different sites on Earth. This paper de-scribes the chi-square and cross-correlation methods used to estimate such delays and their uncertainties. Since a direct comparison of the experimental curves is used, these methods do not rely on any model for the alert algorithm implementation.
Supernova emission consists of a neutrino flux with a different light-curve for each flavour. Therefore, the compared experimental light-curves should consist of events detected through the same neutrino interaction process. This is the case for most detectors sensitive to the electron anti-neutrino interactions via inverse beta decay in water. The method can be extended to detectors sensitive to different channels accounting for systematic biases, which can be estimated on a model dependent basis.
The methods are tested using a time dependent parametrisation of the neutrino luminosity. The detectors are described with two parameters, the effective mass and the background rate. The signal and background are sampled assuming Poisson distributions. These simulations may be further improved by the respective collaborations using accurate detector simulations. Detailed CCSN fluxes can be used in the simulation to estimate the performance of the methods for any given model.
Merging the time delay information between several detector sites allows to infer the supernova localisation. The 90% confidence areas are about 300-400 deg 2 when combining IceCube, Hyper-Kamiokande, JUNO and KM3NeT/ARCA detectors. Such analysis can be performed in real-time within the framework of the SNEWS alert system. This location uncertainty can be reduced further by intersecting this area with the CCSN progenitors distribution in the Milky Way and combining with the confidence areas from the Super-Kamiokande and JUNO detectors performing standalone source localisation.
Using a bootstrapping strategy, the algorithm parameters can be optimised directly on data to improve the time resolution and to localise the source with better accuracy.