On the robustness of IceCube’s bound on sterile neutrinos in the presence of non-standard interactions

The mixing parameters of sterile neutrino(s) preferred by the MiniBooNE and LNSD experiments are in strong tension with the exclusion limit from the IceCube experiment. Recently it has been claimed that by considering the non-standard neutrino interactions (NSI) in addition to the sterile neutrino, the IceCube’s limit can be relaxed and the tension can be reconciled; a baroque scenario as it has been called. We will show that this claim is just an artifact originating from the energy cuts of the chosen datasets. Contrary to the claim, by turning on the NSI and fixing the NSI parameters to the proposed values, not only the IceCube’s limit on sterile neutrino cannot be alleviated, but in fact the tension will be aggravated (or at least keeps its strength). To this aim, an analysis of the IceCube’s atmospheric neutrino data in the full energy range is crucial.


Introduction
Currently, almost all the neutrino data can be explained consistently in the 3ν formalism, consisting of three active neutrino flavors and the corresponding mixing parameters (for a global fit to all the available data see [1][2][3]). However, there are still some anomalies, coming from the electron neutrino appearance experiments, namely, the LSND [4] and MiniBooNE [5] experiments, which indicate the existence of extra neutrino states, the so-called sterile neutrinos, with the mass ∼ O(1) eV. In particular, the recent update from the MiniBooNE experiment [6], which combines the ν e andν e appearance data, reports an excess of 4.8σ in the low energy range that can be increased to 6.1σ if combined with the LSND data. This excess can be interpreted in the 3 + 1 scenario (3 active + 1 sterile neutrino state) with Δm 2 41 ∼ O(1) eV 2 and sin 2 2θ eμ 10 −2 . The a e-mail: arman@puc-rio.br b e-mail: nunokawa@puc-rio.br allowed region in the (sin 2 2θ eμ , Δm 2 41 ) plane can be found in [6].
The main obstacle in a happy interpretation of the LSND/MiniBooNE excess in terms of the sterile neutrinos is the strong tension with the disappearance data including MINOS/MINOS+ [7] and IceCube [8] experiments (for a global status of the sterile neutrino mixing from various experiments see [9,10]). Among these the IceCube's limit has a different nature: while the MINOS experiment is sensitive to the sterile neutrino mixing through an averaged effect over the long baseline, the IceCube sensitivity originates from a resonance effect, an amplification of theν μ →ν s oscillation probability (or ν μ → ν s for Δm 2 41 < 0) for atmospheric neutrinos crossing the Earth in the ∼ TeV energy range. The possibility of exploring the active-sterile neutrino mixing by looking at the energy and zenith angle distributions of the high energy atmospheric neutrinos has been proposed in [11,12] (see also [13]). By the realization of the IceCube detector, as the first km 3 -volume neutrino telescope which is able to detect ∼ TeV atmospheric neutrinos, this possibility has been studied in detail: the limit on eV-scale sterile neutrinos has been derived using the data collected during the construction phase of the IceCube [14] and it has been shown that few-years worth of the IceCube data can exclude completely the preferred region by LSND and MiniBooNE [15]. The effect of the various mixing parameters in the 3 + 1 scenario has been studied in [15] and the search strategy generalized to the cascade event topology in [16] (see also [17][18][19][20] for further analyses). Finally the IceCube collaboration published the result of the sterile neutrino analysis [8], considering one-year muon-track data with the reconstructed muon energy proxy in the range 400 GeV − 20 TeV, demonstrating the exclusion of the preferred parameter space region by the appearance experiments (to be precise, constraining the angle θ 24 in the 3 + 1 scenario).
Another new physics scenario that can be probed by the high energy atmospheric neutrinos observed by the IceCube is the non-standard neutrino interactions. This possibility has been proposed and used to derive the most stringent bound on the NSI parameters (the ε μτ and ε μμ − ε τ τ ) in [21], with consistent results in [22,23] (see [24] for bounds on the NSI parameters from the global analysis of oscillation data). Recently it has been proposed in [25] that the addition of non-standard neutrino interaction to the 3 + 1 picture can relax the limit of IceCube on sterile neutrinos and reconcile the appearance and disappearance discrepant results. The same claim has been repeated in a more recent work [26] which analyses the data of MINOS+ and IceCube, including the DeepCore data, and concludes that a combination of the charged-current and neutral-current NSI can relax the limits of both experiments. In this framework, the charged-current NSI is required for the relaxation of the MINOS+ limit (a nonzero detector NSI, ε D μμ , has been assumed) while the neutral-current NSI parameters (nonzero ε μμ , ε τ τ and ε ss ) will, pretendedly, loosen the IceCube's limits. The IceCube data analyzed in [26] consists of the publicly available sterile search data 1 [8] with the muon energy 2 ∈ [501 GeV, 10 TeV] and the DeepCore oscillation data 3 [28] with the muon energy range [6,56] GeV.
In this paper we will argue that in fact the apparent relaxation of the IceCube's limit on the sterile neutrino in the (3 + 1)+NSI scenario originates from the negligence of the [56, 500] GeV energy range and has no physical rationale behind. It should be emphasized that although the recent Ice-Cube publicization of data on atmospheric neutrinos does not include the [56, 500] GeV energy range, older data are available, such as the IC-79 [29] and IC-40 [30], which cover this energy range and do not show any significant deviation from the 3ν framework expectation. Considering this energy range, contrary to the claim in [25,26], the limit on sterile neutrino from IceCube will be stronger in the (3 + 1)+NSI scenario or will keep the strength, depending on the quality of the data in the [56, 500] GeV energy range. To this aim, actually, a simple oscillation probability calculation is enough to manifest the argument.
The paper is organized as follows: In Sect. 2 we summarize the main features of the atmospheric muon (anti-)neutrino oscillation in the (3 + 1)+NSI scenario and describe our assumptions. In Sect. 3 we discuss in detail, based mainly on oscillation probabilities, why the addition of NSI to the 3 + 1 model cannot help to the reconciliation of the LSND/MiniBooNE and IceCube tension. Finally, in Sect. 4, we provide our conclusions. 1 https://icecube.wisc.edu/science/data/IC86-sterile-neutrino. 2 The public data is available in the muon energy range [400 GeV, 20 TeV]. 3 https://icecube.wisc.edu/science/data/2018nuosc.

