Atmospheric neutrinos, $nu_e-nu_s$ oscillations, and a novel neutrino evolution equation

If a sterile neutrino nu_s with an eV-scale mass and a sizeable mixing to the electron neutrino exists, as indicated by the reactor and gallium neutrino anomalies, a strong resonance enhancement of nu_e-nu_s oscillations of atmospheric neutrinos should occur in the TeV energy range. At these energies neutrino flavour transitions in the 3+1 scheme depend on just one neutrino mass squared difference and are fully described within a 3-flavour oscillation framework. We demonstrate that the flavour transitions of atmospheric nu_e can actually be very accurately described in a 2-flavour framework, with neutrino flavour evolution governed by an inhomogeneous Schroedinger-like equation. Evolution equations of this type have not been previously considered in the theory of neutrino oscillations.


INTRODUCTION
A number of anomalies observed in short baseline neutrino experiments can be explained through the existence of a sterile neutrino ν s with mass in the eV range [1][2][3][4][5][6][7][8][9][10][11][12] (see Refs. [13][14][15][16][17] for reviews). In particular, the so-called reactor neutrino anomaly [4][5][6] and gallium anomaly [7][8][9][10][11][12] indicate that such a neutrino may have a sizeable mixing with electron neutrinos, |U e4 | ∼ 0.1. Though the evidence for the existence of an eV-mass sterile neutrino is not compelling, this is a very exciting possibility, and a considerable amount of work, both theoretical and experimental, has been and is being done to explore it. In particular, a number of experiments focussed on searches of an eV-scale sterile neutrino are currently being performed or planned.
One possible way of looking for a light sterile neutrino is to study the earth matter effects on neutrino oscillations in atmospheric and long baseline accelerator experiments. In particular, while ν µ ↔ ν τ atmospheric neutrino oscillations are in general only very weakly affected by the earth matter, ν µ ↔ ν s oscillations may be strongly influenced by it, leading to characteristic energy and zenith angle dependent distortions of the observed muon neutrino flux (see  for an incomplete list of references). Oscillations between ν e and ν s may also be strongly affected by the earth's matter.
The simplest framework to describe active-sterile neutrino oscillations is the so-called 3+1 scheme, with three light mass-eigenstate neutrinos consisting predominantly of the usual active neutrinos ν e , ν µ , ν τ and a heavier (or much lighter) state ν 4 , consisting mainly of a sterile state ν s . The short baseline anomalies mentioned above together with constraints from cosmology suggest that ν 4 is the heavier state, with a mass m 4 = O(1) eV. In this case in the energy range of a few TeV one can expect resonantly enhancedν µ ↔ν s and ν e ↔ ν s oscillations of atmospheric neutrinos that come to detectors from the lower hemisphere and therefore pass through the matter of the earth. Disappearance of atmosphericν µ due to oscillations into a sterile state driven by ∆m 2 ∼ 1 eV 2 was considered in [34,36,[38][39][40][41][42][43][44], and the recent analysis of the IceCube data [45] has put stringent constraints on the corresponding allowed parameter space. The matter-enhanced ν e -ν s oscillations of atmospheric neutrinos have been considered in [42], where it was pointed out that the available IceCube data posed only rather weak constraints on the parameters governing this oscillation channel.
In the present paper we concentrate on ν e -ν s oscillations of atmospheric neutrinos. As mentioned above, they are expected to be strongly enhanced by the earth's matter in the TeV energy range, making them an interesting candidate for discovery of sterile neutrinos. For baselines limited by the diameter of the earth and TeV-scale neutrtino energies, neutrino flavour transitions in the 3+1 scheme are governed by just one neutrino mass squared difference, ∆m 2 41 , and are fully described within 3-flavour oscillation frameworks [38,43]. We demonstrate that the flavour transitions of atmospheric ν e can actually be very accurately described in a 2-flavour approach, where neutrino oscillations are governed by an inhomogeneous Schrödinger-like equation. Evolution equations of this type have never been previously used for describing neutrino oscillations.
While numerical integration of the full 3-flavour neutrino evolution equation does not in general pose any calculational difficulties, the 2-flavour approach is more advantageous in several aspects. In particular, it admits simple analytical solutions for a number of matter density profiles (such as e.g. constant matter density, several layers of constant densities and an arbitrary density profile in the adiabatic approximation) which are much more transparent and more easily tractable than the corresponding 3-flavour solutions.

FORMALISM
Neutrino oscillations in matter are described in the 3+1 framework by the evolution equation where ν = (ν e , ν µ , ν τ , ν s ) T , U is the leptonic mixing matrix, ∆m 2 ik = m 2 i − m 2 k , G F is the Fermi constant, and the neutrino potentials are with N e and N n being the electron and neutron number densities of matter, respectively.
For the TeV energy range and terrestrial baselines the energy splittings ∆m 2 21 /2E and ∆m 2 31 /2E can be neglected. In this case neutrino oscillations are governed by just one mass squared difference, and CP violating effects in neutrino oscillations are unobservable. One can therefore choose the mixing matrix U to be real. The evolution equation (1) then takes the form From experiment we know that |U µ4 |, |U τ 4 | 0.2. Unitarity then implies 1.
It will be convenient for us to perform a rotation ν µ , ν τ → ν µ , ν τ , that is, to introduce a new basis for the neutrino amplitudes according to where c β ≡ cos β, s β ≡ sin β. The neutrino evolution equation in the new basis reads i(d/dx)ν = V HV † ν . Note that the rotation with V does not affect the matrix of neutrino potentials, and therefore the evolution equation in the primed basis has the same form as eq. (4), except that one has to replace ν µ and ν τ by and U µ4 and U τ 4 in the effective Hamiltonian by, respectively, It is easy to see now that one can reduce the 3+1 neutrino evolution equation to a 3-flavour form [38,43]. Indeed, we can choose the angle β by requiring, e.g., U τ 4 = 0, which gives The condition U τ 4 = 0 results in vanishing the third line and the third column of the effective Hamiltonian in the rotated basis V HV † . This means that the state completely decouples from the rest of the neutrino system and does not evolve, whereas the orthogonal combination of ν µ and ν τ , is mixed with ν e and ν s . Thus, in the limit ∆m 2 21 , ∆m 2 31 → 0 neutrino flavour transitions take place only between ν e , ν µ and ν s . 1 The corresponding 3-flavour evolution equation is where We are interested in the ν e ↔ ν s oscillations inside the earth at the energies close to the MSW resonance energy in this channel. The transitions in the ν µ −ν e and ν µ −ν s channels are then non-resonant, so that the probabilities of ν µ oscillations remain small. However, at the energies of interest the original (unoscillated) flux of ν µ is significantly larger than the original ν e flux (F e ), and the same holds true for the ν µ flux, 2 barring the possibility |U µ4 | |U τ 4 |. Therefore, in calculating the ν e flux at the detector site, one 1 Obviously, had we chosen U µ4 = 0 instead of U τ 4 = 0, the decoupled state would have been ν µ rather than ν τ . However, the physical content of the oscillating and non-oscillating neutrino states in both cases is, of course, the same, and is given by the right-hand sides of eqs. (10) and (9), respectively. 2 By which we mean the weighted sum of the original νµ and ντ fluxes, Note that at the TeV energy scale the flux of atmospheric ντ is still rather small, F cannot in general neglect the transitions from ν µ , unless the mixing parameter U µ4 satisfies |U µ4 | |U e4 |. Still, a significant simplification is possible if one takes into account that, though the transitions ν µ → ν e and ν µ → ν s have to be taken into account because of the large initial ν µ flux, the inverse processes, i.e. back reaction of ν e and ν s on the ν µ flux, can be neglected. This means that the flux of ν µ can be considered as fixed, serving as an external source for the evolution equations for the ν e and ν s amplitudes. The 3-flavour oscillation problem of eqs. (11), (12) can therefore be reduced to a much simpler 2-flavour one.
The argument proceeds as follows. First, since the initial value |ν µ (0)| is much larger than |ν e (0)| and ν s (0) = 0, we retain on the right-hand side of the equation for the amplitude ν µ (x) in eq. (11) only the term H µµ ν µ (x). The equation can then be immediately solved, giving Substituting this into the equations for the amplitudes ν e (x) and ν s (x) yields an inhomogeneous 2-flavour evolution equation: Here the external sources are Eq. (11) should be solved with the initial conditions This would be equivalent to solving the original system (11) with the initial conditions 3 and with the back reaction on ν µ (x) neglected. Such a calculation would have been correct if the initial neutrino state were a coherent superposition of ν e and ν µ , whereas in reality the original unoscillated atmospheric neutrino flux consists of an incoherent sum of the ν e and ν µ fluxes. This can be taken into account by introducing a random phase for ν µ (0) in eq. (17) according to ν µ (0) → e iϕ ν µ (0) and averaging over the random phase ϕ at the probabilities level. Such an averaging can be readily done, as will be discussed below.
Let us now demonstrate how one can find the solution (ν e , ν s ) T of eq. (14) provided that the solution (ν The amplitudes ν e (x) and ν s (x) should be rephased similarly, but this would not affect the oscillation probabilities and therefore we keep for these amplitudes the old (unprimed) notation.
with the initial condition (16).
is actually the evolution matrix of the homogeneous equation. Indeed, it satisfies the equation is the Hamiltonian of eq. (19)) and obeys the initial condition S(0, 0) = 1. Note that the matrix S is unitary, as it should be.
The solution of the full inhomogeneous equation (19) is then readily found: One can check by direct substitution that satisfies eq. (19) with the initial condition (16). This leads to very simple expressions for ν e (x) and ν s (x). In terms of the parameters α and β we have These formulas can be further simplified by noting that the expressions in the square brackets in (23) and (24) are actually the elements of the matrix S(x, x ) (this can also be seen directly from eq. (22)). Indeed, from the well known properties of the evolution matrix it follows that S(x, x ) = S(x, 0)S(0, x ) = S(x, 0)S −1 (x , 0). By making use of eq. (21) one then finds This is our final result (though eqs. (23) and (24) may also be useful). Note that the sources f e (x) and f s (x) have the same coordinate dependence, 5 which simplifies the calculation of ν e (x) and ν s (x). Very simple analytic expressions for the parameters α and β can be obtained, e.g., for constant-density matter or constant-density layers model of the earth density profile [47].
It is also useful to present expressions (23) and (24) for ν e (x) and ν s (x) in a more compact form: where In order to obtain the survival probability of electron neutrinosP  14)), over the random phase ϕ of the initial amplitude ν µ (0). The quantities a(x) and b(x) are proportional to ν µ (0) (because so are f e (x) and f s (x)), and therefore the averaging over the random phase of ν µ (0) is achieved by discarding the contribution to |ν e (x)| 2 coming from the interference of the first term on the right-hand side of eq. (27) with the remaining two terms: 6 If ν e (x) is found from eq. (25), the averaging would correspond to neglecting the interference between the first and the second terms on the right-hand side of this equation when calculating |ν e (x)| 2 . The oscillated ν e flux at the detector site (corresponding to the baseline L) is then The dependence on the original muon neutrino flux F At the same time, when the oscillated neutrino flux is calculated directly from the 3-flavour neutrino evolution, without approximating it by a 2-flavour one, the electron neutrino flux at the detector site is given by where the oscillation probabilities P ee (L) and P µ e (L) should be obtained from the evolution equation (11).

