Sterile neutrinos in the light of IceCube

We determine constraints on parameters of a single eV-scale light neutrino using IceCube-59 data. Particular emphasis is put on the question whether such an analysis can rule out sterile neutrino hints. While important complementary information is provided, the different dependence on the various sterile neutrino mixing angles makes it currently not possible to fully exclude short baseline appearance results or sterile neutrinos in general.


Introduction
Though the standard 3-neutrino mixing paradigm has been well established [1], there are still several short-baseline anomalies, most notably in LSND [2], MiniBooNE [3] and reactor neutrino flux measurements [4,5]. These anomalies could be explained by introducing light sterile neutrinos with the mass-squared difference ∆m 2 at 0.1 ∼ 1 eV 2 scale and mixing matrix elements around 0.1. If such light sterile neutrinos would indeed exist, the theoretical implications would be profound. Therefore several experiments are running or in construction in order to confirm or refuse the existence of light sterile neutrinos. See the recent reviews [6,7] for an overview of the hints, consequences and tests of light sterile neutrinos.
In this work we focus on effects of eV-scale sterile neutrinos in atmospheric neutrino oscillations at high energies as measured in the IceCube experiment. We use here an IceCube-59 data set from ref. [8], where a search for diffuse astrophysical neutrinos was performed. The fact that sterile neutrinos would have an impact in IceCube is easily understood by noting that a mass-squared difference of order eV 2 corresponds to maximal oscillations at energies E ν ∼ 10 3 GeV and a baseline around Earth radius R ⊕ 6.4 × 10 3 km. Indeed, atmospheric neutrinos observed in IceCube have energies ranging from 10 2 GeV to 10 6 GeV, peaked at about 10 3 GeV. Several papers have in the past analyzed the effect of light sterile neutrinos at high energies as a potential test of the sterile neutrino hypothesis [9][10][11][12][13][14][15][16].
We will perform here a χ 2 -analysis on the IceCube-59 data within a 3+ 1 scheme to see how significant the constraint on sterile neutrinos is. 1 We are particularly interested in the interplay of this constraint with the results of short baseline appearance and other experiments. We work in the minimal framework of only one sterile neutrino with a unitary 4 × 4 matrix and no additional interactions. Even in this simple situation the dependence on the various lepton mixing angles is different in each experiment, and in addition, muon neutrino 1 Here the 3 + 1 scheme refers to the case that there are 3 active neutrinos and 1 heavier sterile neutrino.
In principle other cases such as 1 sterile neutrino + 3 heavier active neutrinos, 3 active neutrinos + 2 heavier sterile neutrinos, are also possible. However, since in these cases the total neutrino mass Σimi is much larger, they are disfavored by cosmological constraints.

JHEP01(2016)124
disappearance in IceCube depends on the weakly constrained angle θ 34 , that should not be set to zero in analyses. For scenarios with more than one sterile state the larger number of mixing angles will complicate further the direct comparison of oscillation probabilities in IceCube and other experiments. Moreover, there is dependence on the sterile neutrino global fit results to which one compares the IceCube sensitivity. All in all, a full exclusion of short baseline appearance and/or other results is currently not possible, though of course important complementary constraints on sterile neutrinos are provided by IceCube data.
The paper is build up as follows: in section 2 we discuss the procedure to obtain the oscillation probability of high energy neutrinos including matter effects, before discussing event numbers in IceCube in section 3. The numerical analysis of IceCube data is performed in section 4, where we also discuss the comparison of the parameters crucial for IceCube with the ones for short baseline appearance and other experiments.

Atmospheric neutrino disappearance induced by sterile neutrinos
Neutrino oscillation can be described by the following Schrödinger equation in flavor space, where H is the effective Hamiltonian and |ν(L) denotes the flavor state of the neutrino at a distance of L from the source. For the standard 3+1 neutrino framework with 3 active neutrinos (ν e , ν µ , ν τ ) plus one sterile neutrino ν s , including matter effects, the effective Hamiltonian H has the following matrix form Here E is the neutrino energy, m i and U are the neutrino masses and mixing matrix; G F is the Fermi constant, N e is the electron number density of matter and κ is a ratio defined as where N n is the neutron number density. For anti-neutrinos, we need to replace U → U * and G F N e → −G F N e in eq. (2.2). For constant matter density, H does not vary with L so the solution is simply |ν(L) = e −iHL |ν(0) . But the matter density of Earth [17] varies significantly from 1 g/cm 3 (at the surface) to 13 g/cm 3 (at the inner core), so that to obtain accurate results one has to either numerically solve the full differential equation eq. (2.1), or divide the full neutrino path into many segments with approximated constant densities so that