The (3 + 1)+NSI: a concert of new physics scenarios
The phenomenology of high energy atmospheric neutrino oscillation in the 3 + 1 and NSI frameworks have been studied separately in [15,21]. The characteristics of oscillation in the (3 + 1)+NSI is just a simultaneous consideration of both frameworks. In the following we will summarize the evolution equation of neutrinos in the (3 + 1)+NSI scenario, in order to fix the notation and remind the main features.
The evolution equation of the neutrinos, in the flavor basis, in the presence of matter effect, in the (3 + 1)+NSI scenario can be written as: where H 0 corresponds to the Hamiltonian in vacuum, and is given by Here E ν is the neutrino energy, 3,4) are the mass-squared differences, and U 3+1 is the mixing matrix for the 3 + 1 model, parameterized as where R i j (θ i j ) is the rotation matrix with the angle θ i j in the i-j plane; while the rotation-like matrices R i j (θ i j , δ i j ) can be obtained from R i j (θ i j ) by the following replacements: The V m term in Eq. (1) represents the matter potential induced by the standard (coherent) interaction as well as by the NSI, which can be parameterized generally as where G F is the Fermi constant, n e is the electron number density of the propagation medium (in our case is the Earth), the ε αβ are the dimensionless parameters characterizing the strength of the NSIs, 4 and κ ≡ n n /(2n e ) ∼ 1/2 (n n being the neutron number density) for the Earth's matter. For antineutrinos, the overall sign of the matter potential in Eq. (1) as well as that of the CP violating phases in Eq. (3) must be flipped.
The oscillation pattern of the high energy atmospheric neutrinos in 3+1 scenario is well-known: assuming Δm 2 41 > 0, where otherwise a strong tension with the cosmological data will appear, a resonance enhancement of thē ν μ →ν s oscillation occurs as neutrinos pass the Earth's core and/or mantle as follows. For core-crossing trajectories, that is cos θ z −0.8 where θ z is the zenith angle, the enhancement is a parametric resonance [31,32] at the energy (2.5 TeV) cos 2θ 24 (Δm 2 41 /eV 2 ); while for the mantlecrossing trajectories (cos θ z −0.8) the enhancement is an MSW resonance at the energy (4 TeV) cos 2θ 24 (Δm 2 41 / eV 2 ). See [33] for an elucidation of these effects.
The resonance enhancement ofν μ →ν s conversion is independent of the θ 14 , and as shown in [15] the enhancement becomes more effective with the increase of θ 34 . In what follows, we set sin 2 θ 14 = 0.02 and θ 34 = 0 where the latter choice leads to the most conservative limit on sterile neutrino mixing. For the standard 3ν-oscillation parameters we use the best-fit values of the global fit [1] for the normal mass ordering, Δm 2 31 > 0, while our discussion will be equally valid also for the inverted mass ordering. It should be noticed that in the high energy range (> 100 GeV) no standard oscillation effect is important and by decreasing the energy to ∼ 20 GeV the atmospheric oscillation parameters start to play a role. All the CP-violating phases also will be set to zero, as recommended in [15] for a conservative limit.
In the presence of the NSI the resonance energies will be modified. To simplify the discussion, and also to consider the same setup as in [25,26], we will assume nonzero (ε μμ , ε τ τ , ε ss ) and setting all the other NSI parameters to zero. By neglecting the ν e -flavor oscillation, see [15], through a simple calculation of the resonance condition in the (3+1)+NSI scenario we can obtain the resonance energy of μ − s conversion. For the mantle-crossing trajectories the MSW resonance energy will be modified to E res ν 4 cos 2θ 24 Δm 2 41 1 eV 2 TeV. (5) For the core-crossing trajectories, where cos θ z −1, the numerical factor of Eq. (5) should be replaced by 2.4 TeV. As can be seen, the nonzero ε μμ and/or ε ss shift the resonance energy. The resonance energy is independent of the ε τ τ . Thus, neglecting for the moment the effect of the NSI parameters in the low energy, that is E ν ∼ (10 − 100) GeV, by turning on the neutral-current NSI the resonance will occur at some shifted energy, but does not disappear. 5 The effect of the NSI Table 1 The two sets of parameters in the (3 + 1)+NSI scenario considered in this work, which come from a scan of the parameter space performed in [26]. All the other parameters of the 3+1 model and/or NSI parameters not indicated in the table are assumed to be zero. parameters on the high energy atmospheric neutrino, including the low energy (10 − 100) GeV range, has been studied in detail in [21] and we will not repeat it here. The atmospheric oscillation probabilities in low energy range depend dramatically on the ε μμ − ε τ τ . Thus, obviously, a large value inserted for ε μμ in Eq. (5) should be compensated with a large value for ε τ τ such that their difference remains small and the fit to the DeepCore data does not deteriorate. This is exactly the case chosen in [25,26] and we will take it also in this paper.

