1 Introduction

Searching for sterile neutrinos in ice

Soebur Razzaque***Present address: Space Science Division, U.S. Naval Research Laboratory, 4555 Overlook Ave, SW, Washington, DC 20375, USAE-mail: srazzaqu@gmu.edu and A. Yu. SmirnovE-mail: smirnov@ictp.it

(1) College of Science, George Mason University, Fairfax, Virginia 22030, USA

(2) The Abdus Salam International Centre for Theoretical Physics, I-34100 Trieste, Italy


Oscillation interpretation of the results from the LSND, MiniBooNE and some other experiments requires existence of sterile neutrino with mass eV and mixing with the active neutrinos . It has been realized some time ago that existence of such a neutrino affects significantly the fluxes of atmospheric neutrinos in the TeV range which can be tested by the IceCube Neutrino Observatory. In view of the first IceCube data release we have revisited the oscillations of high energy atmospheric neutrinos in the presence of one sterile neutrino. Properties of the oscillation probabilities are studied in details for various mixing schemes both analytically and numerically. The energy spectra and angular distributions of the events have been computed for the simplest mass, and mixing schemes and confronted with the IceCube data. An illustrative statistical analysis of the present data shows that in the mass mixing case the sterile neutrinos with parameters required by LSND/MiniBooNE can be excluded at about level. The mixing scheme, however, can not be ruled out with currently available IceCube data.

1 Introduction

There are several experimental results which could be interpreted as due to oscillations related to existence of sterile neutrinos with mass eV and rather large mixing with or/and . This includes the LSND result [1], the MiniBooNE excess of events in neutrino and antineutrino channels [2], the reactor antineutrino anomaly [3] and the results of the solar calibration experiments [4] (see [5] for recent interpretation). Global analysis of the short-baseline oscillation experiments shows certain consistency of different evidences in the two sterile neutrinos context [6]. Furthermore, the analysis of CMB data indicates an existence of additional radiation in the Universe [7] with sterile neutrino being one of the plausible candidates. The effective number of neutrino species, -5, looks preferable. The bound on from the Big Bang Nucleosynthesis (BBN) has been relaxed recently allowing for 1-2 additional neutrinos with the best fit value above 3 species [8].

At the same time, the recent global cosmological analysis which includes the CMB data, large scale structure and BBN results, shows that existence of new neutrino species does not relax significantly the bound on mass of the sterile neutrino [9]. For one obtains eV or eV which is smaller than the LSND-required value.

It has been observed some time ago that mixing of sterile neutrinos with eV, and therefore eV, strongly affects the atmospheric neutrino fluxes in the energy range 500 GeV-few TeV. In this energy range the MSW resonance in matter of the Earth is realized in the or channel [10]. The resonance enhancement of oscillations leads to appearance of a dip in the energy spectrum and to distortion of the angular dependence of tracking (induced) events. These effects can be studied in the IceCube detector [10]. Later in [11] an extended study of the oscillation probabilities has been performed in the presence of one or two sterile neutrinos. As an experimental test it has been proposed to measure the ratio of the tracking and cascade (induced by ) events.

Recently AMANDA [12] and IceCube [13] have published the first high statistics data on the atmospheric neutrinos in the TeV range. (See also results from SuperKamiokande [14]). In this connection we present both analytical and numerical study of properties of the relevant oscillation probabilities for different mixing schemes. We compute the energy spectra and angular distributions of events in IceCube. Results of these computations are confronted with the IceCube data and bounds on the parameters of sterile neutrinos have been obtained. We show that observational results substantially depend in the mixing scheme.

The paper is organized as follows. In Sec. 2 we describe the simplest mixing scheme for sterile neutrino (the mass mixing) for which dynamics of evolution is reduced to the evolution. We obtain the analytical expressions for the oscillation probabilities and present results of numerical computations of the probabilities. In Sec. 3 we study modifications of the atmospheric neutrino fluxes due to mixing with sterile neutrinos. We compute the number of events for IceCube and confront them with experimental data. In Sec. 4 we perform an illustrative statistical analysis of the data and obtain bounds on the mixing of sterile neutrinos depending on the sterile neutrino mass. In Sec. 5 the oscillation effects are considered in the mixing scheme. We compute the probabilities and zenith angle distributions of the events, and perform the analysis. In Sec. 6 we study dependence of the oscillation effects on the mixing scheme in the leading order approximation (valid at high energies). Conclusions are given in Sec. 7. In the Appendix we present explicit expressions for the probabilities in the constant density case.

