SN1987A neutrino burst: limits on flavor conversion

In this paper, we revisit the SN1987A neutrino data to see its constraints on flavor conversion. We are motivated by the fact that most works that analyze this data consider a specific conversion mechanism, such as the MSW (Mikheyev-Smirnov-Wolfenstein) effect, although flavor conversion is still an open question in supernovae due to the presence of neutrino-neutrino interactions. In our analysis, instead of considering a specific conversion mechanism, we let the electron antineutrino survival probability $P_{\overline{e}\overline{e}}$ be a free parameter. We fit the data from Kamiokande-II, Baksan, and IMB detected spectrum with two classes of models: time-integrated and time-dependent. For the time-integrated model, it is not possible to put limits above $1\sigma$ (68% confidence level) on the survival probability. The same happens for the time-dependent model when cooling is the only mechanism of antineutrino emission. However, for models considering an accretion phase, $P_{\overline{e}\overline{e}}\sim0$ is strongly rejected, showing a preference for the existence of an accretion component in the detected antineutrino flux, and a preference for normal mass ordering when only the MSW is present.

One of the main questions regarding supernova neutrinos today is the flavor conversion mechanism.It is expected for the supernova neutrinos to suffer MSW conversion [12][13][14] and a substantial number of works were done considering this as the only conversion mechanism in action, including the ones that analyze the SN1987A data [6,7].However, today it is expected that neutrino-neutrino interactions (forward scattering) become relevant in a supernova environment leading the neutrinos to a non-linear collective evolution [15].Due to the complications that emerge from this type of evolution, there is not a conclusive picture of neutrino conversion in the supernova environment.
Nevertheless, given the equal amount of non-electron antineutrinos ν x = (ν µ , ν τ ) emitted from the supernova, it is possible to write the flavor conversion in terms of only the electron antineutrino survival probability P ee .Therefore, we treat this probability as a free parameter to see how SN1987A data can constrain it.Something similar was done by F. Vissani in [16].However, it seems that the influence of the survival probability is analyzed only for the MSW normal hierarchy scenario (P ee = 0.64) against the no oscillation one (P ee = 0).Here we take a more complete analysis for P ee , allowing it to range from 0 to 1.
In section 2 we describe our model for the detected event rate in each detector (KII,IMB, Baksan) based on two different neutrino emission models, the flavor conversion mechanism, and the detection properties.In section 3 we describe our statistical analysis of the SN1987A data.In section 4 we show our results and discuss them, and finally, in section 5 we present our conclusions.

Model for the neutrino signal
In this section, we describe the model for the expected neutrino event rate in each of the detectors, which is used to fit the SN1987A data.First, we describe the two neutrino emission models considered in this paper: a time-dependent and a time-integrated.In sequence, we describe the flavor conversion in the flux, which depends only on P ee , and, in the end, we discuss the detection features of this analysis.Given that the most relevant cross-section for the considered detectors is the IBD, we will restrict our model to the antineutrino sector ( νe , νµ , ντ )
Time-dependent Given that the neutrino emission evolves in time, a time-dependent model should be at least considered in data analysis.This approach can be found in the famous paper of Lamb and Loredo [6] and some other works [7].In this approach, the antineutrino emission can be divided into two phases: the accretion and cooling phases.Here we will follow the path of [6,7] and model each phase by its most relevant mechanism of emission.
In this case, the accretion phase can be modeled as a positron thermal flux with temperature T a incident in a neutron target, that composes the mass in accretion in the protoneutron star.Therefore, as in [6,7], we consider that only electron antineutrinos are emitted in this phase and the flux is given by: with , where N n (t) is the number of neutrons as a function of the time, σ e + n (E ν ) the positron-neutron cross-section, and g e + (E e+ , T a ) the thermal distribution of positrons with energy E e+ in a temperature T a .The number of neutrons is given by the initial accreting mass M a with a fraction of neutrons Y n , and its time behavior is given by the factor j k (t) = exp − (t/τ a ) k , with τ a being the characteristic time of the accretion phase and the parameter k = 2 following the parametrization in [7] 1 .The denominator 1 + t/0.5s, as in [6,7], is used to mimic the behavior from supernova simulations, where we have a constant flux within the first 0.5 s followed by a fast decrease.The cooling phase, which is dominated by neutrinos and antineutrinos of all flavors emitted by the cooling neutron star, is modeled by a thermal distribution of fermions with temperature T c (t), with characteristic time τ c , emitted from a sphere with fixed radius R c and is given by with the cooling temperature being a function of time As already pointed out, different from the accretion component, the cooling one is composed of antineutrinos of all flavors.However, the non-electron antineutrinos ν x are emitted from deeper regions in the supernova, which can be effectively implemented by considering that they are emitted with higher initial temperatures T c, νx .In fact, during the rest of the paper, we will talk about the ratio between the flavors temperatures τ = T νx /T νe .
To combine the fluxes of both phases of emission, we follow [7] where the cooling phase starts after the accretion one.As argued in the cited work, if the accretion and cooling phases were contemporaneous the first seconds would be composed of two different spectra, given the different temperatures of each of these phases.As numerical simulations of supernovae do not show this feature, we assume that the different emission phases are separated in time.We do this using the following parameterization: where the accretion flux is only composed of electrons antineutrinos φ 0 a, νe , while the cooling flux contains an electronic φ 0 c, νe and non-electronic component φ 0 c, νx .

