IceCube bounds on sterile neutrinos above 10 eV

We study the capabilities of IceCube to search for sterile neutrinos with masses above 10 eV by analyzing its νμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\mu $$\end{document} disappearance atmospheric neutrino sample. We find that IceCube is not only sensitive to the mixing of sterile neutrinos to muon neutrinos, but also to the more elusive mixing with tau neutrinos through matter effects. The currently released 1-year data shows a mild preference, between 0.75 and 3σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} depending on the binning and flux adopted, for non-zero sterile mixing. This hint overlaps with the favored region for the sterile neutrino interpretation of the ANITA upward shower although the null results from CHORUS and NOMAD on νμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\mu $$\end{document} to ντ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\tau $$\end{document} oscillations in vacuum exclude this interpretation, while through a different channel and using a different energy range. At the 99% C.L. an upper bound is obtained that improves over the present Super-Kamiokande and DeepCore constraints in some parts of the parameter space. We also investigate the physics reach of the roughly 8 years of data that is already on tape as well as a forecast of 20 years data to probe the present hint or improve upon current constraints.


Introduction
Over the last 20 years neutrino oscillations have been established as the explanation of the experimental evidence for neutrino flavor transitions [1,2] with the mixing angles and mass squared differences measured to high accuracy (see  Table 1 for recent global fit results of the mass and mixing parameters). The simplest extension of the Standard Model accommodating neutrino masses is the addition of sterile (right-handed) neutrinos to its content. The mass of these extra singlets, unlike all other fermions, would not be related to the Higgs mechanism and the electroweak scale due to their singlet nature. Therefore, it is vital to probe their existence experimentally at all possible scales.
For instance, short baseline (SBL) experiments like LSND [3] and MiniBOONE [4,5] as well as reactor experiments combined with a recent reevaluation of their expected fluxes [6,7] and Gallium source experiments [8][9][10] have reported oscillation results that are consistent with a mass squared difference with a possible fourth neutrino mass eigenstate of m 2 41 ≈ 1 eV 2 , although this interpretation is in strong tension with other searches [11,12]. With an even higher mass, around the keV scale, sterile neutrinos are a viable dark matter candidate [13,14] that can be probed via their decay to light neutrinos and X-rays with both stringent constraints [15] and a possible hint at 3.5 keV [16,17], which could come from the decay of such neutrinos. At larger masses, sterile neutrinos would leave their imprint altering the kinematics of beta decays and meson decays and can also be probed for at beam dump and collider experiments [18,19]. Even when beyond the reach of collider searches, sterile neutrino mixing can be tested indirectly via precision electroweak and flavor observables [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].
In this study we will consider larger sterile mass squared differences and investigate the sensitivity to the mixing of the sterile neutrinos to the μ and τ flavors of the presently released 1 year data, as well as forecasts for 8 years and 20 years, of atmospheric muon neutrino disappearance data at IceCube. In particular, we will study mass squared differences large enough for the sterile-neutrino-driven oscillations to be averaged out at IceCube energies ( m 2 41 100 eV 2 ) for atmospheric neutrinos traveling through the Earth (L 12000 km). Such mass squared differences are too big to explain the SBL anomalies but are compatible with the sterile neutrino interpretation [50] of the upward directed cosmic ray shower observed by ANITA [51]. These mixings are however ruled out by cosmological constraints [52] and some non-standard effect suppressing the production of these sterile neutrinos in the early Universe would be necessary to reconcile the results [53]. Our results apply for sterile neutrino masses m 2 41 100 eV 2 . Note that for sterile masses above 10 MeV stronger bounds on the active-heavy mixing with muon and tau neutrinos are present from laboratory experiments where the sterile neutrino could be detected directly [18]. 1 Super-Kamiokande [40] and DeepCore [54] have already published constraints on the sterile mixing parameters in the averaged out regime. It should be noted that the energy threshold of Super-Kamiokande is lower than the IceCube one and the averaged out regime for Super-Kamiokande therefore starts at smaller mass squared differences ( m 2 41 > 10 −1 eV 2 ). The same parameter space has also been probed by experiments like CHORUS [55] and NOMAD [56] although, instead of analyzing the disappearance of atmospheric ν μ and the effect of the matter potential from neutral current interactions in presence of steriles, they searched for the appearance of ν τ in a ν μ beam through vacuum oscillations. We will compare our results to these current experimental bounds.
This paper is organized as follows: In Sect. 2 we give an overview of the muon neutrino survival probability when the oscillations driven by the new mass eigenvalues are averaged out, in Sect. 3 we analyze one year of though-going muon data in IceCube and give forecasts for the 8 years and 20 years sensitivities. Finally, in Sect. 4 we summarize our results and give our concluding remarks.