RESULTS
In this section we compare the results of the complete 3-flavour calculations with those obtained in the approximate 2-flavour framework. In both cases constant density layers model was used for the matter 6 Indeed, since the last two terms on the right hand side of eq.  density profile of the earth (see [47] for the numerical values of the parameters used and for the explicit formula for the 2-flavour evolution matrix S). For completeness, we quote here the 2-flavour expressions for the parameters α(x, 0) and β(x, 0) for the three-layer model of the earth's density profile for neutrinos crossing the core of the earth (mantle-core-mantle trajectories) [47]: where Here θ 1 and θ 2 are the in-matter mixing angles in the mantle and in the core of the earth, respectively; s 1 (c 1 ) is the sine (cosine) of the oscillation phase in each mantle layer, while s 2 (c 2 ) is the sine (cosine) of the oscillation phase in earth's core. For mantle-only crossing neutrino trajectories one has to put θ 2 = θ 1 and replace c 1 c 2 − s 1 s 2 and s 1 c 2 + s 2 c 2 by, respectively, the cosine and sine of the total oscillation phase.
In figs. 1 and 2 the solid (red) curves show the results of the full 3-flavour calculations, whereas the dashed (blue) curves represent the results obtained in the 2-flavour approach with averaging over the random phase of the initial ν µ state, as discussed above. We observe a very good agreement between the results of these two approaches -the oscillated ν e fluxes obtained in the 3-flavour and 2-flavour frameworks are almost indistinguishable. With increasing |U µ4 | the difference between the 2-flavour and 3-flavour results becomes more pronounced, compare figs. 1 and 2. This is because the approximation of neglecting the back reaction of ν e and ν s on the ν µ flux becomes less accurate in this case. Also the dip at E 8 TeV becomes less deep due to the increased ν µ → ν e transitions.
For illustration, we also show the ν e flux obtained from the 2-flavour calculation without averaging over the random phase of ν µ (0) (dot-dashed green curves). As expected, this gives a wrong result.  Atmospheric neutrino oscillations in the 3+1 scheme can in principle both enhance and suppress the ν e flux compared to the unoscillated flux F (0) e . Interestingly, our calculations presented in figs. 1 and 2 demonstrate only the reduction of the ν e flux. This happens because in the considered TeV energy region the ν e → ν s disappearance is strongly enhanced by matter effects and dominates over the ν e appearance due to the ν µ → ν e transitions.