JHEP01(2016)124
Here H i is the Hamiltonian in the i-th segment and L i is the corresponding baseline. Actually this is the main method used to compute oscillation probabilities in the GLoBES package [18] (see also [19]). In this work, we will adopt the same method, i.e. eq. (2.4), to compute probabilities. Defining the S-matrix as where α, β are flavor indices (e, µ, τ and s), the survival probability of ν µ is given by For the energy range of IceCube, ν µ and ν µ dominate the atmospheric neutrino flux while ν e and ν e are negligible; 2 thus only P νµ→νµ and P νµ→νµ will be used in this work. The 4 × 4 mixing matrix is where R ij is a 4×4 matrix whose (i, i), (i, j), (j, i) and (j, j) elements are cos θ ij , sin θ ij e −iδ ij , − sin θ ij e iδ ij and cos θ ij , respectively. The other elements are the same as for a 4×4 identity matrix. Actually there are only three independent CP-violating phases relevant for neutrino oscillations and one possible convention is to set δ 34 = δ 23 = δ 12 = 0. 3 Note that in the limit of vanishing atmospheric and solar mass-squared differences there is no CP effect, and since these two mass-squared differences are much smaller than the one corresponding to sterile neutrinos, the effect of CP phases is suppressed. Therefore in our analysis, we will neglect them. If all CP-violating phases are zero, the last column of U , denoted as u 4 , which we will need later in this analysis, has the following form: Let us now discuss matter effects. For eV-scale neutrinos with TeV-scale energies, the first and second terms in eq. (2.2) are comparable, about m 2 4 /2E ∼ √ 2G F N e ∼ 10 −13 eV. Therefore one expects that matter effects have a significant influence on the probability. As usual for matter effects, a resonance can appear for certain values of energy and baseline. For the case under study, this can only happen for the ν µ → ν µ channel, as has been recently studied for instance in ref. [15]. As a result of the resonance, the survival probability can become zero; anti-muon-neutrinos would completely disappear, even if the active-sterile mixing angles are small. Hence matter effects are crucial for the sensitivity of IceCube to light sterile neutrinos.
However, a less noticed point is that though the matter effect contribution to the effective Hamiltonian is large, in some case, ν µ and ν µ may oscillate as if they are in vacuum. The condition for this case is that U s4 is zero or small, or equivalently that θ 34 is π/2 or JHEP01(2016)124 large. We first show this analytically with the single mass-squared-difference approximation in constant density matter and then numerically verify it by taking into account all masssquared differences and also density variation. Considering that ∆m 2 41 ∆m 2 31 , ∆m 2 21 we can neglect the effect of ∆m 2 31 , ∆m 2 21 and set them to zero. Defining we can write H as follows: Here M is a dimensionless matrix, with u 4 being the last column of U , see eq. (2.7). Note that u T 4 above would be u † 4 if u 4 was complex. But since rephasing M by M → QM Q † with Q = diag(e iα 1 , e iα 2 , e iα 3 , e iα 4 ) does not have any physical effect, we can always make u 4 real by such a transformation. This also implies that all CP-violating phases in the mixing matrix are negligible if ∆m 2 31 , ∆m 2

21
are negligible, as it should be. The S-matrix, assuming constant density, is does not affect the probability so it can be ignored. Considering the case U s4 = 0, we can write u 4 as where e 4 , µ 4 and τ 4 are short for U e4 , U µ4 and U τ 4 . Now M is a block-diagonal 4 × 4 matrix: where M 3 is a 3 × 3 matrix, 14) The ν µ survival probability P νµ→νµ = |S µµ | 2 and S µµ = (e −itM 3 ) µµ can be computed as follows: according to the Cayley-Hamilton theorem, e −itM 3 can be written as Here the coefficients s 0 , s 1 and s 2 are functions of the three eigenvalues (see e.g. [23]),

