Constraining Absolute Neutrino Masses via Detection of Galactic Supernova Neutrinos at JUNO

A high-statistics measurement of the neutrinos from a galactic core-collapse supernova is extremely important for understanding the explosion mechanism, and studying the intrinsic properties of neutrinos themselves. In this paper, we explore the possibility to constrain the absolute scale of neutrino masses $m^{}_\nu$ via the detection of galactic supernova neutrinos at the Jiangmen Underground Neutrino Observatory (JUNO) with a 20 kiloton liquid-scintillator detector. In assumption of a nearly-degenerate neutrino mass spectrum and a normal mass ordering, the upper bound on the absolute neutrino mass is found to be $m^{}_\nu<(0.83 \pm 0.24)~{\rm eV}$ at the 95% confidence level for a typical galactic supernova at a distance of 10 kpc, where the mean value and standard deviation are shown to account for statistical fluctuations. For comparison, we find that the bound in the Super-Kamiokande experiment is $m^{}_\nu<(0.94 \pm 0.28)~{\rm eV}$ at the same confidence level. However, the upper bound will be relaxed when the model parameters characterizing the time structure of supernova neutrino fluxes are not exactly known, and when the neutrino mass ordering is inverted.


I. INTRODUCTION
The neutrino burst from Supernova (SN) 1987A in the Large Magellanic Cloud was clearly recorded in the Kamiokande-II [1], Irvine-Michigan-Brookhaven (IMB) [2], and Baksan [3] experiments. Although only twenty-four neutrino events in total were observed, we have acquired very useful information about the explosion mechanism of core-collapse SNe and the intrinsic properties of neutrinos themselves [4]. As pointed out by Zatsepin long time ago, neutrinos from SNe can be used to constrain the absolute neutrino masses through their delayed flight time [5]. In fact, the observed neutrino events associated with SN 1987A have been reanalyzed by Loredo and Lamb in great detail [6], and they have obtained a restrictive limit m ν < 5.7 eV at 95% confidence level (CL) 1 , significantly improving the even earlier analyses in the literature [7][8][9][10][11][12].
Recently, the upper limit on neutrino masses from the neutrino detection of a galactic SN has been derived in Ref. [13] for a large water Cherenkov detector, such as the Super-Kamiokande experiment (SK) [14]. It has been shown that m ν < 0.8 eV at the 95% CL can be achieved for a typical galactic SN at a distance of 10 kpc. An upper bound in the sub-eV region has also been claimed in Refs. [15] and [16] for future large water Cherenkov and scintillator detectors (e.g., Hyper-Kamiokande [17] and LENA [18]). In this paper, we investigate the sensitivity of JUNO to absolute neutrino masses via the detection of galactic SN neutrinos. The motivation for such an investigation is as follows: • For a galactic SN neutrino burst, the JUNO detector is expected to collect about 10 4 events, mainly in the inverse beta decay channel ν e + p → e + + n (IBD). The total number of IBD events is comparable to that in the SK. Furthermore, the JUNO detector has a lower energy threshold and a better energy resolution, which are crucial advantages in measuring the time-delay effects of low-energy SN neutrinos and should lead to an improvement of the neutrino mass bound.
• An independent probe of absolute neutrino masses in the sub-eV region is necessary and desirable. Thus far, the most stringent constraint on the effective neutrino 1 For simplicity, we hereafter assume the neutrino mass spectrum to be quasi-degenerate, i.e., m 1 ≈ m 2 ≈ m 3 ≡ m ν . Hence all the experimental bounds on absolute neutrino masses can be translated to those on the neutrino mass scale m ν . In particular, this assumption is well justified for SN neutrinos, since current experiments are only sensitive to neutrino masses lying in the quasi-degenerate region. mass in tritium beta decays is m ν < 2.1 eV at 95% CL from the Troitsk Collaboration [19]. The future KATRIN [20] and Project 8 [21] experiments are planning to improve the present limit by one order of magnitude, reaching the sub-eV level.
Useful information on neutrino masses can also be obtained from the detection of neutrinoless double-beta decays, if neutrinos are Majorana partciles, i.e., neutrinos are their own antiparticles. Given the current bound on the effective neutrino mass [22], where U ei (for i = 1, 2, 3) stand for the three elements in the first row of neutrino mixing matrix U, one can derive m ν < ∼ (0.4 · · · 1.5) eV, depending on the neutrino mass hierarchy, neutrino mixing angles and two unknown Majorana CP-violating phases.
• Cosmological observations give rise to a very tight bound on the sum of three neutrino masses, namely Σ ≡ m 1 + m 2 + m 3 < 0.23 eV at the 95% CL [23]. However, it should be noted that the cosmological bound can be quite different if the systematic uncertainties for some data sets are taken into account or if an extension of the minimal standard model of cosmology is considered [24]. Therefore, it is interesting to investigate whether a high-statistics observation of galactic SN neutrinos could provide an independent and competitive bound on m ν .
Although following closely the general approach proposed in Ref. [13], our analysis differs from previous works in several aspects. First, we compare the mass bound from JUNO with that from SK, and point out that JUNO with a lower energy threshold and a higher energy resolution can indeed improve the bound. Due to the limited statistics in the low-energy and early-time region, however, such an improvement is moderate. Second, we have simulated a large number of experiments to study the statistical uncertainties of neutrino mass bound.
The upper bound at the 95% CL for JUNO turns out be m ν < (0.83 ± 0.24) eV, while that for SK m ν < (0.94 ± 0.28) eV. Third, we point out that both the starting time t s and the rising-time interval τ r for neutrino emission have crucial impact on the neutrino mass bound.
If these two SN model parameters are set to be free, the upper bound will be relaxed to m ν < (1.12 ± 0.33) eV for JUNO, and m ν < (1.49 ± 0.42) eV for SK.
The remaining part of this work is organized as follows. In Sec. II, we briefly review a simple model of the SN ν e flux, and discuss the general approach to constrain absolute neutrino mass scale by observing the time-delay effects of SN neutrinos. Sec. III is devoted to our simulation results for JUNO, and a comparison is made with SK. Furthermore, we investigate the impact of model parameters on the mass bound, and consider the statistical uncertainty by simulating a large number of experiments. Finally we conclude and summarize our main results in Sec. IV.