SUMMARY AND DISCUSSION
If a sterile neutrino with an eV-scale mass exists and has a sizeable mixing to ν e , as the reactor and gallium neutrino anomalies suggest, the flux of atmospheric ν e in the TeV energy region can be significantly affected by ν e ↔ ν s oscillations. Strong enhancement of the oscillation probability in this channel can occur due to the MSW and parametric resonances of neutrino oscillations in the earth, just like it occurs in the GeV energy range for the usual ν e ↔ ν µ and ν e ↔ ν τ oscillations driven by the atmospheric mass squared difference ∆m 2 31 and the mixing angle θ 13 (see, e.g., [48]). The simplest framework to consider active-sterile neutrino oscillations is the 3+1 scheme. For TeVscale neutrinos and terrestrial baselines, oscillations in this scheme are governed by just one mass squared difference, ∆m 2 41 , whereas the effects of ∆m 2 31 and ∆m 2 21 can be neglected. In this approximation neutrino oscillations in the 3+1 scheme can be reduced to pure 3-flavour oscillations between ν e , ν s , and ν µ , which is a linear combination of ν µ and ν τ . We have demonstrated that in the TeV energy range flavour transitions of atmospheric neutrinos can actually be very well described by a 2-flavour evolution equation. This proved to be possible because the probabilities of oscillations of the muon neutrinos are rather small in this energy region. This does not mean that the oscillations of the atmospheric ν µ can be completely ignored, as their inital flux is about a factor of 20 larger than the ν e flux in the TeV energy region. Therefore, the transitions ν µ → ν e and ν µ → ν s should be taken into account. However, the back reaction of ν e and ν s on the ν µ flux can to a good accuracy be neglected. As a result, the flux of the muon neutrinos remains practically unchanged and serves as an external source for a 2-flavour evolution equation.
It can be seen from figs. 1 and 2 that our effective 2-flavour approach provides a very good approximation for the full 3-flavour evolution of atmospheric neutrinos in the TeV energy range. The main advantage of the 2-flavour framework is that for many matter density profiles of practical interest it allows analytical solutions which are much simpler and much more transparent than the corresponding 3-flavour expressions. The 2-flavour neutrino evolution equation derived here, eq. (14) [or eq. (19)], has a form of an inhomogeneous Schrödinger-like equation (equation with external sources). To the best of the present author's knowledge, evolution equations of this type have never been previously used for describing neutrino oscillations.