JHEP01(2016)124
Combine all these result we get Note that for a typical matter density of ρ = 6.5 g/cm 3 and ∆m 2 = 1 eV 2 , A is about 0.5E/TeV. So in the energy range from (10 2 ∼ 10 4 ) GeV, it is quite typical for A to be 1 (or close to 1) and the denominator in eq. (2.22) would be 0 (or close to 0). However, in this case eq. (2.22) is still valid and accurate since the coefficient before µ 2 4 e 2 4 will not blow up when A → 1, as one can check directly. Actually the singularity here corresponds to a branch cut singularity and as it has been proved in ref. [23], all branch cut singularities should cancel out in the S-matrix. This is the deeper reason of the good behavior of eq. (2.22) at A → 1. Therefore, the coefficient before µ 2 4 e 2 4 can be regarded as an O(1) number that varies with t (i.e. with L/E).
If we take the vacuum limit A = 0 we obtain and therefore eq. (2.22) can be written as Eq. (2.24) has an important implication. Since e 2 4 = |U e4 | 2 has been constrained by reactor neutrino experiments to be small, typically less than s 2 13 , the difference between S µµ and S vac µµ is small. We thus reach the conclusion that if U s4 = 0, ν µ (ν µ ) will oscillate as if they would propagate in vacuum. Recalling that the resonance of the matter effect is crucial for the constraints on sterile neutrinos, it implies that the value of U s4 or θ 34 is important for the constraints, and eventually on the ability of IceCube data to rule out sterile neutrino hints.
To further verify the importance of θ 34 , we generate numerical plots without any approximation. We set ∆m 2 21 = 7.5 × 10 −5 eV 2 , ∆m 2 31 = 2.4 × 10 −3 eV 2 , θ 23 = 45 • , θ 13 = JHEP01(2016)124 The plot with θ 34 = 10 • on the right panel (i.e. antineutrino survival probability) shows that the probability reaches 0 at lg(E/GeV) 3.4 and cos θ z −0.9 even though all active-sterile mixing angles are small. This is due to the matter effect resonance in the ν µ -channel. When U s4 is reduced, we can see that the resonance becomes weaker and for U s4 = 0 (i.e. θ 34 = 90 • ), the resonance completely disappears, and the result is indistinguishable from the vacuum case. Note that as θ 34 increases the effect at E ν = 10 2 GeV and cos θ z = −1 becomes significant, which was previously pointed out in [15].
Let us recall here that short baseline disappearance results are essentially electron neutrino appearance results, with an oscillation amplitude of sin 2 2θ µe = sin 2 2θ 14 sin 2 θ 24 . (2.25)