II. SN NEUTRINOS AND MASS BOUND
Since neutrinos are massive, their flight time from the SN to the detector at the Earth will be delayed, compared to that of massless particles [5]. For the neutrinos from a galactic core-collapse SN, the time delay can be written as where E ν is the neutrino energy, and D is the distance between the SN and the detector.
Thus, a time delay at the mini-second level is expected for the neutrinos of eV masses from a typical galactic SN at D = 10 kpc. If the neutrino energies and arrival time can be precisely measured, the shifts of neutrino events in the time distribution will signify nonzero neutrino masses. However, the average energy and luminosity of SN neutrinos evolve in time, which complicates the extraction of neutrino mass information from experimental observations. In order to derive a mass bound, one has to model the time evolution of neutrino energies and fluxes, and take into account neutrino flavor conversions in the propagation from the SN core to the detector.

A. Parametrized Neutrino Flux
In a core-collapse SN, neutrinos and antineutrinos of three different flavors can be produced [25]. Along with the development of SN explosion, neutrino emission can be described As the cross section of IBD is much larger than those of other channels for SN neutrinos of typical energies, we only concentrate on the ν e flux Φ 0 ν e in the accretion and cooling phases. In Ref. [26], a simple parametrization of Φ 0 ν e is presented to capture essential physics of SN neutrino production and the main features of numerical simulations. In the cooling phase, the ν e flux Φ 0 c is parametrized by three model parameters: the initial temperature T c , the radius of neutrino sphere R c , and the cooling time scale τ c . In the accretion phase, one has to model the time evolution of neutron number and positron temperature in order to figure out the ν e flux Φ 0 a . This can be done by introducing the accretion time scale τ a , and requiring the resultant neutrino energy and luminosity to approximately follow numerical simulations. In addition, the initial number of neutrons depends on an initial accreting mass M a , and a thermal energy spectrum of positrons depending on an initial temperature T a is assumed. Put all together, the total flux is [26,27] where f r (t) = 1 − exp(−t/τ r ) with the rising time scale τ r further introduces an early-time fine structure, and j k (t) = exp[−(t/τ a ) k ] (with k being an integer) is the time function interpolating the accretion and cooling phases of neutrino emission. In our calculations, k = 2 will be chosen, and the analytical expressions of parametrized fluxes Φ 0 a (t, E ν ) and Φ 0 c (t, E ν ), which can be found in Ref. [26], will be implemented.