Time-integrated
In this model, we consider that the timeintegrated flux can be described by the following pinched spectrum [17]: where, for a specific neutrino flavor β , L β is the total energy (time-integrated luminosity), E 0β the mean energy, and α β the pinching parameter.We are mainly motivated to use this model due to a collection of works that only use the energy information from the SN1987A [8][9][10].Although the time data could bring new information, it is interesting to check if the energy alone can say something about the flavor conversion.

Flavor Conversion
From emission until detection, the neutrino may suffer flavor conversion.It is still an open question for supernova neutrinos which is the complete mechanism of flavor conversion, given the complications that arise with neutrinoneutrino interactions.However, due to unitarity and the equal initial flux of non-electron antineutrinos φ 0 ν µ = φ 0 ν τ = φ 0 ν x , the equations for flavor conversion can be simplified so that it will only depend on the electron antineutrino survival probability P ee and initial fluxes [18], such that Therefore, we can explore the survival probability P ee as a free parameter representing the flavor conversion occurring during the neutrino propagation.In this paper, we want to see how strong the SN1987A data can constrain P ee in the fitted models, given that the flavor conversion mechanism is still an open question in a supernova environment.Although this probability may be time and/or energy-dependent, we will consider it independent of these variables, given that we do not want to use a specific model.
We will also consider the MSW-only conversion scenario in order to compare it to our free P ee model.In this scenario, the electron antineutrino is created as a ν1 for normal mass hierarchy (NH) and ν3 for inverted mass hierarchy (IH).Therefore, the survival probability for each mass ordering can be written as follows [19] where we have considered an adiabatic evolution, with a flipping probability equal to zero at the high and low-density resonances.The vacuum mixing parameters are taken from the update values published for the global fit analysis in [20].Although this energy dependence of P ee is negligible in the standard MSW effect, other possible effects associated with collective effects, such as spectral split among different neutrino flavors lead to a strong energy dependency, changing drastically this scenario [15].However, given the unknowns associated with such collective effects nowadays, we limit our analysis to consider a P ee that is uniform in energy, leaving the spectral split analysis for a future work.

Detection
In the case of the SN1987A, we have data from three detectors: Kamiokande-II, IMB, and Baksan.In all of them, the dominant channel for electron antineutrino detection is the Inverse Beta-decay (IBD), which is the only one that we will consider.Therefore, the event rate R IBD νe as a function of the positron measured energy E e + , the angle between the incoming neutrino and the scattered positron θ and time (for the time-dependent model) can be calculated as follows where N p is the number of free protons, φ νe (E ν ,t) the electron antineutrino flux at the detector, dσ IBD νe (E ν )/d cos θ the differential cross-section for IBD, and η d (E e + ) the detector intrinsic efficiency.For the IBD, the incoming neutrino energy E ν is related to the created positron energy by E e + ≈ E ν − 1.293MeV , due to the mass difference between the initial proton and the final neutron.The energy threshold for the IBD is E th ν = 1.806MeV [21].