Sterile neutrino mixing
Upon the addition of several sterile neutrinos, the flavor eigenstates of the weak interactions |ν α (α = e, μ, τ, s 1 , s 2 , . . .) are related to the neutrino mass eigenstates |ν i , with masses m i (i = 1, 2, 3, 4, 5, . . .) via the elements U αi of the lepton mixing matrix according to In general, the mixing matrix for n neutrino flavors can be decomposed as the product of n(n −1)/2 rotations with mixing angles θ i j , with (n −1)(n −2)/2 physical phases δ i j . The usual parametrization is through a series of unitary rotations V i j in the i-j-plane given by with and where U 0 has the usual PMNS matrix, U ν , as the upper left 3 × 3 block. Note that we have not included rotations in the purely sterile sector, e.g., V 45 , as such rotations are unphysical. Written in this fashion, the full mixing matrix takes the block form where 1 − α, , X , and Y represent the different blocks of U, and U 0 is assumed block diagonal with no active-sterile mixing. Here, if the rotations are performed in the order given by Eq. (2), α is a lower triangular matrix of the form [58][59][60][61][62] whose components to leading order in the active-heavy mixing elements are given by The ν μ → ν μ oscillation probability, P μμ , will here be derived and discussed in the case where the active-heavy mixing angles are small, the corresponding mass squared differences are large enough for the oscillations to average out, and where the electron neutrinos do not participate in the oscillations 2 (i.e., m 2 21 L/2E 1 and θ 1i = 0). We do so by considering a basis that is rotated by U relative to the flavor basis. In this basis, the neutrino oscillation Hamiltonian in matter takes the form where H 0 is the standard Hamiltonian for μ-τ oscillations in vacuum, H 1 is a diagonal matrix containing large entries, and V NC = ∓G F N n / √ 2 (with the upper sign for neutrinos and the lower sign for anti-neutrinos). The upper left 2 × 2 block describing the ν μ -ν τ oscillations (not including the electron neutrino states) can be treated separately, leading to the 2 × 2 effective Hamiltoniañ where the represents equality up to a matrix proportional to unity and to leading order in α. This can be rewritten as where 2 The oscillation of ν μ to ν e at the energies and baselines that characterize the IceCube data are strongly suppressed. Indeed, θ 13 is small and the solar mass squared difference m 2 21 is too small for the oscillations with θ 12 to develop. As a simplifying assumption, we also set θ 14 = 0, noting that there are strong bounds on |U e4 | 2 10 −3 for m 4 and λ is a phase factor of modulus one. Rotating back to the flavor basis, the muon neutrino survival probability is given by where the last term is a constant leaking term [63]. Note that, except for the leaking term, all the sterile neutrino effects are encoded in the matrix α, in particular in the elements α μμ , α τ τ , and α τ μ , regardless of how many sterile neutrinos are considered as long as they are all in the averaged out regime [62]. However, in our analysis of IceCube data we will allow a free normalization of the events, given the large uncertainties in the atmospheric neutrino fluxes, thus we expect limited or no sensitivity to the normalization factor (1 − α μμ ) 4 and to the leaking term, which does not depend on energy nor baseline. The simulations presented in the next section have been performed with the full numerical oscillation probability and confirm this expectation. At leading order in α, and neglecting m 2 31 whose effect is negligible at the energies of the IceCube data sample, the following probability is obtained where the overall normalization has also been dropped since we allow a free normalization in the analysis. In order to ease the comparison with existing constraints from Super-Kamiokande [40] and DeepCore [54] and to make use of the The baseline has been set to the diameter of the Earth nuSQuIDS software [64,65] for numerical calculations without approximations, we will now particularize these expressions for the addition of a single sterile neutrino. With our given parametrization, we find that U μ4 = s 24 e −iδ 24 and U τ 4 = c 24 s 34 , so that and thus Therefore, the bounds will essentially follow a hyperbola in the |U μ4 | 2 -|U τ 4 | 2 -plane. Note that, in contrast to IceCube, for the Super-Kamiokande and DeepCore energies the atmospheric oscillation driven by m 2 31 is relevant. Thus, the approximate Eq. (15) is not valid and the sensitivity mainly stems from the interference between the standard and sterile oscillations in Eq. (10). Therefore, the phase of α τ μ , i.e., δ 24 in the one extra sterile neutrino scenario, has an impact on the oscillation probability. Specifically, it can change the sign of the interference term between the atmospheric and the sterile terms in the expression for the energy and the mixing angle in matter. As an example of the impact of the phase, in Fig. 1 the muon neutrino survival probability as a function of the energy for |U μ4 | 2 = 10 −1.6 , |U τ 4 | 2 = 0.15, L = 1.2×10 4 km, and two different values of the phase, δ 24 = 0 (solid line) and δ 24 = π (dashed line) is shown. As comparison, the muon neutrino disappearance oscillation probability for zero sterile mixing is also shown. These values of the sterile matrix elements are at the border of the 90% C.L. region of Super-Kamiokande. The sign of the interference term can also be changed by changing the mass ordering (i.e., the sign of m 2 31 ) or by switching between neutrinos and antineutrinos (i.e., changing the sign of V NC ). However, neither IceCube nor Super-Kamiokande or DeepCore can distinguish between neutrinos and antineutrinos so this dependence is diluted in their data.
Conversely, experiments such as CHORUS and NOMAD explored the same parameter space but instead exploiting the ν μ to ν τ appearance channel with negligible matter effects leading to