JHEP01(2016)124
We will use later the global fit results from ref. [25] for sin 2 2θ µe (see figure 8 therein), which have been obtained by a fit of available appearance and disappearance results. Also used for comparison will be fit results to sin 2 2θ µe from ref. [7] (see figure 4 therein) that includes ν µ → ν e (ν µ → ν e ) appearance results in combination with various ν e (ν e ) and ν µ (ν µ ) disappearance results (excluding the MiniBooNE low energy excess). Another fit result is from ref. [26] (see figure 4 therein). The data used in those fits is not always the same, as is the treatment of the data, so differences arise. However, the analyses of [25] and [26], using very similar data sets, are giving results in approximate agreement with each other. Hence statements regarding ruling out sterile neutrino hints will depend on the fit result one compares to. Less differences arise for fit results of only appearance data, and we will compare to the results from ref. [25] on ν µ → ν e (ν µ → ν e ) appearance data. Let us note that the LSND results are crucial for the hints for sterile neutrinos, excluding them from global fits reduces the significance dramatically [7]. Another bound of interest is from Super-Kamiokande [27], which found |U µ4 | 2 < 0.054 at 99% C.L., though with assuming U e4 = 0 and ∆m 2 41 > 0.1 eV 2 , and a limit on |U e4 | 2 of about 0.09. Finally, we should mention the 90% C.L. constraint θ 34 ≤ 25 • , obtained from an analysis of muon neutrino disappearance in the MINOS experiment [28]. Regarding electron (anti)neutrino disappearance results, severe tension with various appearance results exists [7,25,26], resembling situations in which inconsistent data sets are combined. Anyway, later we will often take the example values θ 14 = 4 • and θ 14 = 10 • , which are compatible with the 3σ ranges of a most recent global appearance and disappearance fit from ref. [7].
Note that from eq. (2.10) one can show that the presence of matter effects will in general make muon survival probabilities depend on U µ4 and U τ 4 , hence on θ 14 , θ 24 and θ 34 [15]. Indeed, assuming U e4 = 0 (this matrix element has little influence on the final result, as we have essentially a two flavor oscillation case) and following the same calculation as the one leading to eq. (2.21), leads to where C 2 = 1 + Aκ 2 − 4(µ 2 4 + τ 2 4 ) + Aκ and C 3 = 2µ 2 4 + 2τ 2 4 − 1. The muonneutrino survival probability is therefore to good precision a function of |U µ4 | 2 /|U τ 4 | 2 = tan 2 θ 24 / sin 2 θ 34 . Again, we see that the comparison of IceCube atmospheric neutrino results with short baseline disappearance experiments depends on θ 34 .

Neutrino event numbers in IceCube
The neutrino event numbers depend on the neutrino energy E and the zenith angle c z ≡ cos θ z . It can be computed via where A eff (Ā eff ) is the effective area of IceCube for ν µ (ν µ ), Φ (Φ) is the flux of ν µ (ν µ ) and P (P ) is the survival probability of ν µ (ν µ ). All three quantities (A eff , Φ, P ) are functions

JHEP01(2016)124
of the neutrino energy E and the zenith angle c z . For the IceCube-59 data, T = 348.1 days [8]. The factor 2π is due to integration over the azimuthal angle. For Φ andΦ, we use the data from ref. [20]. Since the flux has been computed only up to 10 4 GeV while a small part of events in the IceCube-59 data have energies above 10 4 GeV (most events are in the energy range from 10 2 GeV to 10 4 GeV), we need to extrapolate the data to 10 6 GeV to cover the full data. The final result should be insensitive to the extrapolation because the high energy part has little contribution to the total event number.
The effective area can be extracted from [8], where in figure 1 the simulation results for the event numbers as functions of energy and of zenith angle (without neutrino oscillation) are shown for the conventional atmospheric neutrino flux. So A eff Φ can be obtained from that figure and A eff can be extracted, provided that Φ is known. In practice, a more detailed procedure is adopted by us to take into account the difference between ν and ν and the zenith angle dependence of A eff : we assume an energy-dependent ratio of A eff toĀ eff , where the ratio λ(E) can be taken from figure 2 of ref. [29]. We also assume that the dependence of A eff on the zenith angle is mainly due to the detection efficiency of photons generated by the muon tracks, since the cross section σ(E) of neutrinos with nuclei of water molecules should only depend on E. Under this assumption, we have where N is the number of water molecules and f (c z ) is the detection efficiency. So without neutrino oscillation, eq. (3.1) becomes ∂E∂cz dc z , correspondingly. So we can solve for f and g from the two curves to obtain A eff andĀ eff . Note that the detection efficiency function of Cherenkov photons f in eq. (3.3) (which should mainly depend on c z because it is essentially a geometrical effect) may also have weak dependence on energy. However, events in IceCube are not uniformly distributed from 10 2 GeV to 10 6 GeV, but are rather concentrated around 10 3 GeV. Therefore, the effective integration region of In the actual measurement E can only be partially reconstructed from the muon track; only a lower bound on E can be obtained from the truncated energy loss of the muon (see figure 4 in ref. [8]). Besides, the correlation of the true energy E of the neutrino and the truncated energy loss of the muon is very difficult for us to handle. So in this paper we will not use the energy spectrum information of the data and simply integrate over E in eq. (3.1), though this will somewhat reduce the sensitivity on sterile neutrinos. A full analysis involving the energy spectrum information should be implemented by the IceCube collaboration.