B. Neutrino Flavor Conversions
When propagating from the SN core to the envelope, neutrinos experience three different stages of flavor oscillations. First, the matter density in the SN core is so high that the coherent flavor conversions will be interrupted by the frequent scattering of neutrinos on matter particles. Therefore, the lepton flavors are indeed conserved in the dense core [28].
Even for the heavy-lepton flavors ν µ and ν τ , the one-loop corrections to neutrino refractive index are significant enough to enhance matter effects, suppressing the flavor oscillations.
Second, from the neutrino sphere to the radius of several hundred kilometers, neutrino number densities are even higher than or comparable to the electron number density of matter, and thus the neutrino-neutrino refraction may lead to an instability in the flavor space.
As a consequence, collective neutrino oscillations may take place, and cause remarkable changes in the neutrino energy spectra [29,30]. The flavor instability induced by neutrino self-interaction is currently an unresolved problem [31]. It has been argued that the high matter density during the accretion phase could completely suppress collective neutrino oscillations [32][33][34]. More recently, the axial symmetry usually assumed in the previous studies of SN neutrino oscillations has been found to be not satisfied by the solutions to the equations of neutrino flavor evolution [35][36][37][38][39], and the flavor instability needs more dedicated investigations. As for the time-delay effects induced by neutrino masses, only the low-energy neutrinos from the early accretion phase are quite relevant. In this case, we expect that the dense matter will highly suppress the growth of flavor instability, and thus the collective neutrino oscillations do not take place at all. Finally, the resonant flavor conversions corre- [40] will occur in the SN envelope, where the matter density becomes suitable for the Mikheyev-Smirnov-Wolfenstein effects to be crucially important [41,42].
The Mikheyev-Smirnov-Wolfenstein effects on neutrino flavor oscillations in the SN envelope have been examined in detail in Ref. [43]. Taking account of a relatively large θ 13 (e.g., sin 2 θ 13 ≈ 0.024), which has been precisely measured in Daya Bay and other reactor experiments [44][45][46][47], one can find that the transition for the ∆m 2 31 resonance is perfectly adiabatic. Therefore, if the neutrino mass ordering is normal, i.e., ∆m 2 31 > 0 or m 1 < m 2 < m 3 , the flux of ν e on the surface of SN is given by with sin 2 θ 12 ≈ 0.302, where Φ 0 ν e and Φ 0 ν x denote the fluxes of ν e and ν x in the case of no flavor oscillations, respectively. If the neutrino mass ordering is inverted, i.e., ∆m 2 31 < 0 or equivalently m 3 < m 1 < m 2 , we have Φ ν e = Φ 0 ν x . Note that Φ 0 ν e receives the contributions from both accretion and cooling phases, as indicated by Eq. (2), while Φ 0 ν x is assumed to be vanishing in the accretion phase [26]. Thus, the total flux of ν e in the case of inverted neutrino mass ordering is much lower than that in the case of normal mass ordering, leading to a reduced number of IBD events and less restrictive bound on neutrino masses. In the following discussions, we shall focus on the normal neutrino mass ordering, but it should be noticed that the derived bound is the most optimistic one.
Given the flux Φ ν e (t, E ν ), one can obtain the event rate R(t, E e ) by convolving it with the IBD cross section σ IBD (E ν ), which can be found in Ref. [48]. In the IBD reaction, the positron energy is approximately given by E e ≈ E ν − ∆ with ∆ = m n − m p ≈ 1.293 MeV being the neutron-proton mass difference. More explicitly, we have where N p is the number of target protons in the detector, and the detection efficiency factor η(E e ) = 1 is taken for JUNO, while η(E e ) = 0.98 for SK. Note that the angular distribution of the final-state positron is almost isotropic, since the kinetic energy of positron is much lower than the nucleon mass. The visible energy is E ′ = E e + m e ≈ E ν − 0.8 MeV, which will be observed as E in the detector due to a finite energy resolution δE. For simplicity, the Gaussian distribution G(E ′ , E; δE) for the energy smearing is assumed.
Our strategy to generate neutrino events and derive the bound on absolute neutrino masses is as follows. First, the event rate R(t, E e ) is used as a target distribution function to randomly produce N neutrino events, represented by (t i , E i ) for i = 1, 2, · · · , N, where the observed energy E i has further been generated according to the Gaussian smearing function G(E e + m e , E; δE) for each event. Second, we implement the neutrino flux model with six astrophysical parameters (i.e., τ a , T a and M a for the accretion phase, τ c , T c and R c for the cooling phase), the rising time τ r for the early time structure, the absolute starting time t s for neutrino emission and the neutrino mass m ν to fit the previously generated neutrino events. In order to make use of the time information of every single event, we define the following likelihood function [26] where t ′ i = t i − ∆t(m ν , E i ν ) − t s stands for the real time when the corresponding neutrino is emitted. Note that a common time needed for massless particles to travel from the SN to the detector has been subtracted from the measured time, and the neutrino energy is constructed from the detected energy via E i ν = E i + 0.8 MeV. Then, the theoretical model parameters can be estimated by χ 2 = −2 ln L, and the test statistic is defined as min . Third, to take into account the statistical fluctuation, we simulate In the following, we shall perform the numerical simulations for JUNO and SK in order to derive the SN bound on absolute neutrino masses. Furthermore, the impact of astrophysical model parameters on the mass bound and the statistical uncertainties are studied.