Simulation and results
One year of high-energy through-going muons released by the IceCube collaboration [39] for the last IceCube detector stage with 86 strings will be analyzed. The data sample consists of up-going track events so as to avoid the background from cosmic ray muons giving, after all cuts, a sample purity better than 99.9%. Hence, the distances the signal neutrinos travel are of the order of 10 4 km. The selected events have reconstructed energies between 400 GeV and 20 TeV and cosine of the reconstructed zenith angle between −1 and 0.2. The sensitivity that a full 8-year IceCube sample would have as well as the prospects for an exposure equivalent to 20 years of IceCube data will also be forecasted. For our simulations, the neutrino flux computed with the analytic air shower code [66] using the cosmic ray flux from HondaGaisser model with Gaissser-Hillas H3a correction [67] together with the hadronic model QGSJET II-04 [68] have been adopted. We have also verified that our results do not change significantly under the assumption of different fluxes, such as using the cosmic ray flux from the poly-gonato model [69,70] or the Zatsepin-Sokolskaya [71] model updated with measurements by PAMELA [72] together with the hadronic model SIBYLL2.3, RC1, point-like [73] or QGSJET II-04. The propagation of the neutrinos was simulated using the nuSQuIDS software [64,65], where the PREM profile [74] is implemented for the Earth matter density. Since we are interested in the averaged out regime our simulations were performed with a sterile mass squared difference of m 2 41 = 10 3 eV 2 , but we have verified that changing this parameter does not alter the results as long as m 2 41 100 eV 2 as expected.
Since neutrino and antineutrino interactions cannot be distinguished on an event basis, the signal will contain both ν μ andν μ events. After propagating the flux for every value of the sterile neutrino parameter, the Monte Carlo provided with the data release [39] has been used to compute the expected number of events N th,i in every bin of reconstructed zenith angle.
In order to obtain the expected significance of the bounds on the sterile mixing parameters, we adopt a Poisson loglikelihood given by where the N th,i and N d,i are the predicted and observed number of events given a set of parameters in bin i, respectively, and the sum is taken over all the reconstructed zenith angle bins i. The log-likelihood has been maximized for a number of nuisance parameters to include the effect of possible systematic errors. In particular, the uncertainty in the pion-kaon ratio of the initial flux (π/k), the efficiency of the digital optical modules (DOMs), and the overall flux normalization have been considered. Since the observable is energy independent for large values of the sterile neutrino mass (see Eq. (15)), only one energy bin has been considered and the uncertainty in the energy spectrum slope has been neglected, while 40 bins for the reconstructed zenith angle have been adopted. For the pion-kaon ratio a Gaussian prior with σ π/k = 0.05 has been adopted and no prior for the DOM efficiency or the overall flux normalization has been assumed. The standard oscillation parameters used in the simulations were set to their respective best-fit values from Table 1. To find the confidence regions from the log-likelihood differences we assume that the prerequisites for Wilks' theorem [75] holds so that likelihood ratios can be directly converted to a confidence level.
In the left panel of Fig. 2, the 90% C.L. constraints (for 2 degrees of freedom) obtained for the public 1-year data (pink contours) in the |U μ4 | 2 -|U τ 4 | 2 -plane is presented. The existing bounds from Super-Kamiokande [40] and DeepCore [54] at the same C.L. are also shown for comparison by the hatched gray area. At 90% C.L. present data prefer some degree of sterile mixing and we find that zero sterile mixing is disfavored at 2.3σ (1 degree of freedom. 3 ) The preference for non-zero sterile mixing is independent on the atmospheric sterile neutrino flux adopted in the analysis but its significance varies between 1.6 and 3.0 σ with the different models tested. Given this preference for non-zero sterile mixing, the current constraints from IceCube do not improve upon the combined bounds from Super-Kamiokande and DeepCore at 90% C.L. In the right panel, the same information is shown at 99% C.L. In this case, the present 1-year data gives an upper bound that already slightly improves upon the present Super-Kamiokande and DeepCore constraints, ruling out the white region in the plot. We have also checked how these results depend on the binning in energy and zenith angle. We have seen that when the data is also binned in energy, the case of no active-sterile mixing becomes slightly less disfavored. Using different combinations of fluxes and binning, the no-mixing scenario is disfavored at between 0.74 and 3.1 σ , depending on the combination. In particular, for the case of 10 energy bins and 21 bins in the zenith angle, as presented in Ref. [39], the significance varies between 0.75 and 3.0 σ . 4 The physics reach of an 8-year run of IceCube data if the present preference for sterile mixing is maintained is also shown in cyan. In particular, the present best-fit value of |U μ4 | 2 = 10 −2 , |U τ 4 | 2 = 0.29 lies in the already disfavored region by DeepCore and Super-Kamiokande. Due to the hyperbola-shaped degeneracy of the oscillation probability in the |U μ4 | 2 -|U τ 4 | 2 -plane, there are values of the sterile oscillation parameters that provide an almost equally good fit without being in tension with the other ν μ disappearance present data. Remarkably, theses values of U τ 4 are also compatible with the sterile neutrino interpretation [50] of the upward directed cosmic ray shower observed by ANITA [51]. Indeed, the sterile neutrino interpretation of the ANITA results requires that the sterile neutrino mass is between ∼ 10 2 and ∼ 10 6 eV, which would also fall in the averaged out regime for IceCube studied here. However, all the parameter space preferred by IceCube at the 90% C.L. is disfavored by NOMAD [56] with the same significance. Indeed, the null results in their ν τ search translates through Eq. (16) into |U μ4 | 2 |U τ 4 | 2 < 8.3 · 10 −5 at the 90% C.L. for m 2 41 100 eV 2 . Nevertheless, the channel and underlying physics explored to obtain the bounds are very different in the two sets of experiments. While Super-Kamiokande, DeepCore and IceCube analyze ν μ disappearance and the steriles are probed via their matter effects as shown in Eq. (15), NOMAD and CHORUS searched for ν τ appearance essentially in vacuum through Eq. (16). Thus, a possible means to reconcile the two results could be through neutrino non-standard interactions (NSI) in matter, whichexcept for the leaking term-are phenomenologically equivalent to averaged-out sterile neutrino oscillations in Ice-Cube [62,63]. The necessary size of the NSI would be at the level of |ε μτ | 2 10 −3 , which is close to the current bound [77]. We therefore simulate 8 years of IceCube data assuming |U μ4 | 2 = 10 −2 , |U τ 4 | 2 = 0.1, and δ 24 = 0 as the true oscillation parameters. As can be seen in Fig. 2, the expected confidence region shrinks significantly with the  [40] and DeepCore [54] (NOMAD [56]) data at the same C.L additional statistics, while keeping its shape. In particular, if the values of the sterile neutrino mixing marked by the star were realized in nature, 8-years of IceCube data would disfavor no sterile mixing around the 5σ level. The capability of larger IceCube samples to improve the present constraints on sterile mixing in absence of sterile neutrinos have been also studied. In Fig. 3, the contours for 90% (left panel) and 99% C.L. (right panel) expected exclusion limits in the |U μ4 | 2 -|U τ 4 | 2 -plane together with the existing bounds from Super-Kamiokande and DeepCore are presented. The bound on |U μ4 | 2 from 8 years of IceCube would improve over present constraints between a factor 1.3 for vanishing values of |U τ 4 | 2 to around an order of magnitude for |U τ 4 | 2 close to 0.1. Similarly, for |U μ4 | 2 ∼ 10 −2 , the constraint on |U τ 4 | 2 would improve around a factor 5. In particular, the present best fit for non-zero sterile mixing would be excluded at high significance (more than 5σ ) and most of the currently preferred parameter space at 90% C.L. (pink area in the left panel of Fig. 2) disfavored. Comparatively, increasing the statistics up to 20-year of IceCube data yields a more modest improvement in sensitivity. Remarkably, not even the 20-year scenario would improve over the present NOMAD limit of |U μ4 | 2 |U τ 4 | 2 < 8.3 · 10 −5 at the 90% C.L. Nevertheless, we consider the two constraints complementary given the different physics probed by each of them.
The effect of the CP-violating phase δ 24 is also shown. In particular, the solid lines correspond to δ 24 = 0 and the dashed lines to δ 24 = π . As can be seen, IceCube is not very sensitive to the sterile phase as oscillations due to the atmospheric mass squared difference at energies above 100 GeV do not have time to develop. Indeed, from Fig. 1 the ν μ survival probability is essentially 1 in absence of sterile mixing for E > 100 GeV.