JHEP01(2016)124
Hence, we integrate over E and compute the event number in each small c z -bin, where ∆ cz is the width of cos θ z -bins, equal to 1/25 for 25 bins and lg E ≡ log 10 (E/GeV) . (3.6) Note that though reconstructed zenith angles are not true zenith angles of neutrinos, for energies relevant to our analysis and track events, the difference is much smaller than the bin width [8], and thus we will not consider the difference between reconstructed and true zenith angles in our analysis. For later use, we also define the no-oscillation event number N 0 , Although the total event number is overall reduced due to neutrino disappearance, in practice this effect is not useful to constrain sterile neutrino parameters. The reason is that the fluxes Φ andΦ have a large uncertainty, e.g. about 25% at 10 3 GeV, indicated in figure 11 in [20]. Compared to the statistical uncertainty (about 3%) and the systematic uncertainty (about 4% [30]) in each c z -bin, the large uncertainty in the normalization factor implies the flux can be almost freely renormalized.
The main observable effect caused by a sterile neutrino in IceCube is tilting the zenith angle (cos θ z ) distribution of events, i.e. the smaller c z the more the event numbers are suppressed [15]. The existence of sterile neutrinos causes disappearance for atmospheric neutrinos going through Earth. For a very small |c z | the corresponding oscillation baseline is very short hence neutrinos do not have enough time to oscillate before they arrive at the detector. So the survival probability is always very close to 1 if c z is small enough, no matter how large the mixing angles are. 5 For a large |c z |, neutrinos may have propagated over enough baseline to oscillate and thus the survival probability could be low. They may also experience several oscillatory periods before they arrive and the survival probability would be a value between the minimum and 1, depending on energy and zenith angle. This qualitative analysis can be verified from figure 1, where P is always close to 1 at c z 0. For a large |c z |, P can be relatively small or still large (close to 1), depending on E. Note that eq. (3.5) is an integral over E so only the average value (roughly) of P is of importance. In this sense, we can say that the disappearance signal is stronger at larger |c z | and weaker at small |c z |. Therefore, the sterile neutrino signal in the zenith angle distribution is mainly a tilting effect. To obtain a qualitative understanding of the sensitivity on sterile neutrino parameters, we plot in figure 2 for illustration N/N 0 at c z = −0.86 for the various sterile neutrino parameters (∆m 2 , θ 14 , θ 24 , θ 34 ), i.e. the ratio of events for the oscillated and unoscillated case at a large zenith angle. The more this quantity deviates from 1, the more the zenith angle distribution tilts. The upper left ∆m 2plot shows that N/N 0 drops down quickly from 1 when ∆m 2 goes from 0 to 0.1 eV 2 . This

χ 2 -fit and result
We perform now a numerical analysis of the IceCube-59 data from ref. [8], adopting a conventional χ 2 -function where N ob,i is the observed event number in each bin with statistical uncertainty σ 2 stat,i and N th,i is the predicted event number for each bin. Since the events number in each bin is large enough, we can simply take σ stat,i = N ob,i . For the systematic uncertainty we take, according to [30], σ sys,i = 0.04N ob,i . The factor 1 + a in front of N th,i is the normalization factor of the full flux with a large uncertainty, for instance σ a = 15% at 10 2 GeV or 25% JHEP01(2016)124 . Also shown are 99% C.L. fit results from the world's appearance and disappearance data in purple (from Kopp et al. [25]), in yellow (from Giunti et al. [7]) and in the white contour (from Conrad et al. [26]). The blue contour line is the result of an appearance data only fit from [25]. Figure 6. Exclusion bounds for (∆m 2 , sin 2 2θ 24 ) after marginalizing over θ 34 ; θ 14 is fixed to 10 • (left) or 4 • (right). Also shown are 99% C.L. fit results from the world's appearance and disappearance data in purple (from Kopp et al. [25]), in yellow (from Giunti et al. [7]) and in the white contour (from Conrad et al. [26]). The blue contour line is the result of an appearance data only fit from [25].