III. RESULTS AND DISCUSSIONS
For JUNO [49][50][51], a fiducial mass of 20 kiloton liquid scintillator with a proton fraction of 12% is adopted, and the energy resolution is assumed to be δE/E = 0.03 MeV/E. Using the SN neutrino flux model in Sec. II, we obtain the total number of IBD events N a = 4516 in the JUNO detector. In each simulation, the actual number of IBD events is determined according to the Poisson distribution with N a as the expectation value.
In order to generate the artificial data, we first set the astrophysical model parameters as ∆χ 2 0 ≡ χ 2 (m ν = 0) − χ 2 min , which is equivalently the logarithm of the likelihood ratio ∆χ 2 0 = −2 ln L(m ν = 0)/L max . Given each artificial data set, the χ 2 function is minimized to calculate ∆χ 2 0 . Then, a large number of simulations are performed to obtain the distribution function f 0 (∆χ 2 0 ) of the test statistic [52]. The hypothesis will be excluded at the 100(1−α)% CL, if the test statistic is larger than λ(α). The latter is fixed by requiring the integration of f 0 (∆χ 2 0 ) over the range ∆χ 2 0 > λ(α) to be equal to α. In Fig. 1, the cumulative distribution function of ∆χ 2 0 has been shown for JUNO in two different cases. In the first case, the parameters τ r and t s are fixed at their true values. In the second one, they are set to be free and the condition ∆χ 2 0 (τ r , t s ) > λ(α) is satisfied for all possible values of τ r and t s .
For SK [14], we assume that the fiducial mass is 22.5 kt and the proton fraction is 11%.
In addition, the threshold of visible energy is taken as 6.5 MeV, while the energy resolution δE/E = 0.023 + 0.41 MeV/E. Using the same SN neutrino fluxes, one can get the total number of IBD events is N a = 4451, which is very close to the number at JUNO. The distribution of ∆χ 2 0 for SK is also calculated by following the above strategy and the result is shown in Fig. 1, where two cases with fixed and free values of τ r and t s are considered.
For comparison, the true χ 2 distribution with one degree of freedom is depicted in Fig. 1.
To draw a neutrino mass bound at the 95% CL, one usually assumes the χ 2 distribution and requires ∆χ 2 0 > 3.84. As one can observe from Fig. 1, the real distribution of ∆χ 2 0 deviates insignificantly from the χ 2 one, for both JUNO and SK. The critical values of ∆χ 2 0 , above which the hypothesis of m ν = 0 will be excluded at the 95% CL, have been indicated in Fig. 1 and summarized in the caption. Therefore, the neutrino mass bound at the 95% CL in the present work should be interpreted as ∆χ 2 0 ≃ 3.28 (or 3.27) for JUNO, and as ∆χ 2 0 ≃ 3.12 (or 3.02) for SK, for the statistical method with fixed (or free) astrophysical time parameters. For a specific simulation, we calculate ∆χ 2 = χ 2 (m ν ) − χ 2 min and set the 95% CL upper bound by requiring ∆χ 2 to follow the corresponding distribution in Fig. 1.