2 Oscillation probabilities in the mass mixing scheme

We will consider mixing of four flavors111 can be treated as the state with zero flavor. of neutrinos which mix in four mass eigenstates , . We assume the neutrino mass hierarchy: , since the opposite situation: , with three active neutrinos in the eV range is strongly disfavored by the cosmological data. The mass squared differences equal

as is required by the LSND/MiniBooNE and fixed by the atmospheric neutrino results.

As we will show, for high energies ( GeV) the electron neutrino mixing can be neglected in the first approximation in consideration of the , oscillations. Therefore the system is reduced to mixing of the three flavor states in three mass eigenstates as , where is the mixing matrix.

In this section we will consider the simplest mixing scheme when mixes in the states and with masses and only. In this case

(1)

Here is the usual 2-3 rotation on the angle and is the rotation of the mass states and on the angle . Explicitly,

(2)

The sterile neutrino mixing is characterized by a single new mixing parameter. In what follows we will refer to (2) as to the mass mixing scheme in contrast to the flavor-mixing scheme which will be discussed in Sec. 5. The simplest mixing scheme allows us to reduce dynamics of the evolution to evolution exactly. Other schemes allow to do this only approximately.

According to (2) mixes with the state

(3)

and there is no mixing of with the orthogonal combination:

(4)

Thus,

In the first approximation at high energies the dominant effect is due to oscillations driven by the largest mass splitting, . Therefore the transitions are described by the flavor mixing in the state. The corresponding elements of mixing matrix equal

(5)

The mass squared difference gives sub-leading effects at high energies. But it produces the leading effects at low energies ( TeV).

Consider evolution of this system in the propagation basis defined as

(6)

It is related to the mass basis as , and therefore the evolution equation for reads

(7)

Here , and is the matrix of the potentials. In we have subtracted the matrix proportional to the unit matrix . In this way we factor out the 2-3 mixing from the evolution of neutrino system. (For earlier work on evolution of 3 and more neutrino states in matter, selection of the propagation basis see [15]). For neutrinos in the electrically neutral medium:

where is the total number density of nucleons, is the number density of neutrons and is the number of electrons per nucleon in the medium. In the electrically and isotopically neutral medium . Therefore , where is the difference of potentials for the system. For antineutrinos: . Explicitly the Hamiltonian is given by

(8)

Here again we have subtracted the matrix proportional to the unit matrix .

The MSW-resonance condition reads

and since the resonance is realized in the antineutrino channel. The resonance energy , (see the level crossing scheme in [17]). The state decouples and is not affected by matter. It evolves independently as

(9)

As follows from the form of the Hamiltonian (8) the evolution matrix (matrix of amplitudes) in the propagation basis can be written as

(10)

(no transitions). From unitarity of we have:

(11)

According to (3) and (4) the states of the propagation basis are related to the flavor states as

(12)

Therefore the matrix in the flavor basis is

Using (12) and (10) we obtain

(13)

Moduli squared of the elements of this matrix give the corresponding oscillation probabilities. In particular, the survival probability, , equals

(14)

and the other oscillation probabilities with participation of are

(15)

These probabilities satisfy the unitarity condition: . Notice that when , the transition probability . In this case . Then for maximal 2-3 mixing we have and , and correspondingly, .

Consider properties of the survival probabilities in the neutrino and antineutrino channels. Using (9) and (14) we obtain

(16)

where is the survival probability. For antineutrinos we have . For TeV the 2-3 phase is small, , so that in the lowest-order approximation

(17)

However, the phase can not be neglected at low energies, TeV. Explicit analytic expression for the amplitude is given in the Appendix. In the absence of mixing with sterile neutrino one has , , and consequently (16) is reduced to usual vacuum oscillation probability due to the 2-3 mixing and 2-3 mass splitting.

Figure 1: The survival probability of the muon antineutrinos (resonance channel) as function of the neutrino energy for different values of the zenith angle () and oscillation parameters (, ).
Figure 2: The same as in Fig. 1 for muon neutrinos.

