Matter vs Vacuum oscillations in Atmospheric Neutrinos

Atmospheric neutrinos travel very long distances through earth matter. It is expected that the matter effects lead to significant changes in the neutrino survival and oscillation probabilities. Initial analysis of atmospheric neutrino data by the Super- Kamiokande collaboration is done using the vacuum oscillation hypothesis, which provided a good fit to the data. In this work, we did a study to differentiate the effects of vacuum oscillations and matter modified oscillations in the atmospheric neutrino data. We find that magnetized iron detector, ICAL at INO, can make a 3 sigma discrimination between vacuum oscillations and matter oscillations, for both normal and inverted hierarchies, in ten years.

and a CP violating phase δ CP . The data from the CHOOZ reactor neutrino experiment leads to the strong constraint θ 13 ≤ 10 • [19,20]. The smallness of this mixing angle implies that the solar and the atmospheric neutrino anomalies can be analyzed as independent problems within three flavour oscillation framework [20]. It also leads to the identification θ sol θ 12 and θ atm θ 23 .
The solar neutrinos, produced at the core of the sun, undergo forward elastic scattering as they travel through the solar matter. This scattering leads to matter effect [21][22][23], which modifies the solar electron neutrino survival probability (P ee ). Super-Kamiokande [12] and SNO [13] have measured P ee as a function of neutrino energy for E > 5 MeV and found it to be of a constant value ≈ 0.3. SNO has also measured [13] the neutral current interaction rate of solar neutrinos to be consistent with predictions of the standard solar model [24].
These results, together with CHOOZ bound on θ 13 , provide a 5 σ evidence for the matter effects in solar neutrinos [25]. The measurements of the 71 Ga experiments [9][10][11] imply that P ee > 0.5 for neutrino energies E < 0.5 MeV. This increase in P ee at lower energies can be explained only if ∆ 21 is positive.
The up going atmospheric neutrinos travel thousands of kilometers through earth, during which they undergo forward elastic scattering with earth matter. They also experience the matter effect [21,22] which modifies their of survival and oscillation probabilities. As in the case of the solar neutrinos, this modification depends on the sign of the mass-square difference which, in this case, is ∆ 31 . Given the different magnitudes of • ∆ 21 and ∆ 31 , • the energies of solar and atmospheric neutrinos and • the solar and earth matter densities, the matter modification of atmospheric neutrino probabilities are of a different mathematical form compared to their solar neutrino counterparts [26,27]. An observation of these matter modifications can establish the sign of ∆ 31 . The case of positive ∆ 31 is labelled normal hierarchy (NH) and that of negative ∆ 31 inverted hierarchy (IH). A number of studies [28][29][30][31] considered matter modified oscillations and explored the sensitivity of future atmospheric neutrino detectors to determine whether hierarchy is NH or IH. Recently, Super-Kamiokande experiment analyzed their data using the hypothesis of matter modified oscillations. Their results indicate that the vacuum oscillations are disfavoured at 1.6 σ only [32]. They prefer NH and disfavour IH at 1.7 σ.
At present, the most precise determination of |∆ 31 | and sin 2 2θ 23 comes from the muon neutrino disappearance data of the accelerator neutrino experiments, such as MINOS [33], T2K [34] and NOνA [35]. For baselines less than 1000 km, the matter effects lead to negligibly small changes in ν µ /ν µ survival probabilities [37,38]. Thus, the ν µ /ν µ disappearance data of accelerator neutrino experiments lead to essentially the same values of |∆ 31 | and sin 2 2θ 23 for the three cases: (a) vacuum oscillations, (b) matter oscillations with NH and (c) matter oscillations with IH. It is important to develop methods to make a distinction between these three hypotheses. Without such a distinction, it will be impossible to measure the CP violating phase δ CP in neutrino oscillations because matter effects mimic CP violation [36,39]. In the present work, we explore how this distinction can be made with future atmospheric neutrino data. Disentangling the changes induced in the oscillation probabilities by the matter effects and by δ CP is non-trivial in general. In an atmospheric neutrino detector, the interaction rate of ν µ (ν µ ) depends on both the survival probability P µµ (Pμμ) and the oscillation probability P eµ (Pēμ). However, it has been shown that the sensitivity of this rate to matter effects does not depend on the value of δ CP [29]. Hence, this data can lead to a distinction between the three hypotheses, independent of the value of δ CP .
The ν e /ν e appearance data in long-baseline accelerator neutrino experiments is sensitive to matter effects [40,41]. A precise measurement of the oscillation probabilities, P µe and Pμē, in principle, can make a distinction among the three possibilities. However, this data is also sensitive to δ CP which at present is poorly determined. A given measured values of P µe and Pμē can have three solutions [42,43]: • NH matter oscillations with δ 2 CP and • IH matter oscillations with δ 3 CP .
For each type of oscillating hypothesis, the value of δ CP determined turns out to be different.
Since the determination of δ CP is one of the important goals of future long-baseline accelerator neutrino experiments, making a distinction between the three hypotheses through independent data is crucial.
At present, the data of T2K and NOνA are analyzed using the matter modified oscillation hypothesis with both NH and IH. T2K prefers NH and disfavours IH at 2 σ [44]. It prefers δ CP close to −π/2 for both NH and IH [45]. NOνA also prefers NH but it disfavours IH only at 1 σ [46]. In the case of NH, NOνA allows the full range of δ CP within 1 σ, though it prefers −π/2 for IH. Very recently, NOνA updated its results on hierarchy and δ CP [47]. It still prefers NH with a wide allowed range of δ CP but also has a large 1 σ allowed region for IH around δ CP −π/2. A combined analysis of the latest data of T2K and NOνA [48] mildly prefers IH over NH (with a ∆χ 2 = 1.8) whereas the combined analysis of T2K and NOνA data along with that of Super-Kamiokande shows a mild preference for NH over IH (with a ∆χ 2 = 2.2). The combined data of T2K and NOνA do not show any discrimination between vacuum and matter modified oscillations with NH [49].