Efficiency
As pointed out by [16], when calculating the differential event rate in equation 9, one should use the detector intrinsic efficiency η d (E e + ).However, when integrating the event rate to get the total number of detected events, one should account for the threshold energy considered when selecting the events.This is achieved by multiplying the intrinsic efficiency by a function g(E e + , E min ) resulting in a total efficiency g(E e + , E min ) = in which the error function Erf accounts for the threshold energy E min and the uncertainty σ (E e + ) on the energy.This distinction between intrinsic and total efficiency is relevant when talking about the ones reported by the experiments, which are total efficiencies accounting for the threshold energies used during the events selections.This distinction becomes even more relevant in the case of the Kamiokande-II when using the low-energy events (numbers 13-16 nad 6 in table 5) added a posteriori and which are below the energy threshold of 7.5 MeV used in the first published data.To incorporate these events in our analysis, we need to infer the intrinsic efficiency from the published total efficiency and extrapolate the last to lower energies in the case of Kamiokande-II.Following this reasoning, we adopt the same parametrization for the intrinsic efficiency as reported in [16], with E min = 4.5 MeV for Kamiokande-II.Both total and intrinsic efficiencies used in this work are shown in figure 10.

Uncertainties
The uncertainties used in this work are experimental ones shown in tables 5, 6, and 7.Although we have the angle uncertainty, we will not consider it in our analysis, due to its non-significant impact on the likelihood, given that the considered cross-section (IBD) has a weak angular dependency.Also, as pointed out in [7], the relative time between the events is measured with good precision so that we also ignore the time uncertainty.As for the energy uncertainty, in addition to reported values for the energy of the events, to implement it in the efficiency expressions, such as equation in 10b, we need to estimate the uncertainty for other values of energy.For this purpose, we adopt an uncertainty parametrization with a statistical component that goes with the square root of the measured energy E e + and a systematic one that grows linear with the energy.as done in [16]: The values that we used for the coefficients are shown in table 4 corresponding to the ones that best adjust the function to the reported uncertainties.

Cross-section
The exclusive interaction considered in the analysis was the inverse beta decay, given the high cross-section compared to other possible channels of KII, IMB, and Baksan.We adopted the differential cross section (in the scattering angle) calculated by Vogel and Beacom in [22].

Off-set time
Another thing that we have to be careful of is to not confuse the time of the first detected neutrino t 1 with the time t 0 = t = 0 which indicates the time that the first neutrino arrives at the detector, even if it was not detected.Not considering this may force that the first detected neutrino is originated from the initial accretion phase, which may not be the case.As we will discuss later, for the MSW conversion in the inverted mass hierarchy scenario (IH), the initial νe flux contributes only to 2% of the detected flux, which makes it probable that the first detected neutrino came from the cooling phase and then t 1 ̸ = t 0 .To get around this problem, it is usual to introduce an offset time t d off = t 1 − t 0 between the first detected neutrino and the time of arrival of the first neutrino, which may be different for each detector given that they do not have an equal absolute time.

Background Modeling
In a realistic approach, we have to consider that detected events may come from background sources.The background rate is considered to be constant over the time of exposure, and also uniform over space, i.e., it depends only on the positron energy of the event B = B(E i ) = d2 N B /dtdE.The independence regarding the spatial position is an approximation, given that there is more background at the wall of the detector, due to the surrounding material.
The background can be measured and it is published by the collaborations.As argued in [23], there is no need to do a convolution of these measured background rates with a Gaussian uncertainty in the energy, as done in [6], given that the background curve adjusted to the data already accounts for the uncertainty in the measurement.Therefore, one only needs to take the background rate from the experimental curve without doing a posteriori uncertainty convolution, which would double count the uncertainty effect.In our case, we use the background rate from [16] for both Kamiokande-II and Baksan, whereas the background is irrelevant for the IMB detector.In the case of the Time-Integrated analysis, we have to integrate the background rate in time to get the event rate per energy B = B(E i ) = dN B /dE.The integration has to be done on the time of exposure to the supernova signal, i.e., the data-taking duration (∼ 30s).

