On the Possibility to Determine Neutrino Mass Hierarchy via Supernova Neutrinos with Short-Time Characteristics

In this paper, we investigate whether it is possible to determine the neutrino mass hierarchy via a high-statistics and real-time observation of supernova neutrinos with short-time characteristics. The essential idea is to utilize distinct times-of-flight for different neutrino mass eigenstates from a core-collapse supernova to the Earth, which may significantly change the time distribution of neutrino events in the future huge water-Cherenkov and liquid-scintillator detectors. For illustration, we consider two different scenarios. The first case is the neutronization burst of $\nu^{}_e$ emitted in the first tens of milliseconds of a core-collapse supernova, while the second case is the black hole formation during the accretion phase for which neutrino signals are expected to be abruptly terminated. In the latter scenario, it turns out only when the supernova is at a distance of a few Mpc and the fiducial mass of the detector is at the level of gigaton, might we be able to discriminate between normal and inverted neutrino mass hierarchies. In the former scenario, the probability for such a discrimination is even less due to a poor statistics.

To determine the neutrino mass hierarchy (MH), many methods by implementing accelerator neutrinos [13], reactor neutrinos [14][15][16][17][18][19], atmospheric neutrinos [20,21] or supernova neutrinos [22][23][24][25][26][27][28][29][30] have been proposed. Among them, those methods using supernova (SN) neutrinos are particularly interesting because of the interplay between intrinsic properties of massive neutrinos and the long-sought mechanism of SN explosions. Historically, the detection of neutrinos from SN1987A in the Large Magellanic Cloud [31,32] has motivated a huge amount of theoretical works in both SN physics and neutrino physics. Due to these potentials and gains from past experience, all modern running and forthcoming neutrino observatories, such as IceCube [33], Hyper-Kamiokande [34], JUNO [35] and KM3NeT [36] have set SN neutrino detection as one of their main physics targets. Now that neutrinos are massive particles, it is quite evident that for two neutrino mass eigenstates with masses m i and m j respectively, the time difference ∆t ij ≡ t i −t j for travelling from the SN to the detector is [37] ∆t ij = 5.15 ms · ∆m 2 ij /eV 2 ( E /10 MeV) 2 · D 10 kpc , where ∆m 2 ij ≡ m 2 i − m 2 j denotes the neutrino mass-squared difference, E is the average neutrino energy and D is the distance of the SN. Note that the energies of two neutrino mass eigenstates in Eq. (1) have been assumed to be the same and represented by the average neutrino energy E ≫ m i to estimate the difference in the time of flight (ToF). In this work, we show that SN neutrino emission with a characteristic time smaller than the difference in the ToF for different neutrino mass eigenstates can be utilized to determine the neutrino MH, if the time resolution of the future SN neutrino detector is good enough as expected.
For illustration, we consider two distinct scenarios in which these features of SN neutrino signals could be satisfied.
The first scenario is the short burst of electron neutrinos in the early-time process of neutronization (i.e., e − + p → n + ν e ) during supernova explosion, which takes place right after the prompt shock wave forms at the surface of inner core and starts to disassociate heavy nuclei in the surroundings. The existence of such a neutronization burst has been found to be a robust feature of SN neutrino emission in all the numerical simulations of SN explosions using different equations of state, treatment of gravity, and numerical approaches for hydrodynamics [38][39][40][41][42]. Depending on the details of simulations, the time duration of ν e burst, characterized by its full width at half maximum of the luminosity and denoted by ∆t N , can range from 3 ms to 20 ms. A typical value of ∆t N = 5 ms, which is close to the median of all the possible values that we have surveyed [27,[43][44][45][46][47][48][49][50][51], will be adopted in the following discussions.
The second scenario is the abrupt termination of neutrino emission due to the formation of a black hole (BH) in the accretion phase of a core-collapse SN. The BH formation is expected to occur when the mass of the progenitor star is between about 25 and 40 solar masses [52] and the shock wave cannot manage to successfully propagate out of the heavy outer core even with neutrino heating. Although there will be no final SN explosion, this scenario does have the advantage of an even shorter characteristic time, i.e., ∆t B ≈ 2R/c ∼ 0.1 ms, where the radius of active region for neutrino emission R ≈ 10 km and the speed of light c ≈ 3 × 10 10 cm s −1 have been used. At this point, we should mention that the ToF of massive neutrinos from the BH forming SNe have been considered previously by Beacom et al. in Ref. [53], where an upper bound on absolute neutrino masses has been obtained.
Different from previous works, the present paper will concentrate on the impact of ToF on the time distribution of SN neutrino events in future detectors. In particular, we explore the distinct features in the cases of NH and IH. In Ref. [27], the rising time of SN neutrino event rate in IceCube has been implemented to discriminate between NH and IH. Other possibilities have recently been summarized in Ref. [29]. The analysis presented in our paper can be regarded an additional effort in the same direction. The remaining part is organized as follows. In Sec. II, the time distribution of SN ν e -burst events in the water-Cherenkov and liquid-scintillator detectors is predicted by taking account of the difference in ToF of different neutrino mass eigenstates, while Sec. III is devoted to the tail of neutrino events in the case of BH formation. Finally, we conclude in Sec. IV.

II. THE NEUTRONIZATION BURST
For the neutronization burst, the neutrino flavor state |ν e is produced and then propagating from the SN to the detector on the Earth. During the propagation of such a long distance (e.g., D = 51 kpc for SN1987A), the coherence will be lost and the narrow peak of |ν e will split into those of three neutrino mass eigenstates |ν i (for i = 1, 2, 3) due to the difference in arrival time according to Eq. (1). Because the separation of these peaks depends on the MH in a simple and definitive way, we can immediately recognize the MH once SN neutrinos are detected with high enough statistics.

A. General Remarks
First of all, we give some general remarks on the time separation of three possible peaks and its dependence on the SN distance D, neutrino MH and neutrino flavor conversions.
Since the time resolution of most modern neutrino detectors are at the O(10) ns level, which is accurate enough to reconstruct the peak with a width about ∆t N ≈ 5 ms, the only requirement to temporally resolve the peaks of different neutrino mass eigenstates is that the difference in arrival time is larger than the peak width, i.e., |∆t ij | ≥ ∆t N . Given the neutrino energy, this inequality can be translated into a lower bound on the SN distance where we define D Nij to be the minimal resolvable distance. Then, corresponding to ∆m 2 21 and |∆m 2 32 |, there will be two independent distances D N21 and D N32 , satisfying the condition D N21 ≫ D N32 due to ∆m 2 21 ≪ |∆m 2 32 | from neutrino oscillation data. Some discussions on the SN distance are in order. In the case of large distances D > D N21 , for both neutrino MHs, there will be three well separated peaks with different ToFs and event rates, each corresponding to one of neutrino mass eigenstates. For the NH, there should first come into the detector two peaks of |ν 1 and |ν 2 at time t 1 and t 2 , respectively, with a difference of ∆t 21 = t 2 − t 1 > 0. Then after a relatively longer time ∆t 32 , which according to Eq. (1) should be ∆m 2 32 /∆m 2 21 times longer than ∆t 21 , there appears the third peak, corresponding to |ν 3 . For the IH, the peak corresponding to |ν 3 will reach the detector earlier than the |ν 1 and |ν 2 peaks do, with the order of the latter two being the same as that in the NH. Therefore, even though the size of each of the time differences ∆t 21 and |∆t 32 | are the same for two neutrino MHs, their temporal order appearing in the detector leads to a clear and unique signature of the MH. Once this ordering is experimentally observed, the MH will be determined completely. In contrast, if the SN distances turn out to be small D < D N32 , then all three peaks become indistinguishable, and it will be impossible to deduce the neutrino MH from the appearing order of the peaks.
If the distance is lying in between D N32 and D N21 , the peak corresponding to the mass eigenstate |ν 3 will be discriminable from the peaks of |ν 2 and |ν 1 , while the latter two are non-separable and we will effectively see only two peaks. Therefore, the method utilizing the temporal order of the observed peaks does not work in this case. However, we can make use of another piece of information, namely, the magnitude of event rate R i at those peaks corresponding to each neutrino mass eigenstate |ν i (for i = 1, 2, 3). To study the event rates of these peaks at the detector, one needs to know not only the initial |ν e spectrum and the reaction for detection, but also how the |ν e spectrum evolves when neutrinos propagating from the production region to the surface of the SN.

B. Flavor Conversions
The short-time burst of |ν e is generated when the first strike of the outer materials onto the inner SN core is bounced back and the heavy nuclei are disintegrated into free nucleons by the energetic shock wave [54,55]. During the outward propagation of SN neutrinos from the dense region to the surface, it is generally expected that both the ordinary Mikheyev-Smirnov-Wolfenstein (MSW) matter effects [56,57] and neutrino self-induced collective oscillations [58][59][60][61] will be crucially important for neutrino flavor conversions. See, e.g., Refs. [54,62,63], for recent reviews on collective oscillations of SN neutrinos.
For the progenitor stars of more than 10 solar masses, they will normally develop heavy iron cores before collapse. In this case, it is natural to expect that during explosion, collective neutrino oscillations occur within a few hundred kilometers above SN core, whereas the MSW effects will come into play far away in the SN envelope. However, for the stars of 8 to 10 solar masses, which finally evolve to the O-Ne-Mg core-collapse supernovae, the matter density profile above the core is so steep that the MSW resonances could happen within the region of collective oscillations [64][65][66]. For supernovae that allow only MSW effect for the neutronization burst, from Ref. [22], we see that in the NH case the initially generated neutronization |ν e burst is equivalent to the heaviest mass eigenstate |ν m 3 in matter. Due to a relatively large θ 13 [9], the flavor conversion proceeds adiabatically when |ν e passes through the density regions of ∆m 2 32 -and ∆m 2 21 -driven resonances and finally becomes a mass eigenstate |ν 3 in vacuum. For the IH case, the initial |ν e burst is almost a mass eigenstate |ν m 2 in matter. Again this neutrino state will also traverse the entire density profile adiabatically and become the mass eigenstate |ν 2 in vacuum after emerging.
On the other hand, for supernovae that allow not only the MSW but also the collective oscillations, the initial neutronization |ν e burst usually evolves to all three mass eigenstates after emerging with different probabilities. Ref. [65] have given transition probabilities of |ν e → |ν i as a function of energy in its Fig. 2. We will use this P ei (E) for the neutronization |ν e burst in supernova allowing the collective oscillation.
Although the self-induced collective oscillation has been proposed, it remains unclear to what extent will this occur in a real SN environment due to the large uncertainties such as progenitor mass, neutrino luminosity and simulation details including dimensionality and multi-angular/single-angular technique in dealing with collective oscillation. Also noticing that the researches on this topic are rapidly advancing, therefore we will remain conservative in this work by considering two different cases for the conversion probability from initial |ν e to |ν i after emerging: • Case (A) -For the O-Ne-Mg core-collapse SNe, both MSW effects and neutrino selfinduced collective oscillations play an important role, so the initial neutronization |ν e burst usually evolves to all three mass eigenstates in vacuum but with different probabilities. As already demonstrated in Ref. [65], the transition probabilities for |ν e → |ν i (for i = 1, 2, 3) are actually functions of the neutrino energy, i.e., P ei (E), where the spectral splits are found and explained analytically. In the following calculations, we use the probabilities from Ref. [65] as the first example.
• Case (B) -Second, as a simple working assumption, we neglect the energy dependence and specify the transition probabilities for |ν e → |ν i as P e1 = 1/6, P e2 = 1/3 and P e3 = 1/2, representing a class of scenarios in which the probabilities are comparable in magnitude. Although the exact values of those probabilities are not important, the discriminating power for MH will be lost if any one of them becomes negligibly small.
It is worthwhile to mention that there are large uncertainties in the progenitor mass, neutrino luminosity and the details of numerical simulations, such as dimensionality, equations of state and neutrino transport, so it is a complicated situation to deal with collective oscillations. In the following, we will compute the event rates of three neutrino mass eigenstates only in the above two cases for illustration.

C. Neutrino Event Rates
To numerically check whether the neutronization burst method is feasible in currently running and future detectors, we calculate neutrino event rates. Starting with the neutrino spectrum n ν e , we will adopt the quasi-thermal spectrum approximated by a Gamma distribution [67][68][69] n ν e (t, where the spectral index α is given by E(t) is the average energy, and E rms (t) is the root-mean-square, while the luminosity L(t) is given by the SN simulation data [69]. Since the burst of |ν e will finally evolve into an incoherent superposition of three neutrino mass eigenstates |ν i at the SN surface, the number densities of the latter are given by where we have dropped a common ToF of |ν 1 and identified the detection time t d as the the emission time of |ν 1 without loss of any generality.
For SN neutrinos of energies about a few tens of MeV, both elastic neutrino-proton and neutrino-electron scatterings can be observed in the liquid-scintillator detectors due to a low energy threshold, while only the elastic neutrino-electron scattering is observable in the water-Cherenkov detectors. Moreover, for liquid-scintillator detector, it allows chargecurrent interactions and neutral-current interactions with its 12 C nuclei. The differential cross section for elastic neutrino-proton scattering is universal for all neutrino flavors at the lowest order and given by [70,71] dσ where N p is the total number of protons in the target. For the liquid-scintillator detectors, the recoil energy of the final-state proton will be quenched significantly and one can establish the relationship between the original recoil energy and the observed energy as done in Refs. [71][72][73]. Furthermore, when the observed energy falls below 0.2 MeV, the radioactive backgrounds will be dominant. Therefore, we have placed an energy cut at T c = 0.2 MeV, corresponding to an original recoil energy of T p = 1.0 MeV, for which a minimal neutrino As we have mentioned, the recoil energy of proton under discussions is at most a few MeV and thus there will be no Cherenkov light, so it is impossible to observe any signals from neutrino-proton scattering in the water-Cherenkov detectors. For the elastic neutrinoelectron scattering, the total cross sections for ν e and ν x , where the latter collectively denotes ν µ and ν τ and their antiparticles, are well known in the standard model. At the tree level, the explicit expressions are [74,75] where the kinetic energy of the final-state electron T e ≡ E ′ e − m e is lying below T max For electron neutrinos, the coefficients are given by ǫ − = −1/2 − sin 2 θ W and ǫ + = − sin 2 θ W ; while ǫ − = 1/2−sin 2 θ W and ǫ + = − sin 2 θ W for muon and tau neutrinos.
Thus, the event spectrum of neutrino mass eigenstate |ν i can be simply written as where N e is the total number of electrons in the target and |U βi | 2 is the probability for the projection of |ν i to |ν β with U βi being the lepton flavor mixing matrix. For the neutrino reaction with 12 C, the charge-current interactions and neutral current interactions have been well established both theoretically and experimentally. In this work, for liquidscintillator neutrino detectors, we will consider these reaction too. We will directly use the cross-section tabulated in Tab. 1 of Ref. [76]. Denoting these cross-section as respectively, we can similarly compute the event spectrum of neutrino mass eigenstate |ν i caused by these four reactions with 12 C as where the subscript X stands for the four reactions in (12) and (13).
Further integrating the event spectrum (9) over the neutrino energy, we obtain the total event rate of elastic neutrino-proton scattering and likewise R i,ES for the neutrino-electron scattering (11) and neutrino-12 C reaction (15).
In the case of elastic neutrino-electron scattering, the observed recoil energy should also be larger than 0.2 MeV for the scintillator detector, but 3.5 MeV for the water-Cherenkov detector. This implies that the minimal neutrino energy E min ν will be different for these two types of detectors. If the target is composed of water, the fiducial mass of 2.5 megaton corresponds to N p ≈ 2 × 10 35 and N e ≈ 1 × 10 36 . For the scintillator detector of the same fiducial mass, the proton and electron numbers are quite similar. The number of 12 C, assuming that the liquid scintillator is chosen as linear alkyl bencene (C 18 H 30 ) as in Juno [35], can be calculated as N12 C ≃ 3N e /23.
In Fig. 1, we present the numerical results of the individual event rates R i,PS (i = 1, 2, 3) for three neutrino mass eigenstates and also the total rate (the thick and black curve) for either NH or IH. Similarly, the numerical results of R i,ES are shown in Fig. 2 and that of Fig. 3. For all these reaction channels, three representative distances, D 1 = 1 Mpc for the small distance, D 2 = 20 Mpc for the intermediate distance and D 3 = 400 Mpc for the long distance, have been considered. In addition, the transition probabilities P ei (E) given in Fig. 2 of Ref. [65] have been adopted in the calculations. Some comments on the results are in order.
First of all, let us recapitulate the fractions of neutrino mass eigenstates after the action of both MSW matter effects and collective oscillations on an initial flux of pure |ν e in a SN model from Ref. [65]. The fractions depend crucially on the neutrino energy and their main features can be summarized as follows: (1)  why |ν 1 dominates the contributions to the total event rate. As an immediate consequence of this observation, there will never appear two or three peaks, which is evident from all the plots in Fig. 1. Therefore, it is difficult to tell the difference between NH and IH for the small and intermediate distances. However, for the large distance D 3 = 400 Mpc > D N21 , the |ν 1 peak in the IH case is broadened significantly compared to that in the NH case.
The reason is simply that the lightest neutrino mass has been set to be vanishing in the numerical calculations, namely, the absolute mass m 1 ≈ 49 meV of |ν 1 in the IH case is much larger than that m 1 = 0 in the NH case.
Now we turn to the event rates of elastic neutrino-electron scattering depicted in Fig. 2.
Since the recoil energy of the final-state electron is not quenched in the liquid scintillator, all the neutrino mass eigenstates can contribute to the event rates. However, the |ν 3 peaks in the IH case are highly suppressed due to a tiny transition probability P e3 . Some discussions about this reaction channel are helpful.
• As neutrinos in the entire energy range contribute, the magnitudes of the event rates of neutrino-electron scattering in all cases of MH and distances are comparable to or cross section of neutrino-electron scattering in Eq. (10) compared to that in Eq. (8).
In addition, the cross section for electron neutrinos is about six times larger than that of other neutrino flavors, so the neutrino mass eigenstate that has the largest component of |ν e is most important. This non-universality of the cross sections results in different total event rates for the two MHs. Such a difference is most apparently seen by comparing the heights of the total rate peaks in the two subplots in the small distance case of D 1 = 1 Mpc < D N32 . However, the SN distance is too short for three peaks to be well separated in the NH case. Indeed, even if the probabilities have to be changed, as long as they are known to a good accuracy such that the heights and shapes of the |ν 3 peak and |ν 1 + |ν 2 peak for the NH and IH can be calculated, the MH can always be deduced by comparing them with the observed ones.
• For the large distance (D 3 = 400 Mpc > D N21 ), one can observe that in the NH case, the peaks due to |ν 1 and |ν 2 are already separated while the |ν 3 peak (not shown in order to make the |ν 1 and |ν 2 peaks clear) lags behind them by about 0.4 second.
In the IH case, the |ν 1 peak is also separated from the |ν 2 peak although the latter is now flattened too much to be seen as a peak. Therefore, the difference between the NH and IH event rates is also quite obvious to determine the MH in the large distance case. Compared to the intermediate distance case, however, the large distance case suffers the drawback of a much lower total event rate. For the neutrino reactions with 12 C that are allowed by liquid-scintillator detectors, the combined event rates of all four reactions in (12)-(13) are shown in Fig. 3 for three typical distances. It is seen by comparing with the event rates of proton scattering in Fig. 1 that generally, the shape of the event rates of the combined neutrino-12 C reaction is quite similar to that of the neutrino-proton scattering, although the magnitude of the former is only ∼ 1/3− ∼ 1/2 of the latter. This is understandable because the cross-section of neutrinoproton scattering is about 1.2∼1.5 times the total cross-section of neutrino-12 C reactions, and the proton density is 1.5 times the 12 C density. Moreover, while the neutrino-proton cross-section is universal to all neutrino flavors, the neutrino-12 C cross-sections (in particular the charge-current interaction (12)) are not. Also similar to the case of proton scattering, the main contribution from the total event rate is also from |ν 1 for the same reason as in Fig. 1. Moreover, this also implies that there will not appear two or three peaks in the event rate of the neutrino-12 C reaction channel and therefore difficult to tell the neutrino mass ordering.

D. Further Discussions
Finally, we discuss how the variations of a few important parameters affect the event rates. These key parameters include the transition probabilities P ei (E), the initial neutrino spectrum n ν e (t, E) and the absolute mass of the lightest neutrino. Here we give some brief remarks on the impact of those parameters and leave a full analysis for future works.

Transition Probabilities
For comparison, we adopt the transition probabilities in Case (B) of Sec. II B and present the event rates in Fig. 4, where the other input parameters are taken to be exactly the same as in Fig. 2. The main purpose for such a comparison is to show that the feasibility of our method does not depend too critically on these conversion probabilities. For the small distance, the total event rates for both NH and IH are almost identical and thus cannot be used to pin down the MH, as expected. Now the difference is that the probability P e3 in the IH case is not small anymore and therefore the |ν 3 peaks are no longer suppressed. For the intermediate distance, the shapes and arriving order of the |ν 1 + |ν 2 peak and the |ν 3 peak are clearly seen different for the two MHs, and therefore can also be used to tell the MH. For the large distance case, since this time in the IH case the |ν 3 is not suppressed by the probability, the |ν 3 peak is behind and ahead of the |ν 1 peak by about 0.4 second for the NH and IH, respectively.

Initial Spectrum
The initial spectrum Eq. (4) is expected to describe the true neutrino spectrum quite well [69]. Substituting it into Eq. (11) we can obtain the event spectrum and then study the impact of two pivotal parameters, namely, E(t) and E rms (t). If the average energy is increased, while other quantities are kept unchanged, the central energy of spectrum will shift to a larger value while its peak value will decrease due to the constraint from the fixed total luminosity. However, considering the energy dependence of the cross section, we find that a higher average energy will in general enhance the total event rate R i .
The parameter E rms (t) controls how much the spectrum is concentrated on the central energy. Its influence on both the neutrino spectrum and the event rate in Eq. (11) are secondary compared to E(t) . It turns out that as E rms (t) decreases, the total event rate R i also drops down but slowly. Besides its impact on the total event rate, a better concentration of the neutrinos in energy also means a narrower spread in the arrival time after including the time delay effect. Given that the main difficulty of our method here comes from the resolution of different peaks, a narrower distribution of the initial neutrino spectrum in energy should be more favorable.
The luminosity L(t) affects the resolvability of the peaks of neutrino mass eigenstate through its appearance in the spectrum in Eq. (4). As we have seen earlier, if the neutronization peak duration is prolonged by a factor of N, then we need the SN distance to be increased roughly by the same factor in order to achieve the same temporal separation of the peaks. Consequently, the total event rate will be reduced by a factor of N 2 . As demonstrated in Ref. [45,47], a more complete inclusion of relevant weak interactions and the full general relativistic treatment tend to increase the width of the neutronization burst. On the other hand, Ref. [27,49,50] suggest that the SN for a smaller progenitor mass (e.g., around 10 to 15 solar masses) will do the opposite, i.e., producing a neutronization neutrino burst narrower in time. The latter observation indicates that the SNe with smaller progenitor masses are favored for our method to work effectively.

Absolute Neutrino Masses
In the numerical results presented in Figs. 1, 2 and 4, we have assumed the lightest neutrino to be massless. Therefore, the shape of the event spectrum of the lightest neutrino resembles exactly the feature of neutronization burst from the numerical simulation data [69].
A nonzero mass of the lightest neutrino mass eigenstate leads to a time delay, and accordingly the luminosity peak will be widened in time and its maximum will be lowered, since the total neutrino number is kept unchanged. On the other hand, it is now known from the observational data of Baryon Acoustic Oscillation and Planck [77] that the sum of three neutrino masses has an upper limit, namely, m 1 + m 2 + m 3 < 0.176 eV, which after taking into account the measured neutrino mass-squared differences leads to an upper bound of m 1 < 0.052 eV for NH and m 3 < 0.044 eV for IH.
We have checked that even if the upper bounds on the lightest neutrino are saturated and the largest distance of D 3 = 400 Mpc is assumed, the width and height of the luminosity peak are changed only by about a factor of 1.5. There exists even more stringent constraints on the sum of three neutrino [78], for which we find that the luminosity peak is practically unchanged in both its width and strength compared to those in Figs. 1, 2 and 4. Moreover, the essential idea is to probe the relative shift in time among different mass eigenstates, which are little affected by the total peak shape and strength. Therefore, the impact of absolute neutrino masses on our method can be safely ignored.

III. BLACK HOLE FORMATION
The method using neutronization neutrino burst to determine the MH is not useful when the SN distance D is not large enough to resolve three peaks of neutrino mass eigenstates.
One can observe from Eq. (2) that D is mainly limited by the relatively large duration ∆t N ≈ 5 ms of the neutronization burst. Fortunately, in the scenario of failed SNe, there exists another characteristic process with an even shorter time span: the termination of neutrino signals due to the BH formation during the accretion phase.
The rate of BH formation core-collapse SNe and the formation mechanism are still uncertain, and numerical simulations crucially depend on the initial progenitor mass and the details of models. In many concrete SN models with BH formation, the neutrino signal will be abruptly terminated when the neutrino flux is still measurably high. After the energy-dependent ToF is taken into account, the sharp cutoff on the neutrino flux will cause characteristic signals at the detector with different descending rates, corresponding to different neutrino mass eigenstates. More importantly, the BH formation is a phase transition process that takes very short time ∆t B , which can be estimated as ∆t B 2R/c ≃ 0.1 ms, where R is the radius of the active region of neutrino emission. In comparison with ∆t N for the neutronization burst, ∆t B is much shorter and therefore allows for more practical SN distances and higher statistics.
For the resolution for the cutoff edges in the signals of three mass eigenstates, a similar criterion to Eq. (2) yields the requirement on the SN distance Corresponding to ∆m 2 21 and |∆m 2 32 |, there will also be two characteristic distances D B21 and D B32 , which take the typical values D B21 ≈ 2.58 Mpc and D B32 ≈ 79.3 kpc for an average neutrino energy E = 10 MeV and ∆t B ≃ 0.1 ms. Apparently, because ∆t B is about 50 times smaller than ∆t N , the minimal distances for given ∆m 2 ij and E to resolve the neutrino tail due to the BH formation should be about 2% of that for the neutronization burst. As a consequence, when the distance is the same, there will be a much larger event rate and statistical significance for the MH discrimination.
For clarity, we define the time derivative of the event rate as Γ i ≡ −dR i /dt, for which the time t max i to reach its maximum for the neutrino mass eigenstate |ν i will be different from one another. For the BH forming SN at a distance larger than D B21 , we can find that t max i for the three mass eigenstates will be temporally separated by a duration larger than ∆t B .
For NH, t max 1 and t max 2 of |ν 1 and |ν 2 appear temporally close to each other, but earlier than that of |ν 3 . For IH, the opposite is true, namely, t max 3 appears first and is further separated will not be resolvable but are still separated from t max 3 . If the relative strength of the signals for three mass eigenstates is known, then the MH can still be deduced from the shapes of the descending event rates. Similar to the case of neutronization neutrino burst, for D < D B32 , again the MH will not be deducible in this way. Taking these distances for example, we perform numerical computations of the event rates in the scenario of BH formation.
Just before the BH formation, all the neutrinos and antineutrinos of three flavors can be produced. For the detection of antineutrinos, it is obvious that we should first consider the inverse beta decay (IBD) of ν e , whose cross section for SN neutrinos is much larger than those of other reactions in water or liquid scintillator. For neutrinos, the main contributions are from elastic scattering on protons and electrons. Therefore, we will only take account of these three processes in the our calculations. The IBD event rate of the i-th mass eigenstates of antineutrinos |ν i is given by while the total rate can be calculated by summing over all possible contributions, namely, where the part of elastic proton scattering should be omitted for a water-Cherenkov detector. A simple approximation to the IBD cross section was given in Ref. [79]. However here we will use a more accurate differential cross-section given by Ref. [80] which takes account into the weak magnetism effect and nucleon recoil effect dσ IBD d cos θ where θ is the angle between the antineutrino and positron direction in the lab frame, σ 0 , f, g, f 2 are some constants, E  (19) over θ allows us to find the total IBD crosssection σ IBD (E ν ), which we will direct use in this work. It is worth mentioning that R i,PS and R i,ES in this section include the contributions from both neutrinos and antineutrinos, for which the cross sections of neutrino-proton scattering in Eq. (9) are equal while those for neutrino-electron scattering in Eq. (10) are different for electron and non-electron flavors. Furthermore, the cross section in Eq. (10) can be adapted for antineutrinos by exchanging The fluxes F ν i and F ν i at the time t d on the Earth can be computed from the spectra n νe (t, E), n ν e (t, E) and n ν x (t, E) where ν x collectively denote ν µ and ν τ and their anitparticles, multiplying them by the conversion probabilities P βi = |U βi | 2 and taking into account the time delay due to Eq. (1). Explicitly, the formula for |ν i can be expressed as which is very similar to Eq. (7) except that the contributions from all flavors are summed up and the conversion probabilities P βi (for β = e, µ, τ ) are energy independent. Sucn an independency is expected because the BH formation normally takes place in the accretion phase, when the self-induced collective oscillations are found to be suppressed by the large matter density [81][82][83]. Hence the neutrino spectra experience only adiabatic conversions inside the SN, and then the mass eigenstates freely stream from the SN surface to the Earth.
In our numerical calculations, we use the n ν e , n ν e and n ν x spectra obtained from the SN simulation in Ref. [84] for the progenitor star of 30 solar masses. After converting them into the spectra of neutrino mass eigenstates, we find the energy integrated spectra for each ν i and similarly for each ν i the corresponding N ν i (t), which together with their corresponding average energies are given respectively in the left and right panels of Fig. 5. Since n ν e is slightly larger than n ν e , we have the integrated spectra N ν i also larger than N ν i before BH formation. Moreover, the time structure of neutrino signals is mainly determined by the neutrino emission in a very short period right before the BH formation, we can verify that the number spectra are constant in this short time slot before BH forms [53]. Some observations from Fig. 5 are summarized below: 1. The integrated spectra before BH formation decreases in the order of |ν 1 , |ν 2 and |ν 3 , and that of |ν 1 is roughly 1.5 times that of |ν 3 . In addition, the integrated spectra for neutrinos |ν i are almost in the same order as those of the corresponding antineutrino |ν i , except that the |ν 1 spectrum is slightly larger than that of |ν 1 due to the higher initial luminosity of |ν e than |ν e . In Fig. 6 we show the event rates around the time of the termination of the neutrino signal due to BH formation, where the SN distance is D = 5 Mpc and a 2.5 megaton liquidscintillator detector is assumed. The results for NH and IH are presented in the left and right column, respectively. Some comments on the numerical results are in order.
• From the first row of Fig. 6, we present the IBD events for both hierarchies, where one can see that the |ν 1 rate roughly doubles that of |ν 2 which is much higher than the In contrast, the average neutrino energy of the neutronization burst is only about 12 MeV.
• The total event rates, and the its time derivatives for each mass eigenstate are given respectively in the last two rows of Fig. 6 for both NH and IH. We can see that the total rates at the time long before and after the BH formation are equal for both hierarchies. It should be noticed that the signal decays for heavy mass eigenstates (i.e., |ν 1 and |ν 2 and their antiparticle states for IH while |ν 3 and |ν 3 for NH) are much slower than those for light mass eigenstates. This is the key feature and can be used to distinguish one neutrino MH from another if enough neutrino events can be detected. Such a difference between NH and IH can also be clearly seen from the time derivative of event rates shown in the last row. For the NH, the peak of |ν 3 + |ν 3 is low in height and slowly decreasing, and comes later than the higher peaks of the first and second mass eigenstates. For the IH, the opposite happens: the |ν 3 + |ν 3 peak becomes sharp and decreases fast while located on the right of the peaks of other mass eigenstates.
In the above numerical computations, the SN distance D is set to 5 Mpc, which is larger than D B32 ≃ 1 Mpc, the distance required for the maximal descending rates of |ν 3 and |ν 2 to separate by a gap larger than ∆t B , obtained by using E = 24 MeV, as shown in Fig. 5. With the total event rate in Fig. 6

B. Further Discussions
Now we discuss how the input parameters affect the tail of neutrino signals from the BH forming SN. Compared to the case of neutronization burst, the situation for the BH case is much simpler. First, the variation of neutrino spectra with time is practically constant because the time window under discussion is narrow enough. Second, the flavor conversions of neutrinos and antineutrinos are adiabatic, and thus it is straightforward to figure out the transition probabilities. Third, the impact of the absolute mass of the lightest mass eigenstate, as argued in the case of neutronization burst, is expected to be small if the cosmological bound on the sum of neutrino masses is applied. Therefore, we analyze the effect of the distance D on the total event rate in the following.
The first and most important effect of the SN distance is that the total event rate during the decay of neutrino signals is inversely proportional to D 2 . Moreover, the distance also has an influence on the shape of the total event rate, which turns out to be crucial for the MH discrimination. In Fig. 7, we plot the total event rates as functions of time for a series of distances ranging from 200 kpc to 10 Mpc for the two MHs. In fact, it is more practical to determine the distance by using other methods, e.g., optical observations, and then the MH can be determined by comparing the measurement with the predicted signal shape. For the previous example, once the distance is determined (e.g., D = 1 Mpc), then the neutrino MH can be deduced by a comparison of the descending rates dR/dt of both hierarchies.

IV. CONCLUDING REMARKS
In this work, we have carried out a phenomenological analysis of the ToF effects of massive neutrinos, and applied them to the SN neutrinos with short-time characteristics. For the neutronization ν e burst of core-collapse SNe, the peaks corresponding to different neutrino mass eigenstates for the NH will appear in a temporal order different from that for the IH.
A clear discrimination between two MHs requires a large distance for the SN, whereas a high statistics favors SNe at small distances. For this reason, we have found that it seems impossible to determine the neutrino MH via the neutronization burst for a typical corecollapse SN in the currently running and near future SN neutrino detectors. However, for the BH forming SNe, the abrupt termination of neutrino emission is shown to allow for a determination of neutrino MH with an enough statistics. In this case, a liquid-scintillator detector of one gigaton will register about 60 events for a SN located at 1 Mpc during the period of BH formation.
The impact of neutrino flavor conversions (particularly those in the case of neutronization burst), absolute neutrino masses and the supernova distance (particularly in the case of BH forming SNe) on the discriminating power is also discussed. It is shown that neutrino flavor conversions generally will not destroy the feasibility of the basic idea as long as the flavor conversions do not completely suppress the peaks in such a way that only one mass eigenstate is left. As for the absolute mass of the lightest neutrino, it has been numerically verified that if the cosmological upper bound on the sum of neutrino masses is applied, the shapes of the neutrino event rates in both neutrinonization burst and BH formation cases will have little changes. The SN distance does affect the shape of the neutrino signal decay in the case of the BH forming SN. In reality, one can fix the SN distance by other means and then extract the information of neutrino MH from observations.
There are a few other factors that are related to the applicability of our conclusions.
The most important one is the uncertainty in the luminosities and spectra of SN neutrinos.
For the neutronization burst, the uncertainties involved in the numerical simulations of SN explosions can introduce noticeable differences in the duration time of the burst and its luminosities. Moreover, the simulations do not take into account the full neutrinoneutrino interactions which severely affect the emerging neutrino spectra. If this burst is narrower, with a higher luminosity and more |ν 3 component, then the approach of using the neutronization burst to probe the MH becomes more promising. Secondly, the fraction of the BH forming SNe among all core-collapse SNe can vary from a few percent to a sizable value (see, e.g., Ref. [53]), and the neutrino luminosities and spectra before the BH formation could also be quite different [51]. The value that we have used for the luminosity of each flavor is on the order of 10 52 erg s −1 , which comes from a SN with BH formation at its early stage [53]. If the BH is formed at a later stage, then the luminosity and thus the total event rate can be one order of magnitude lower. One more factor associated with SN neutrino detection is about the low supernova probability within applicable distance. It was estimated that the core collapse SN rate will be only 3.2 +7.3 −2.6 in the Galaxy per century [85] and that within 1 Mpc is only about 0.04∼0.08 per year [86]. These numbers suggest that other experimental methods might potentially discriminate the neutrino mass hierarchies earlier than any new supernova neutrino observation.
Besides the neutronization burst and the BH termination of neutrinos, we emphasize that the basic idea of using the ToF difference to distinguish the MH is also applicable to any neutrino sources with short-time characteristics. Here we considered only two such characteristics, i.e., the neutronization neutrino |ν e burst and neutrino spectrum tail during black hole formation, because these are the only two theoretically known features with short time duration in SN neutrino spectrum. Given that there was effectively only one SN neutrino observation (SN1987A) and the large uncertainty in simulations, the true detailed SN neutrino spectrum with high statistics is still experimentally unknown. Therefore there remains possibility that some other features with similar or even shorter time duration might exist in the spectrum. If so, these features will also be usable to determine the neutrino MH in the same way. However, until then, to more accurately validate the idea in this work for SN neutrinos, we will have to perform more careful studies of the initial neutrino fluxes from realistic SN simulations including a full treatment of neutrino flavor conversions, which will be left for future works.