In Fig. 1 we show the probability as a function of neutrino energy for different values of the zenith angle () and the oscillation parameters. (In our computations we use the PREM model for the Earth density profile [18].) The typical energy-dependent feature of is the resonance dip in the range determined by the resonance energies in the core and in the mantle. For there is a single dip at TeV which corresponds to the MSW resonance in the mantle of the Earth. For (core crossing trajectories) the dependence of the probability on is more complicated. The dip between the resonance energies in the core and mantle is due to the parametric enhancement of oscillations, i.e. due to an interplay between the oscillation effects in three layers with nearly constant density (mantle-core-mantle) [19]. The width of this dip is larger than the width of the MSW dip in constant density medium. There is also the parameteric enhancement of the oscillations at energies above the resonance energy in the mantle [19].

For the (non-resonance) channel, the peaks are absent (see Fig. 2), but another feature related to the matter effect is realized: enhanced transition at low ( TeV) energies. The survival probability decreases with energy in contrast to the channel where increases with energy. The reason can be understood from consideration in the case of constant density (Appendix). At energies below 0.5 TeV the oscillations induced by the 2-3 mixing and mass splitting become important. The dependence of probabilities on energy is given by the oscillatory curve with low frequency in the energy scale and the depth (see analytic expression in (51)). This curve is modulated by high frequency oscillations driven by with small depth. At low energies the phase of the low frequency oscillations is given (see (52) in the Appendix) by

where in last expression the first term is due to the matter effect; the plus sign corresponds to neutrinos and the minus sign to antineutrinos (according to (48), ). In the energy interval TeV the two contributions are comparable. Thus, the matter effect produces an opposite change of the phase velocity: increasing the velocity in the neutrino channel and decreasing it in the antineutrino channel. Consequently the oscillations due to the 2-3 mass splitting and 2-3 mixing develop in the channel at higher energies. Notice that the phase shift is proportional to , and at low energies (see Fig. 2, the upper panel).

Let us consider more general situation when mixes also in the state. We introduce an additional rotation in the - subspace, so that the propagation basis becomes: . Explicitly the mixing matrix in the propagation basis becomes

(18)

Here , , etc.. Now mixes in all three mass eigenstates:

The Hamiltonian in the propagation basis equals , and it can be represented as

(19)

where is the Hamiltonian without rotation (8). The correction is proportional to a small quantity which produces even smaller (suppressed by ) phase than considered in the simplest case above. (The matrix in (19) is symmetric and elements denoted by dots equal to the corresponding transponent elements.) So, the effects of mixing in can be neglected in the first approximation.

The mixing matrix in the flavor basis is given by . The elements of this matrix which describe oscillations with large mass split (dominant at high energies) are the same as in our simplest mixing case (5). They do not depend on .

In what follows we present predictions for IceCube in the simplest mass mixing case. Consideration of the flavor mixing schemes is given in Sec. 5, where we show that, in fact, the probabilities and observables substantially depend on the mixing scheme.

3 Fluxes and numbers of events

The flux at the detector equals

(20)

where and are the original fluxes of and without oscillations. Similar expression holds for the antineutrinos. The effect of oscillations can be neglected (the last equality in (20)). The reason is two fold: at high energies , with ratio for TeV. Furthermore, the transition probability and can be mostly converted to .

Let us consider oscillations in some details. At high energies the mixing of and is strongly suppressed: , where GeV is the resonance energy associated to the “solar” mass splitting . The mixing is absent in the limit , but if non-zero, the 1-3 mixing in matter is also suppressed in the TeV energy range as , where GeV is the energy of 1-3 resonance. Consider the whole scheme with admixture, , in the state . Since for the potential we have in the isotopically neutral medium, the level crossing is in the neutrino channel. The corresponding resonance energy . The depth of oscillations driven by equals

(21)

where and are the mixing parameters in matter. In vacuum: . The mixing and the depth can be enhanced in resonances. In the resonance the mixing is enhanced, , whereas the mixing is suppressed: . As a result, . In the resonance, inversely, the mixing is suppressed , and mixing is enhanced . So that the depth of oscillations equals . Therefore , and the contribution of the original flux to flux at a detector, , is smaller than .

The rate of events in a detector such as IceCube is given by

(22)