Statistical Analysis
For the statistical analysis, we use the method of maximum unbinned likelihood, due to the low number of events.Our expression for the likelihood is similar to the one adopted in [7] Here we made implicitly the dependency of L in the parameters of our models.In this equation, i is the index of each event, R(t, E, cos θ ) is the expected event rate from equation (9), R(t) the event rate integrated in the angle and energy, and B the background rate 2 discussed in section 2.8.
Here we differ from [7] in the definition of R(t), in which we consider the total efficiency to calculate the event rate integrated in the energy, as discussed in section 2.4.The integration in the positron energy E e is made considering a Gaussian distribution L i (E e ) around the measured value E e,i with standard deviation given by the measurement uncertainty.As already discussed, we consider that the time and angle uncertainties are irrelevant.We also consider the dead time τ d for each detector (d = K, B, I), where f d is the live-time fraction [7].In the case of the time-independent model, we only have to consider a time integration in the event rate for the signal R(t i , E e,i , cos θ i ) and for the background B(E i ).
To find the set of parameters that best adjusts our model to the data, we only have to maximize the likelihood L or minimize −2 log(L ).The last one is useful because it transforms multiplication into a sum and has a straightforward connection to confidence intervals.Given that we have a set of parameters ⃗ θ , taking their the best-fit ⃗ θ we can define the likelihood ratio as follows.
so that −2 log λ ( ⃗ θ ) follows a χ 2 distribution in the asymptotic limit of large samples N → ∞, with m degrees of freedom representing the number of parameters not constrained to be in its best-fit value.With this procedure, we can estimate the best-fit values for the parameters and their confidence interval, given a confidence level.However, we have to note that our data is not a large sample so our confidence level is an approximation.In any case, in this paper, we consider that it is an acceptable approximation given the allowed region for the astrophysical parameters to be comparable to previous works [6] that use other approaches to set the confidence levels, as we discuss in Appendix A.

Time-dependent model
For the time-dependent model, following the references [6,7], we consider two possible cases, one with just cooling emission and the other with an initial accretion phase.For the cooling component, we have four astrophysical parameters, the initial cooling temperature T c , the time constant of the phase τ c , the radius of the neutrinosphere R c , and the ratio between the initial temperatures of the electronic and non-electronic antineutrinos τ = T νx /T νe .Previous works [7] fix this temperature ratio based on supernova simulations.Here, we check the impact of changing this ratio given that it has strong implications in how similar the initial spectra are, which reflects how well we can identify flavor conversion in the detected spectrum.Nevertheless, we limit ourselves to the range of temperature ratio expected from supernova simulations [17].When considering the accretion phase, we introduce three new astrophysical parameters: the initial accretion temperature T a , the time constant of the phase τ a , and the accretion mass M a .In addition to the astrophysical parameters, there is the offset time for each detector and the survival probability, resulting in a total of 8 parameters for the cooling model and 11 for the cooling plus accretion.
To analyze how the SN1987A data can put limits on P ee , we can do a marginal analysis, as described in section 3. Figures 1 and 2 show the marginal plot of P ee for the models with only cooling component and for the one with cooling and accretion, respectively.For the model with just cooling, we can see that it is not possible to put limits on P ee up to the 1σ for τ values considered.This probably happens because both initial fluxes φ 0 ν e and φ 0 ν x come from the same mechanism, resulting in almost indistinguishable spectra, even allowing the temperatures to be different.
When we consider the accretion phase, we have a different scenario, where P ee ∼ 0 is strongly rejected, as we can see in Figure 2.This stronger constraint in P ee happens because in the accretion mechanism only electrons antineutrinos are emitted, making their initial flux φ 0 ν e more distinguishable from the non-electronic one φ 0 ν x , which in turns facilitates the identification of flavor conversion.Given that, the excluded region of P ee ∼ 0 corresponds to the case where the detected flux is composed only by the initial φ 0 ν x , i.e., a flux with no accretion component.This shows us that the detected electron antineutrinos are better described by a flux with an accretion component coming from φ 0 ν e , as already found by [6].However, in [6] they do not consider the role of flavor conversion, while here we can see that the existence of an accretion component has strong implications on the conversion mechanism.If we consider only the MSW effect with adiabatic propagation, this implies that the normal hierarchy scenario is favored over the inverted.Comparing them with the best-fit of free P ee , the normal hierarchy scenario is not significantly rejected, while the inverted one is rejected by ∼ 3σ of significance.
It is also possible to see in figure 2 some kind of discrete transition to a lower ∆ χ 2 at P ee ∼ 0.5.This happens because there is a preference for a non-zero off-set time in the IMB data, as can be seen in the best-fit value of t I o f f in table 2 if the accretion component is strong enough (MSW-NH or free P ee ).However, if we go to lower values of P ee , such as in the MSW-IH, it becomes preferable to describe some of the first events of IMB as coming from the cooling, i.e. t I o f f = 0.This transition can be seen in figure 3 in which we plot the ∆ χ 2 profile for t I o f f = 0.5 and 0 s.We have also tested the implications of considering the cooling and accretion components as contemporaneous.As argued by [7], there is no evidence of a composed spectrum in supernova simulations, so the two mechanisms with dif- ferent mean energies should occur at different times.However, from supernovae physics, we may expect that the PNS starts to cool down by neutrino emission soon after its formation, simultaneously with the accretion mechanism [24].Therefore, we decide to test the implications of that hypothesis in our analysis.As we can see in Figure 4 there is no significant modification on P ee limits.The only modification appears on the best-fit of t IMB off , which can be seen in Appendix A.