A. Impact of Model Parameters
To examine the impact of model parameters on the neutrino mass bound at JUNO, we study a single experiment and the fit to the artificial data in more detail. In Fig. 2 Based on the same set of artificial data, we also examine the experimental sensitivity to the starting time of neutrino emission at JUNO. In Fig. 3, the data are fitted by t s in three cases, where the neutrino mass m ν and the rising time τ r are fixed or marginalized over. It is worthwhile to note that t s can never exceed the arrival time of the first observed neutrino event, so the curves in Fig. 3 show a common boundary to the right. In the most optimistic situation, the starting time can be determined within a few ms at 95% CL, which is comparable to the sensitivity at IceCube in reconstructing the bounce time [53]. Another promising way to determine the absolute emission time of SN neutrinos is to observe the gravitational waves associated with the SN explosion [13,54]. If an independent determination of the absolute emission time can be achieved, we may realize the optimistic scenario with fixed t s . However, both scenarios of free and fixed time parameters are considered in our studies.
One can also consider to detect the neutronization burst of ν e from a core-collapse SN, for which the sharply-rising time can be used to probe the time-delay effects. However,   for the current and future huge scintillator detectors, the statistics is limited and the reconstruction of neutrino energies from the elastic neutrino-electron scattering is subject to large uncertainties, because the direction of the final-state electron is not determinable. The charged-current interaction of ν e on the carbon target can be implemented as well, but the energy threshold for this process to occur is as high as 17 MeV, leading to a negligible time delay. Another interesting approach to probe absolute neutrino masses based on the time structure of SN neutrinos is to observe the abrupt halt of neutrino signals when a black hole forms during the accretion phase [55,56].