JHEP01(2016)124
at 10 3 GeV (see figure 11 in [20]). The major constraint on sterile neutrinos comes from the second term in eq. (4.1). For minimization of the χ 2 -function, it is easy to compute the value of a at the minimum analytically: Here d i ≡ σ 2 stat,i + σ 2 sys,i . In practice, we will use eq. (4.2) for the minimization of the χ 2 -function.
Let us first simply fix θ 14 and θ 34 to certain values so χ 2 is a function of θ 24 , ∆m 2 and a. The conventional treatment is then to replace a in the χ 2 -function with a min given by eq. (4.2). After that, χ 2 will be a function only of θ 24 and ∆m 2 . We perform the χ 2 -fit for several cases. For instance, for (θ 14 , θ 34 ) = (10 • , 0 • ) we find the best-fit point is θ 24 = 8.3 • , ∆m 2 = 0.014 eV 2 with χ 2 min = 9.19. This is no significant effect as for θ 24 = 0 • (i.e., no sterile neutrinos) χ 2 = 9.43 is essentially an equivalently good fit to the data. In figure 4 we show the event distributions of some cases where we can see that the best-fit curve (red, dashed) almost overlaps with the green curve (no sterile neutrinos).
We see from figure 5 that the value of θ 34 is very important. As one increases θ 34 from zero, the constraint becomes stronger until θ 34 is large enough. After that, the constraint becomes weaker when θ 34 is increased. This can be understood via figure 2, where in the JHEP01(2016)124 bottom right panel we can see the event number generally drops down (which implies the disappearance of neutrinos) when θ 34 is increased from 0 • . The event number reaches its minimum at about 30 • and then increases because for θ 34 = 90 • , as we have discussed, the signal of sterile neutrinos is weak since there is no matter effect enhancement. θ 14 , however, has little influence on the IceCube result for ∆m 2 and sin 2 2θ 24 . But the angle is important to compare the outcome to the short baseline appearance data. Note that with θ 14 = 0 there would be no short baseline appearance oscillations, see eq. (2.25). The larger θ 14 , the smaller the value of θ 24 necessary to generate the same value of sin 2 θ µe , hence it becomes more difficult for IceCube to rule it out. Marginalizing over θ 34 results in figure 6. We see that for θ 14 = 10 • the appearance signal is partially compatible within 90%, but ruled out at 99.73% if θ 14 = 4 • . Regarding the appearance plus disappearance data, strong dependence on the fit result exists. For θ 14 = 10 • , the signal from [7] is fully consistent with IceCube-59 at 90%, whereas the full regime from [25] and a part of the regime from from [26] are excluded at 90%. At θ 14 = 4 • , the signal from [7] still cannot be ruled out at about 90%, whereas the full regime from [25] and the one from [26] are ruled out at 95%. Such strong dependence on fit results and parameters not appearing in IceCube expressions will also be present in fits of IceCube-86 data, not yet available, and should be taken into account when interpreting the results.

Conclusion
We have performed a χ 2 -fit to IceCube-59 data in order to constrain sterile neutrino mixing parameters. Special emphasis was put on the question on how muon neutrino disappearance compares to the muon neutrino disappearance in short baseline experiments and other sterile neutrino results. We have stressed the different dependence on the three relevant mixing angles and the important role of the very weakly constrained angle θ 34 , which governs the strength of the matter resonance, and is often set to zero in analyses. The value of θ 14 is also crucial. Fixing this angle implies via the short baseline appearance results the value of θ 24 , on which the analysis of IceCube data is most sensitive to. It was demonstrated that only part of the global parameter space is currently constrained by IceCube, though of course important complementary information is provided. Moreover, there is dependence on which fit result one compares the data to.
The comparison of the various oscillation channels relies on several assumptions, in particular the presence of only one sterile neutrino. In this framework we assume unitarity of the 4 × 4 mixing matrix, and the absence of additional interactions that sterile neutrinos could be sensitive to. As even with these assumptions no final solution to the short baseline disappearance and other problems can be given, it shows that several dedicated oscillation experiments are required to fully settle the issue.