II. VACUUM VS. MATTER MODIFIED OSCILLATIONS
The mixing between the neutrino flavor eigenstates and the mass eigenstates is given by where U is a 3 × 3 unitary PMNS matrix. It is parameterized as For neutrino propagation in vacuum, the oscillation probabilities depend on the six parameters: the two mass-squared differences, ∆ 21 = m 2 2 − m 2 1 and ∆ 31 = m 2 3 − m 2 1 , the three mixing angles and δ CP . At present, ∆ 21 , |∆ 31 |, θ 12 and θ 13 are measured quite precisely. In case of the third mixing angle, sin 2 2θ 23 is measured to be close to 1 but sin 2 θ 23 has a rather large range of (0.4 − 0.64). As mentioned in the introduction, the sign of ∆ 31 is not known at present.
The effect of neutrino propagation in matter is parameterized by the Wolfenstein matter term A = 0.76 × 10 −4 ρ (in gm/cc) E(in GeV) [21,22]. Inclusion of this matter term in neutrino evolution induces a change in the mass-square differences and the mixing angles and hence in the probabilities. In this work, we study how this change can be utilized to make a distinction between vacuum and matter modified oscillations. This change depends on not only the matter term but also on the sign of ∆ 31 . Hence we consider both positive and negative values of ∆ 31 .
We first study the difference between matter and vacuum oscillation probabilities for two representative path-lengths for atmospheric neutrinos, L = 5000 km and L = 8000 km. Fig:1 shows the plots of neutrino and anti-neutrino oscillation probabilities P eµ and Pēμ and survival probabilities P µµ and Pμμ for vacuum oscillations as well as for matter modified oscillations with NH. The matter modified probabilities are calculated numerically using the code nuCraft [50], which uses the earth density profile of the PREM model [51]. From the expressions for the oscillation probabilities, it can be shown that the probabilities for IH can be obtained via the relations P µµ (IH) = Pμμ(N H), Pμμ(IH) = P µµ (N H), P µe (IH) = Pμē(N H) and Pēμ(IH) = P eµ (N H) [38]. From Fig:1, we see that, the matter effects increase the peak value of P µe from 0.05 to 0.2 (0.5) for L = 5000 (8000) km and these peak values occur at nearly the same energy where P µµ also peaks. Conservation of probability implies that either P µµ or P µτ should decrease significantly. The maxima of P µµ generally coincide with the minima of P µτ and it can be shown that the change in the value of P µτ near its minimum is very small [38,52]. Hence, most of the reduction in the probability occurs in P µµ . Thus P µµ for matter oscillations with NH is lower than P µµ for vacuum oscillations over a wide range of energies and path-lengths. But Pμμ is essentially the same for both the cases. In the case of IH, the situation is reversed. Therefore, to study the difference of vacuum oscillations from matter modified oscillations of either sign, it is important to measure neutrino and anti-neutrino event rates separately. In this work, we study the sensitivity of Iron Calorimeter (ICAL) at the India-based Neutrino Observatory (INO) to make a distinction between vacuum and matter modified oscillations using atmospheric neutrino data. The charge identification capability of ICAL leads to a very good sensitivity for this distinction [53].