Oscillation probabilities in the (+ 1)+NSI scenario
The oscillation probabilities P(ν μ → ν μ ) and P(ν μ →ν μ ) can be obtained by the numerical solution of the neutrino evolution equation in the (3 + 1)+NSI scenario; i.e., the Eq. (1) and the corresponding one for the anti-neutrinos. In the numerical solution we will use the PREM model [27] for the Earth matter density profile assuming Y e = 0.5, where Y e denotes the number of electrons per nucleon. By looking at the (3 + 1)+NSI oscillation probabilities, in this section we will argue and show that the relaxation of the IceCube's bound on the sterile neutrinos, as claimed in [25,26], does not happen. To be specific, we will consider two different sets of NSI parameters, motivated by and studied in [26], which come from a scan over the NSI parameter values. These two sets, named case (a) and case (b), are summarized in Table 1. For the case (a) the ε ss has been set to zero while a scan over |ε τ τ | < 6 and |ε τ τ − ε μμ | < 0.5 has been done. 6 For the case (b) the scan is over |ε ss | < 6, |ε τ τ | < 0.5 and |ε τ τ −ε μμ | < 0.5. As these two cases have been reported as the best scenarios of an scan over the NSI parameter values, any other choice of the NSI parameters will be less effective in relaxing the IceCube's bound according to [26].
The sterile neutrino mixing parameters in the cases (a) and (b) have been fixed to the reported best fitted values in Footnote 5 continued anomalies in the cascade-type events of IceCube and can be constrained. We will leave this setup for a future study 6 Although this set of values for the NSI parameters are in conflict with the solar neutrino data [24], we will continue with it as an example of the claimed maximum relaxation of IceCube's bound reported in [26] [26]. Although, the P(ν μ → ν μ ) and P(ν μ →ν μ ) oscillation probabilities do not depend on θ 14 , we will just set it to sin 2 θ 14 = 0.02 in our numerical calculation.
Let us start with the case (a). The upper and lower panels of the Fig. 1 show, respectively, the P(ν μ → ν μ ) and P(ν μ → ν μ ) oscillation probabilities for mantle-crossing trajectory cos θ z = −0.8. The black thick dashed curve shows the oscillation probability in the 3ν framework. The red solid curve shows the oscillation probability in the 3+1 model with sterile mixing parameters as in case (a); i.e., the mixing parameters shown in the first row of Table 1 and setting the NSI parameters to zero. As expected, the Δm 2 41 = 0.32 eV 2 leads to an MSW resonance at E ν 1.2 TeV in the anti-neutrino channel. The IceCube's sensitivity to the mixing parameters of case (a) originates mainly from this dip in the anti-neutrino oscillation probability. By turning on the NSI parameters, the oscillation probability shown by blue dashed curve will be modified. The resonance dip inν μ →ν μ oscillation, in accordance with Eq. (5), shifts to ∼ 200 GeV. In Fig. 1 the green and pink shaded regions show, respectively, the energy ranges of the IceCube's sterile neutrino analysis [8] and the DeepCore oscillation analysis [28]. For comparison, the oscillation probability for 3ν+NSI scenario, with NSI parameter values of case (a), is also shown by the brown dotdashed curve. As can be seen, addition of the sterile neutrino will makes the low energy part of the 3ν+NSI oscillation probability compatible with the 3ν, that otherwise is completely ruled out. Obviously, if one considers just the Deep-Core and IceCube sterile neutrino analyses energy ranges, the limit on sterile neutrino in the (3 + 1)+NSI scenario can be relaxed; in the pink and green shaded regions there is a small difference between the 3ν and (3 + 1)+NSI oscillation probabilities. However, the huge discrepancy between the 3ν and (3 + 1)+NSI scenarios lies in the gap between the two energy ranges.
The Fig. 2 shows the oscillation probabilities for the case (a) and for the core-crossing trajectory with cos θ z = −1, with the same color code and line type as in Fig. 1. For cos θ z = −1 the muon anti-neutrinos experience the parametric resonance in 3 + 1 scenario at ∼ 2 TeV as expected; while in the (3 + 1)+NSI scenario the resonance shifts to lower energies, again in the gap between the energy ranges of DeepCore and IceCube sterile neutrino analyses. The are more deviations in the low energy part for the core-crossing trajectories due to the higher matter density in the propagation path of the neutrino which amplifies the effect of the NSI. In fact this deviation in the low energy part of the Fig. 2 is the origin of the limit obtained in [26] for case (a); otherwise, by just considering the green shaded region, IceCube sterile neutrino analysis is not sensitive to case (a). Figures 3 and 4 show the oscillation probabilities for the case (b), respectively, for the mantle and core crossing trajectories. The discussions presented for the case (a) apply also for the case (b). Again, the principal effect of the NSI is to shift the resonances into the gap in the energy ranges of the DeepCore and IceCube sterile neutrino analyses.
A comment on the energy ranges of the DeepCore and IceCube neutrino sterile analyses is in order: although the respective [6,56] GeV and [501 GeV, 10 TeV] energy ranges are referring to the muon energy produced in the ν μ andν μ charged-current interactions, we are taking it as a proxy of the neutrino energy and show them in the Fig. 1 (and the  subsequent figures) as cuts on the E ν . Needless to say, this is just roughly correct and in a detailed analysis of the data the difference between the muon energy and E ν should be taken into account. However, it is not our goal in this paper to perform an analysis of the data, basically since there are no publicly available data from IceCube collaboration in the energy range of [56, 500] GeV. Instead, we will do a sensitivity analysis, similar to the one done in [15], just to quantify our argument. The main points of our argument are clear enough that hopefully motivate the IceCube collaboration to release the data in the whole energy range and perform an analysis by taking into account all the details.
Even without an elaborate calculation it is clear that a dip in the oscillation probability at ∼ 300 GeV can be excluded more strongly than or at least at the same level of a dip at ∼ 2 TeV; which is what happening in (3 + 1)+NSI scenario (see, for example, the lower panel of Fig. 1). The atmospheric muon neutrino flux drops roughly as ∝ E −3.7 ν with the increase of energy. By taking into account the energy dependence of the cross section and the increase of effective volume by the increase in energy, the statistics is higher roughly by a factor of few at E ν ∼ 500 GeV with respect to E ν ∼ 2 TeV. For example, in the IC-79 dataset the peak of the statistics is at ∼ 500 GeV, and the dataset contains as much events at ∼ 300 GeV as in ∼ 2 TeV (see Fig. 1 of [29]). Other elements, such as the better energy resolution in the low energy due to the shorter muon range, also help toward a better constraining power. Thus, for sure, one can conclude that the (3 + 1)+NSI for both cases (a) and (b) is more strongly excluded than the 3 + 1 model with the same sterile mixing parameters (or pessimistically it is excluded at the same level). This can be seen from the figures in this paper: by adding the NSI to the 3+1 model not only the resonance dip in the muon neutrino survival oscillation probabilities slides to ∼ 300 GeV but also the oscillation probabilities will be modified at the DeepCore energy range. The sum of these two can lead to a strong bound from IceCube data. A detailed quantification of this statement is not possible for us since there is no recent public data from IceCube in the gap energy range. Thus, let us quantify a bit this statement by performing a sensitivity analysis using the information from the older IceCube datasets.
For our sensitivity analysis we use the last public effective area of the IceCube in the whole range of energy, which  Table 1. For completeness, the case where only the NSI effect is added to the standard 3ν oscillation is also shown by the brown dot-dashed curve. The energy ranges used by the IceCube's sterile neutrino analysis [8] and the DeepCore oscillation analysis [28] are indicated by the green and pink shaded regions, respectively goes back to the construction period of the detector: the IC-40 [30] and IC-79 [29] configurations. Of course, at least 10 times more data are available now and so we just increase the data-taking period correspondingly. For fixed values of the (sin 2 2θ 24 , Δm 2 41 ) we calculate the χ 2 value (which is the same as Δχ 2 with respect to the 3ν framework) by marginalizing the following χ 2 function over the parameters α and β: χ 2 (Δm 2 41 , sin 2 θ 24 ; α, β) = Fig. 3 The same as Fig. 1 but for the case (b) shown in Table 1 i, j where α and β are the pull parameters taking into account the correlated uncertainties of the atmospheric neutrino flux normalization and its zenith dependence (tilt), respectively. Using the Honda flux of atmospheric neutrinos [34], these uncertainties are roughly σ α = 0.24 and σ β = 0.04. The  Table 1 N 0 i, j is the number of events in the 3ν scenario in the i th bin of energy and j th bin of cos θ z , which can be calculated by convoluting the effective area with the neutrino flux and oscillation probability. The N i, j is the number of events in the (3 + 1)+NSI scenario with the sterile mixing parameters (sin 2 θ 24 , Δm 2 41 ). The σ i, j,stat = N i, j is the statistical error and we have added an uncorrelated systematic error σ i, j,sys = f N i, j where f quantifies it.
We have verified that by using the [400 GeV, 20 TeV] energy range, four energy bins (uniformly distributed in log), 10 linearly distributed bins of cos θ z in [−1, 0], and an uncorrelated systematic error f = 10% we can reproduce the exclusion plot of the IceCube analysis in [8] almost exactly. Focusing on the case (b), the χ 2 value for the (Δm 2 41 , sin 2 θ 24 ) = (0.63 eV 2 , 0.032) in the 3 + 1 model using the [500 GeV, 10 TeV] energy range is ∼ 12. By turning on the NSI parameters and going to (3+1)+NSI scenario, using the same energy range, the χ 2 value drops to ∼ 1. This is in agreement with the claim in [25,26] that the IceCube's limit on sterile neutrino relaxes in the presence of NSI. However, by extending the energy range to [10 GeV, 10 TeV] the χ 2 value increases to ∼ 26, as we expected. The same pattern, that is a significant increase of the χ 2 value, occurs also for the case (a) if we extend the energy range as done for the case (b).
By fixing the values of NSI parameters to the values of cases (a) and (b) in Table 1, and calculating the exclusion limit for all the (Δm 2 41 , sin 2 θ 24 ) values, there is no point in the (Δm 2 41 , sin 2 θ 24 ) plane with a Δχ 2 value less than ∼ 26 for both the cases (a) and (b), and therefore no allowed region exists in this plane (up to ∼ 5σ C.L.). Of course, this result is expected since the NSI parameters of cases (a) and (b) are large. Even by assuming a very small θ 24 , the (3 + 1)+NSI scenario is excluded just because of the large NSI parameters in cases (a) and (b).