Time-integrated model
For the time-integrated model, we considered a Fermi-Dirac emission (α ν e = α ν x = 2.3), a choice that does not have big impact in the fitting for 2.3 < α < 4 3 .We also con- 3 By letting α νe and α νx run free in this interval, the variation of the likelihood ratio L /L max was not above 1σ (C.L. ≈ 68%).sider a hierarchy for the mean energy E ν x > E ν e , which is physically motivated given that non-electron neutrinos interact less (lack of τ and µ leptons in the environment) and then escape from deeper regions in the supernova with higher temperatures.The best-fit values for the astrophysical parameters are shown in Table 3 considering the 3 different conversion scenarios.As we can see, there is a preference for a detected spectrum φ ν e to be composed mostly by the initial non-electron neutrino spectrum φ 0 ν x , given that there is basically no constraint for the total energy ε ν e , the same behavior was also found in [10].Even in the MSW mechanism with inverted mass hierarchy, where the composition of φ 0 ν x in the final flux is small (P ee ≈ 2.18%, the flavor conversion is compensated by a higher total energy ε ν x .This preference is a combination of the imposed energy hierarchy E ν x > E ν e and the low detection efficiency for lower energies, where the low energy events can be as well described as coming from the background.However, we did not inves- tigate this preference deeply 4 .As we are interested in the flavor conversion parameter P ee , we leave the Appendix A to compare our marginal and contour plots with previous analyses to show the consistency of our method, at least regarding the astrophysical parameters. For the flavor conversion analysis, we again fix the initial temperature ratio (more precisely the mean energy ratio τ = E ν x /E ν e = T ν x /T ν e ) and let the other parameters run freely over the allowed range (Table 3).Figure 5 shows the marginal plot of P ee minimizing over the other model parameters.Again, there is no constraint on the survival probability above 68% of confidence, even for spectra with higher mean energy differences such as τ = 1.4.

Problems with fitting the data with some models
In our numerical implementation, we found some difficulties in working with the two-component model (accretion + cooling).The main one is the existence of different local minima, which make the minimizer algorithm give different best fits depending on the initial conditions.To get around this problem, we used two methods to find the global minimum.In the first method we fit this model multiple times (≈ 1000) fluctuating the initial conditions of parameters uniformly in the ranges shown in Table 2, and taking the minimum value of −2 log L as the initial condition to find the global best-fit.The second method was based on using different minimizers (MINOS, scipy, simplex) 5 to see if this dependency on the initial conditions was algorithm dependent.In the end, we found that all the different minimizers obtained the same best fit given initial conditions around it, and in agreement with the first method.Given the concordance between the two methods and algorithms, we have confidence that the best fit obtained is the most probable one inside the allowed parameter space.

Conclusion
In this paper, we have explored the role of flavor conversion in the SN1987A neutrino data, and how it can impose limits on the flavor conversion mechanism.We found that the time-integrated model, which uses only the energy information, could not put any limit on the electron antineutrino survival probability P ee .The same happens for the timedependent models that consider antineutrino emission only from the cooling mechanism.However, with the existence of an accretion emission of electron antineutrinos, strong limits are imposed on low values of P ee .This is impressive given the low statistics of the SN1987A neutrino data and it is in agreement with the previous work of Lamb and Loredo [6] in which the data shows a strong preference for the existence of an accretion component.
In previous works, such as [19], it was already pointed out that the inverted mass hierarchy was disfavored in MSW adiabatic scenario with a significance of 3σ for some values of θ 13 , which was unknown at that time.Here we confirm this statement, as it can be seen from the figures 2 and 4. Our improvement to their analysis was to use the current well-known neutrino vacuum mixing angles [20] and extend the analysis to the whole spectrum of possible values for the survival probability P ee .
As we discussed, our analysis does not consider any time or energy dependency on P ee , which may happen when we consider collective effects due to neutrino-neutrino forward scattering.We leave the study of time and energy dependency for a future paper.In any case, our results can still be used to constrain conversion models that result in a fixed value for P ee .Fig. 6 T c,0 vs R c contour plots comparing our results with previous ones [6,7].
parameter agrees with the obtained contour plots.Also, the conversion model with free P ee encompasses the MSW-IH and MSW-NH scenarios, as one would expect given that the latter are specific cases from the former, with P ee ≈ 2.18% and P ee ≈ 67.8% respectively.7 Marginal and contour plots for the astrophysical parameters T c , τ c , R c for the only cooling model, keeping the detection off-set times t KII off ,t IMB off ,t Bak off in their best-fit value.For the contour plots, we use color bands for the MSW-IH and MSW-NH scenarios and lines for the free P ē ē, corresponding to confidence levels of 68% (dashed) and 99.7% (solid).Note that minimum χ 2 min = −2 log L max is the absolute one among all the curves/conversion scenarios.

2 |U e1 | 2 |U e3 | 2 MSWτFig. 1 P 2 τFig. 2
Fig. 1 P ee likelihood ratio (∆ χ 2 = −2 log L /L max ) for the SN1987A data considering the time-dependent model with only the cooling component.The horizontal dashed lines correspond to 1, 2 and 3σ of C.L. Note that minimum χ 2 min = −2 log L max is the one absolute regarding all the curves.

Fig. 3 2 τFig. 4
Fig. 3 Same as Fig. 2 but fixing τ 1.2 for two different values of off-set time for the IMB data t I o f f

2 τFig. 5 P
Fig. 5 P ee likelihood ratio for the SN1987A data considering the timeintegrated model.

Fig.
Fig.7Marginal and contour plots for the astrophysical parameters T c , τ c , R c for the only cooling model, keeping the detection off-set times t KII off ,t IMB off ,t Bak off in their best-fit value.For the contour plots, we use color bands for the MSW-IH and MSW-NH scenarios and lines for the free P ē ē, corresponding to confidence levels of 68% (dashed) and 99.7% (solid).Note that minimum χ 2 min = −2 log L max is the absolute one among all the curves/conversion scenarios.

Fig. 8
Fig. 8 Same as figure 7, but including the accretion component with parameters new parameters T c , τ c , R c , T a , τ a , M a .

Fig. 10
Fig.10Intrinsic (solid curves) and total efficiencies (dashed and dotted curves) for each detector.

Table 1
Range and best-fit (BF) for all parameters in the timedependent model Only Cooling.We show the best-fit for three flavor conversion scenarios: MSW with NH, MSW with IH, and a modelindependent free P ee .

Table 4
Characteristics of each detector