III. METHODOLOGY
The atmospheric neutrinos consist of ν µ ,ν µ , ν e andν e . The ICAL at INO is a 50 kt magnetized iron calorimeter whose iron plates are interspersed with the active detector We have used NUANCE event generator [55] to simulate the atmospheric neutrino events used in this study. It generates neutrino events using atmospheric neutrino fluxes and the relevant cross sections. For a generated event, NUANCE gives the information on the particle ID and the momenta of all interacting particles. The information of the final state particles is given as an input to a GEANT4 simulator of ICAL. This simulator mimics the response of ICAL and generates the electronic signals of the detector in the form of a hit bank information as the output. A reconstruction program sifts through the hit bank information of each event and tries to reconstruct a track. Electrons and positrons in the final state produce a shower and quickly lose their energy. Identifying such particles and reconstructing their energy is an extremely difficult problem. Muons, being minimum ionizing particles, pass through many layers of iron, leaving behind localized hits in the RPCs. Using this hit information, the track of the muon can be reconstructed. Because of the magnetic field, this track will be curved and the bending of the track is opposite for negative and positive muons. Thus ICAL can distinguish between the CC interactions of ν µ andν µ . If a track is reconstructed, the event is considered to be a CC interaction of ν µ /ν µ . The charge, the momentum and the initial direction (cos θ track ), of a reconstructed track, are also calculated from the track properties [56].
We have generated un-oscillated atmospheric neutrino events for 500 years of exposure, using NUANCE. In generating these events, the neutrino fluxes at Kamioka are used as input along with ICAL geometry. The ν µ /ν µ CC events are given as input to GEANT4 and the GEANT4 output is processed by the reconstruction code. Events for which one or more tracks are reconstructed are stored along with the charge, the momentum and the initial direction of the track with the largest momentum. These track variables will be used later in the analysis to bin the events. In the case of ν e /ν e CC events, the electron/positron are redefined to be µ − /µ + and the events are processed through GEANT4 and the reconstruction code. Once again the charge, the momentum and the initial direction of the track with the largest momentum are stored. This redefinition of ν e /ν e CC events is done so that the events which undergo ν e (ν e ) → ν µ (ν µ ) oscillation are properly taken into account in our analysis.
We used accept/reject method on the un-oscillated sample to obtain the oscillated event sample. We calculated the vacuum oscillation probabilities, P µµ , Pμμ, P eµ and Pēμ, using the formula for three flavor oscillations. The corresponding matter modified probabilities, for both signs of ∆ 31 , are calculated numerically using the code nuCraft [50]. The accept/reject method is applied to ν µ (ν µ ) CC events using P µµ (Pμμ) to obtain the muon events due to the survival of ν µ /ν µ . The same method is applied to ν e (ν e ) CC events using P eµ (Pēμ) to obtain the muon events due to the oscillation of ν e /ν e . The track information for each of the selected events is taken from the simulation described in the previous paragraph.
We first do our calculation with the input value δ CP = 0. We later show that our sensitivity to matter effects does not depend on the input value of δ CP . The generated sample is We consider the vacuum oscillations as a hypothesis to be tested against the data samples In Fig:2 ij as a function of track momentum. Fig:3 gives the plots of ∆N µ − j = Σ i ∆N µ − ij and ∆N µ + j = Σ i ∆N µ + ij as a function of track direction. We see that ∆N µ ∓ i , in general, are positive because the matter effects suppress the peak values of P µµ . We also see that the magnitude of ∆N µ − i is larger than that of ∆N µ + i for ∆ 31 positive. The situation is reversed when ∆ 31 is negative. This, of course, is a reflection of the fact that the matter effects are more pronounced in P µµ for positive value of ∆ 31 and in Pμμ for negative values of ∆ 31 . A similar pattern is also seen in ∆N µ ∓ j , for the same reasons.
Track momentum (in GeV) where we introduced three systematic errors π k ij (k = 1, 2, 3), each with its pull parameter ξ k . The first of these is the systematic error in flux normalization which is independent of track momentum and track direction. The second systematic error depends on the track momentum and the third one depends on track direction. The method of calculation of the last two systematic errors is described in the appendix. We computed the χ 2 between the data and the test event samples by where and priors on the pull parameters ξ 2 k are added. For each test value of sin 2 θ 23 , the minimum value of χ 2 is computed by varying the pull parameters ξ k in the range (−3, 3) in steps of 0.1. We obtain the χ 2 for a ten year exposure by dividing the minimum χ 2 by 50. This χ 2 is a measure of ICAL sensitivity to distinguish vacuum oscillations from matter oscillations with ∆ 31 positive.
In calculating the vacuum oscillation probabilities, we held ∆ 21 , ∆ 31 , θ 12 , θ 13 and δ CP fixed. However, each of these parameters has an associated uncertainty. Varying ∆ 21 and θ 12 has very little effect on the atmospheric neutrino oscillation probabilities. So these parameters are kept fixed for the whole calculation. The variation in the two parameters, ∆ 31 and sin 2 θ 13 , can lead to a noticeable change in the probabilities and hence in the event rates. Therefore, we varied these two parameters and marginalized over them. The χ 2 is computed with the addition of the following priors where σ sin 2 θ 13 = 0.00066 and σ ∆ 31 = 0.033 × 10 −3 eV 2 [54]. We found that the minimum χ 2 range, these differences become bigger. And the contribution from all these energy values leads to a significant rise in the χ 2 . The uncertainty in flux normalization, which is common to all the bins, is fairly large. Due to the systematic error in the flux normalization, the χ 2 occurring due to the differences near the P µµ minima regions is drastically reduced whereas the χ 2 occurring due to the large difference near the P µµ maxima regions is relatively unaffected. Hence, the χ 2 with systematic errors is relatively flat with respect to all the test values of sin 2 θ 23 .
We see that χ 2 min = 11.8 for ∆ 31 positive and is 9.5 for ∆ 31 negative. Hence, ICAL can rule out vacuum oscillations at better than 3 σ confidence level, if the matter effects, as prescribed by Wolfenstein [21,22], are present. This sensitivity is there for both the signs of ∆ 31 . The sensitivity obtained is for the input values of sin 2 θ 23 = 0.582 and sin 2 2θ 13 ≈ 0.09, which are the current best-fit values from the analysis of global data [54]. This sensitivity is slightly larger than that obtained for hierarchy discrimination for ICAL in previous studies.
For example, fig 5.7 of ref. [53], gives the hierarchy sensitivity of ICAL for different input values of mixing angles and exposures. The ten-year hierarchy sensitivity, for sin 2 θ 23 = 0.6 and sin 2 2θ 13 = 0.1, is shown to be 11.5, independent of whether the true hierarchy is normal or inverted. It should, however, be noted that the method of analysis used here is very different from that used in previous studies. The kinematic variables used here are track momentum and track direction which are obtained from the reconstruction of the GEANT simulation of atmospheric neutrino events [56], whereas in the previous studies the kinematic variables are the NUANCE output information on muon, smeared with resolution functions.
A detailed comparison of the results from the old method and from the new method is beyond the scope of this paper.
The results presented in Fig:4 assumed δ CP to be zero for both matter and vacuum oscillations. However we should check the sensitivity if the test value of δ CP is varied over its full range (0, 360 • ). We performed this calculation where we kept the true value of δ CP = 0 for matter oscillation probabilities and considered the four test values δ CP = 0, 90 • , 180 • , 270 • for vacuum oscillations. The minimum χ 2 , as a function of test δ CP is shown in Fig:6. We see that this marginalization over δ CP has essentially no effect on the minimum χ 2 . The Recently Super-Kamiokande collaboration performed an analysis of atmospheric neutrino oscillation data along with external constraints [32]. They also searched for evidence of matter effects in this data. They parameterized the matter term in the form α * the standard matter effect and varied α in the range (0, 2). Vacuum oscillations correspond to α = 0 and standard matter oscillations correspond to α = 1. The best fit point occurred for ∆ 31 positive and α = 1. Vacuum oscillations were disfavored with ∆ χ 2 = 5. Negative ∆ 31 , for all values of α, was disfavored with ∆ χ 2 in the range (5 − 6). We did a similar analysis for ICAL also. We generated the matter modified oscillation events or both NH and IH and tested this "data" against the hypothesis of partial matter effect oscillations, as parameterized by Super-Kamiokande. Our results are presented in the Fig: 8, which confirms that the vacuum oscillations, given by α = 0, are ruled out with χ 2 = 11.5 for NH (left panel) and with χ 2 = 9.5 for IH (right panel). We also see that ICAL is very effective in ruling out the wrong hierarchy for all values of α.