with the appropriate integrations over the neutrino energy and solid angle. Additional contribution to the muon events comes from the oscillations, producing a flux at the detector. The tau lepton from interaction has probability to decay into muon, which is then recorded as a event. The energy, however, needs to be times higher than the energy to produce muon tracks of the same energy in the detector. Notice that in the mass mixing scheme ’s appear in the oscillation dip, but this will lead to additional events at low energies. In other mixing schemes ’s are transformed mainly into ’s, and production of is suppressed.

In (22) and are the effective areas of the detector for and . They are given by the effective volume from which the events (muons) are collected with an efficiency of detection as

Here is the number density of nucleons in the surrounding medium and is the neutrino-nucleon charge-current cross section. In turn, is determined by the geometry of the detector and the muon range : . The range can be estimated as

where GeV m,  m, is the mean inelasticity and is the minimum muon energy for detection. At low energies and at TeV the linear increase of changes to the logarithmic one (see [13] for details). Since, usually the data are presented using energy bins of equal size in the scale, the relevant quantity which determines the number of events in a given energy bin is (where the originates from the Jacobian). The differential neutrino flux decreases as , and therefore at low energies increases as . It reaches maximum at TeV and then decreases since has only logarithmic increase. The median energy interval TeV is determined by a condition . This interval includes the region of dips in the oscillation probability and therefore IceCube is well optimized to search for sterile neutrinos with eV. The described dependence of on energy allows one to understand various features of the predicted effects.

The effective area is also given by

where is the geometrical area of the detector, is the survival probability of neutrino passing through the Earth at a given trajectory and is the charged current neutrino-nucleon interaction probability in the vicinity of the detector. The survival probability is given by

where is the length of the trajectory, is the matter density at a distance along the trajectory and is the total neutrino cross-section. For TeV, . The interaction probability is given by

where is the average muon range in the medium and is the Avogadro’s number.

Figure 3: The energy spectrum of integrated over the zenith angles in the intervals and with and without oscillations to sterile neutrinos versus IceCube result. We use the mass mixing scheme with and eV. In the top panel, the error bars denoted by the dashed lines include both statistical and systematic errors as reported by Icecube (Table II in [13]). In the lower panels, statistical-only error bars are denoted by solid vertical lines. For the top panel, the statistical error bars are about the same size as the points.
Figure 4: The energy spectrum integrated over the zenith angle with and without oscillations to sterile neutrinos versus the IceCube results. We use and eV.

In Figs. 3 and 4 we show the sum of the and energy spectra integrated over the solid angle for the mass mixing scheme. An estimation of the size of the oscillation effects is rather easy: maximal, , effect is for in the resonance range; summation with (whose flux is about 1.4 times larger) gives effect; averaging over the zenith angle from to produces another factor , and therefore one arrives at the maximal suppression in the dip. Relative effect increases with narrowing the integration region around vertical direction (see Fig. 5). Now the maximal effect can reach and further enhancement would require experimental separation of neutrino and antineutrino signals. With increase of the dip shifts to high energies as . Increase of the size of the dip with is more complicated. Suppression effect extends to low energies due to oscillations in the channel driven by the 2-3 mixing.

Figure 5: The energy spectrum integrated over different intervals of zenith angles for different values of and .

We also compare the predicted neutrino energy spectra in Figs. 3 and 4 with the “unfolded” energy spectra reconstructed by IceCube [13]. Presently, this comparison can be used for illustration only since reconstruction of the unfolded spectra implies significant smearing and in general is not sensitive to the spectral distortion in small energy intervals. Notice, however, that the size of the dip in the energy scale is larger than the size of the bin of the reconstructed spectrum. To have better sensitivity to the distortion one can further decrease the size of the bin.

According to the Fig. 22 of [13] the statistical error in the relevant energy range is about which is substantially smaller than the size of the dip. Continued operation of IceCube in future will reduce this error further. Large errors are due to systematics: mostly due to uncertainties in the total normalization and tilt of the spectrum. To a large extent they can be eliminated when searching for the dip. Indeed, the systematics has smooth dependence on energy, the systematic errors in different bins are strongly correlate. One can parametrize these uncertainties by a few parameters and determine them by fitting data.

The problem of smearing does not exist in the case of the zenith angle distribution, since muons nearly follow neutrinos, and the zenith angle resolution is . We compute the number of events in a given zenith angle bin using (22) and performing integration from the threshold :