B. Statistical Uncertainties
As we have already mentioned, the neutrino mass bound from one single simulation is not robust, since it suffers from statistical fluctuations. Therefore, we have carried out a large number of simulations for both JUNO and SK, and the final results have been given as histograms in Fig. 4 and Fig. 5  the model parameters τ r and t s are free parameters, the upper bounds will be relaxed to m ν < (1.12 ± 0.33) eV and m ν < (1.49 ± 0.42) eV, respectively.
Since the energy resolution is better and the threshold is much lower at JUNO, compared to those at SK, we expect a more restrictive upper bound. This is really the case, but the improvement is moderate. The main reason is the limited statistics of the low-energy and early-time neutrino events, for which the time-delay effects are significant.
Finally, it is worthwhile to make some remarks on the ultimate sensitivity of future large detectors to the absolute neutrino masses. For this purpose, we have performed the numerical simulations for artificial experiments, which are similar to JUNO but with different energy thresholds, energy resolutions, or target masses.
• To investigate the impact of energy resolution, we consider two scenarios based on the JUNO setup: (1) no energy smearing; and (2) δE/E = 0.023 + 0.41/ MeV/E, i.e., the same smearing as for SK. The upper bounds on neutrino masses turn out to be m ν < (0.81 ± 0.24) eV and m ν < (0.85 ± 0.25) eV at the 95% CL, respectively. On the other hand, we maintain the energy smearing of JUNO, but assume the energy threshold is 6.5 MeV. In this case, the upper bound m ν < (0.92 ± 0.27) eV becomes much worse than that in the realistic case. Therefore, it is now clear that the energy threshold is the most important reason for the difference between JUNO and SK.
• Concentrating on JUNO, we now increase its fiducial mass by a factor of two, five and ten, leading to a significant increase in the total number of IBD events. In these three cases, the upper bounds are improved to be m ν < (0.67±0.20) eV, (0.52±0.15) eV, and (0.42 ± 0.13) eV, respectively. Note that we have assumed a perfect energy resolution, and that all the model parameters are exactly known. One can observe that increasing the number of neutrino events is the most efficient way to improve neutrino mass bound. Hence, the last bound from ten times larger statistics can be regarded as the ultimate sensitivity of future large detectors.
Although the above simulations are not aimed at any realistic detectors, they help us better understand the difference between JUNO and SK, and clarify the limitation of the entire approach to probe absolute neutrino masses.

IV. SUMMARY
Recent years have seen great progress in experimental neutrino physics. In particular, the smallest neutrino mixing angle has been precisely measured in the Daya Bay reactor neutrino experiment. The primary goals of future neutrino oscillation experiments will be the determination of neutrino mass ordering and the discovery of leptonic CP violation.
Unfortunately, the oscillation experiments are insensitive to the absolute neutrino masses.
As is well known, the tritium beta decays and neutrinoless double-beta decays could provide us with useful information about neutrino masses in the sub-eV region. Currently the tightest bound on the sum of neutrino masses Σ < 0.23 eV comes from the cosmological observations, which however may suffer from large systematic uncertainties. Therefore, an independent way to probe neutrino masses at the sub-eV scale is desirable and important.
In the present work, we consider the possibility to constrain neutrino masses by observing galactic SN neutrinos at the JUNO detector. As the largest scintillator detector over the world, JUNO will register about five thousand neutrino events in the inverse-beta decay channel for a galactic core-collapse SN at a distance of 10 kpc. Since the arrival time and neutrino energy can be well measured at JUNO, the distortion in the time distribution of SN neutrino events caused by the delay of flight time is sensitive to the absolute scale of neutrino masses. Based on a simple but useful model of SN neutrino fluxes and the maximum likelihood analysis, we have carried out a number of simulations to explore the upper bound on absolute neutrino masses at JUNO. It is found that m ν < (0.83 ± 0.24) eV at the 95% CL can be reached, where the mean value and standard deviation are shown to account for the statistical fluctuation. For comparison, we find that the bound in the Super-Kamiokande experiment is m ν < (0.94 ± 0.28) eV at the same CL. Different from the previous works, the impact of astrophysical model parameters on the neutrino mass bound has been emphasized and studied in more detail. Moreover, the statistical uncertainties of the mass bound have also been taken into account.
Although the neutrino signals from a galactic core-collapse SN explosion depend very much on the intrinsic properties of the progenitor star, such as the distance and the initial mass, the rapidly rising or falling feature in the time structure of SN neutrinos can be used to extract useful information about neutrino masses. For instance, the early-time neutrino burst from neutronization and the abrupt halt of neutrino signals due to the black hole formation will be advantageous for this purpose. We hope to return to those possibilities in the near future.