Summary and conclusions
In this work we have presented the current constraints from the public 1-year IceCube data as well as the expectations of a full 8-year dataset and forecasts for 20 years worth of statistics to the mixing of sterile neutrinos with the μ and τ flavors. In particular, we concentrated for the first time on larger masses for the extra neutrinos ( m 2 > 100 eV 2 ) than usually explored so that their oscillations are averaged out at IceCube. We find that the public 1-year IceCube data presents some preference for non-zero sterile mixing in the averaged out regime that would manifest via neutral-currentinduced matter effects in the ν μ disappearance channel. In particular, values of the squared sterile mixing with the μ flavor of order 10 −2 and with the τ between 10 −1 and 10 −2 are favored at around 2σ with respect to no sterile mixing when no binning in energy and 40 bins in zenith angle are adopted. Interestingly, the large masses assumed in our analysis and the size of the preferred mixing with tau neutrinos correspond to the region of the parameter space that could also explain the upward directed cosmic ray shower observed (δ 24 = π ). The solid (dashed) gray hatched regions are disfavored by Super-Kamiokande [40] and DeepCore [54] (NOMAD [56]) data at the same C.L by ANITA [51] with sterile neutrinos [50]. These mixings are however in strong tension with cosmological constraints [52] and some non-standard effect suppressing the production of these sterile neutrinos in the early Universe would be necessary to reconcile these results [53]. Moreover, these mixings are also in tension with present data from CHORUS [55] and NOMAD [56] which, however, explore a different channel without matter effects. Thus, in presence of non-standard matter effects the two results could be potentially reconciled.
We have also studied the sensitivity that 8 years of Ice-Cube data, close to the data that should be presently available, would have and find that it would be sufficient to either confirm the present preference or exclude it with high significance (more than 5σ ) and set stringent constraints improving around an order of magnitude over Super-Kamiokande and DeepCore present bounds in some parts of the parameter space. Since sterile neutrinos at some mass scale are a general expectation of many extensions of the SM accounting for neutrino masses, it will be very interesting to explore this part of the parameter space with averaged out sterile neutrino oscillations using the full data sample collected by IceCube.