(23)

We then define the suppression factors in the individual bins as

(24)

where are the numbers of events without oscillations which correspond to in (23). In Figs. 6 and 7 we show the zenith angle dependence of the suppression factor for different values of the mixing parameter and (this corresponds to ) and two different thresholds GeV (Fig. 6) and TeV (Fig. 7). Oscillations lead to distortion of the zenith angle distribution. For nearly horizontal direction the effect is mainly due to vacuum oscillations which have enough baseline to develope if TeV. In this case the averaged oscillation effect is given by in agreement with the results of Figs. 6 and  7. The matter effect increases with . According to these figures substantial differences between the energy-integrated distribution with and without sterile mixing are expected in the bins near the vertical direction. For the effect is about and the statistical errors, , are much smaller. For other mixing schemes the distortion can be different. In particular, in the flavor mixing scheme maximal suppression is in the bins (see Sec. 5).

For vertical directions the evaluation of the suppression (integrated over the energy) can be done using the survival probabilities of Figs. 1 and 2. If e.g. , the probabilities averaged over the median energy interval in the neutrino and antineutrino channels are and respectively. Then averaging the contributions of the neutrinos and antineutrinos we obtain , in agreement with results in Figs. 6 and 7. With increase of threshold, the effect of vacuum oscillations in nearly horizontal directions becomes smaller. The effect in the channel increases, whereas in the channel it decreases, thus compensating the overall change.

Figure 6: The zenith angle dependence of the suppression factor of the muon events integrated over the energy from TeV. We use the mass mixing scheme.
Figure 7: The same as in fig. 6 for TeV.
Figure 8: The zenith angle distribution of muons from interactions integrated over the energy with and without oscillations (solid black histograms) to sterile neutrinos. We have renormalized the event distribution according to the best-fit normalization and tilt parameters from the fit (Table 1). Also shown are the IceCube results.

In Fig. 8 we confront the experimental results with the predicted zenith angle distributions computed as

where is taken according to the IceCube simulation (see Fig. 19 from [13]). We have implemented an overall normalization and tilt of the distribution, as we discuss below in (25).

4 Bounds on parameters of sterile neutrinos

To get an idea of the sensitivity of the currently available IceCube data to the sterile neutrino mixing we have performed a fit of the IceCube zenith angle distribution. For a given “model” of mixing characterized by we compute the expected number of muon events in the zenith angle bin . For this we use the IceCube simulation, [13]:

(25)

where is an overall normalization parameter and is a zenith angle tilt parameter. The model without mixing is recovered when . We compare the expected numbers with data and the is defined as

(26)

The variance is calculated by adding in quadrature the statistical and systematic uncertainties as given by IceCube [13]. For our analysis we use the IceCube data in the range of zenith angles (i.e., bins -18), leaving out the last two near horizontal bins where the detector response is not well-understood and contamination of the atmospheric muons is possible. For fixed values we minimize the varying the parameters. The difference

quantifies the rejection significance of the mixing model with respect to the model without mixing.

In Table 1 we show results of our statistical analysis which is reduced to determination of the minimal values of and for given and . We show and the best fit values of and for the case of statistical errors for the individual bins only (see Fig. 8). Also shown is the fit for the “null” hypothesis. Notice that the mixing fits the data better than the model without mixing (“null” model) . Also notice that for , and are below and then they quickly increase with reaching for .

 (eV)
0.005 14.09 0.991 0.0175
0.01 15.29 0.997 0.0086
0.5 0.02 16.50 1.008 -0.0082
0.04 18.88 1.032 -0.0394
0.08 31.73 1.085 -0.1238
0.005 14.56 0.991 0.0217
0.01 15.33 0.997 0.0126
1.0 0.02 16.97 1.010 -0.0052
0.04 20.19 1.033 -0.0347
0.08 39.41 1.092 -0.1344
0.005 14.40 0.991 0.0247
0.01 14.45 0.996 0.0184
2.0 0.02 16.11 1.008 0.0043
0.04 21.87 1.034 -0.0323
0.08 43.29 1.094 -0.1298
0.005 14.03 0.991 0.0246
0.01 14.92 0.996 0.0166
3.0 0.02 15.84 1.008 0.0079
0.04 18.66 1.033 -0.0217
0.08 41.98 1.098 -0.1387
IceCube sim. 14.16 0.982 0.04024
Table 1: Results of the analysis of the IceCube zenith angle distribution. Shown are as well as the best fit values of the normalization parameter and tilt for given values of and in the mass-mixing scheme.