Conclusions
The IceCube's bound on the active-sterile neutrino mixing(s) strongly excludes the parameter space preferred by the appearance experiments LSND and MiniBooNE, such that a global fit of the data in the 3 + 1 scenario shows strong tensions. Recently it has been claimed that by adding nonstandard neutrino interaction to the 3+1 scenario it is possible to relax the IceCube's bound and weaken the tension [25,26].
We have revisited the (3 + 1)+NSI scenario and studied the impact of the NSI on the IceCube's bound on the sterile neutrino. We have shown that in the presence of NSI the resonance enhancement of theν μ →ν s occurs at a shifted energy, but does not disappear. The reason behind the relaxation of the IceCube's bound claimed in [25,26] is that for the chosen NSI parameters the resonance enhancement occurs in the gap [56, 500] GeV between the energy ranges of the two considered datasets (the IceCube sterile neutrino analysis [8] and the DeepCore oscillation analysis [28]). However, although the recent public data of the IceCube do not cover this gap, the older data such as IC-79 [29] and IC-40 [30] include this energy range and do not show any significant deviation from the 3ν oscillation framework. By performing a sensitivity analysis, covering the energy range [10 GeV, 10 TeV], we have shown that in fact the IceCube's bound on the sterile neutrino mixing becomes stronger in the presence of the NSI.
We therefore conclude that the NSI with the parameter values reported in [26] does not relax the IceCube's bound on the sterile neutrino. In the presence of NSI, the tension between the IceCube's limit and the LSND/MiniBooNE preferred region persists or even can get stronger.
A detailed analysis by the IceCube collaboration is required to provide the correct limit on the (3 + 1)+NSI scenario. We hope that this work motivates the IceCube collaboration to publicize the atmospheric neutrino data in the full range of energy to avoid such misinterpretations.