V. CONCLUSIONS
Both atmospheric neutrino data and accelerator neutrino data give a consistent picture of neutrino oscillations driven by a mass-squared difference of about 2.5 × 10 −3 eV 2 and a nearly maximal mixing angle. However, the present data can not effectively distinguish between vacuum oscillations and matter modified oscillations. It is important to make this distinction because our ability to determine the CP violating phase δ CP depends on it. In this paper, we have considered the sensitivity of ICAL at INO to make a distinction between these two types of oscillations. We find that a ten year exposure leads to a better than 3 σ distinction, whether ∆ 31 is positive or negative. The difference between the matter and vacuum oscillations is significant for neutrinos if ∆ 31 is positive and for anti-neutrinos if ∆ 31 is negative. Hence, the charge identification capability of ICAL has an important role in giving rise to such good distinguishing ability. This ability is independent of the true value of δ CP . With no charge identification, the discrimination ability is reduced by half.

VI. APPENDIX
In the appendix, we describe the method by which we computed the two systematic errors, one of which depends on track momentum and the other on track direction. Usually, the atmospheric neutrino fluxes are listed with three systematic errors: • Normalization error (common for all the bins) • Tilt error (which depends on neutrino energy) • Direction error (which depends on neutrino direction).
The normalization error on neutrino flux translates directly into the normalization of the neutrino event rates in all bins and hence is independent of track momentum and track direction also. The tilt error is taken to be π tilt ij = π tilt i = 0.05 * ln(E ν /2) and the direction error is assumed to be π nudir ij = π nudir j = 0.05 * | cos(θ ν )| [57], where E ν and θ ν are the energy and the zenith angle of the neutrino respectively.
We outline the procedure for the construction of the transfer matrices which convert the systematic errors in neutrino energy and direction to systematic errors in track momentum and direction. To keep the calculations tractable, we make two simplifying assumptions. We assume that track momentum dependent systematic error depends only on neutrino energy dependent systematic error and not on neutrino direction dependent systematic error. In computing this error, we consider only the modulus of the track momentum and do not make a distinction between the charges. Similarly we will also assume that track direction dependent systematic error depends only on the neutrino direction dependent systematic error and not on the neutrino energy dependent systematic error. Here also, we consider only modulus of the cosine of the zenith angle and do not make a distinction between up-going and down-going particles.
From the simulation of 500 years of un-oscillated events, the ν µ -CC andν µ -CC events with a well reconstructed track are chosen. These events are classified into ten bins of E ν and the event numbers are labeled N 1 , N 2 , ...
The numbers in the column matrix on the LHS, T 1 , T 2 , ..., T 10 are the number of events distributed into ten bins in the modulus of track momentum. Thus, we construct the transfer matrix which converts the distribution of events in neutrino energy into distribution of events in the magnitude of the track momentum.
Similarly we can compute the transfer matrix and π trkdir j for converting the systematic error in neutrino direction into systematic error in | cos θ trk |. Once these two sets of numbers, π trkmm i and π trkdir j are ready, we can write the theoretical number of events (including the systematic errors), binned in |trkmm| and cos(θ trk ) to be N th ij = N 0 ij 1 + 0.2 * ξ norm + π trkmm i ξ trkmm + π trkdir j ξ trkdir , where N 0 ij is the number of theoretically expected events in the ijth bin, according to trkmm and cos(θ trk ). This expression for N th ij is used to calculate χ 2 .