Fig. 9 (left panel) shows the bounds on the sterile neutrino mixing as function of from the analysis which takes into account statistical uncertainty in each bin as well as the systematic uncertainties due to overal normalization and tilt of the zenith angle distribution. These are the main uncertainties. To illustrate possible effect of other systematics we have taken the extreme case: uncorrelated errors for individual bins (see Fig. 9, right panel). (Although it is expected that other possible uncertainties are smooth functions of the zenith angle and therefore correlate in different bins.) In reality the effect of additional errors should be smaller than that. The parameter space to the right hand side from the lines in Fig. 9 is excluded at the indicated confidence level.

The bounds weakly depend on the , as can be seen from the behavior of the suppression factors (Figs. 6 and 7). The bounds are slightly weaker for smaller since in this case the resonance dip shifts from the energy range where IceCube has the highest sensitivity.

We find that with statistical uncertainties only (Fig. 9 left panel) the upper bound is or at level and eV. At level the bounds are and . At the same time interpretation of the LSND/MiniBooNE results in terms of oscillations in the presence of sterile neutrinos requires for eV, and for eV.

With uncorrelated systematic errors (Fig. 9 right panel) the limits become substantially weaker: , is excluded at C.L. only.

Figure 9: Bounds on the active - sterile mixing angle as function of obtained from the IceCube zenith angle distribution of events. Left panel - statistical errors only, right panel: statistical plus uncorrelated systematic errors in each bin.

5 Oscillation effects in the mixing scheme

Let us consider the mixing only, i.e., the simplest scheme of flavor mixing. The corresponding mixing matrix in the flavor basis equals

(27)

where , etc.. Formally it differs from the mixing in (1) by permutation of with the mixing matrix.

Now the mixing matrix elements, which determine the oscillations with splitting , equal , , . They are reduced to the elements of our simplest case (5), if formally we take and and . Therefore in the leading order approximation for high energies the probabilities can be obtained from the probabilities in mass mixing case by taking . In particular, according to (14)

(28)

For we obtain , whereas in the mass mixing scheme this value gives the minimum of the dip .

It is possible to find relation between the sizes of dips for different mixing schemes. For maximal 2-3 mixing we have from (14) the survival probability in the mass mixing scheme (2):

(29)

In the resonance, the amplitude is approximately real. This can be seen using explicit results for the constant density case. Indeed, according to (48) in resonance , and therefore (49) gives . Then from (29) and (28) we obtain relation between the probabilities:

(30)

Our numerical results in Fig. 15 confirm this relation.

Let us consider corrections to the leading order result due to oscillations driven by the 2-3 mixing and splitting. They are sub-dominant at high energies, but become dominant at low energies. In the flavor mixing case it is convenient to consider oscillations immediately in the flavor basis, i.e. take the flavor basis as the propagation one. Using the mixing matrix (27) we find the Hamiltonian of evolution which can be represented in the following form

(31)

At high energies the evolution is described by the first term of the Hamiltonian (which does not depend on the 2-3 mixing), decouples and the corresponding matrix in the flavor basis can be written as

(32)

So that the survival probability, , is in accordance with (28). Indeed, the first term of the Hamiltonian (31) coincides with the Hamiltonian (8) up to permutation of the 2-3 lines, 2-3 columns and substitution , and therefore in this approximation. With the sub-leading term of the Hamiltonian taken into account, evolution is not reduced to the evolution.

Effect of the 2-3 mixing at low energies ( TeV) can be estimated in the following way. In the basis defined in such a way that the Hamiltonian is given by

or explicitly

(33)

For energies much below the sterile resonance, , one can perform a block diagonalization thus decoupling the heaviest state, or simply neglect the 1-3 terms in the Hamiltonian (33). The latter is equivalent to an approximation of negligible matter effect on the angle . So, the evolution is reduced to problem. Similarly to our consideration in Sec. 2 we find (returning to the flavor basis) that the survival probability equals

(34)

where