Pinning down neutrino oscillation parameters in the 2–3 sector with a magnetised atmospheric neutrino detector: a new study

We determine the sensitivity to neutrino oscillation parameters from a study of atmospheric neutrinos in a magnetised detector such as the ICAL at the proposed India-based Neutrino Observatory. In such a detector, which can separately count νμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\mu $$\end{document} and ν¯μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\nu }_\mu $$\end{document}-induced events, the relatively smaller (about 5%) uncertainties on the neutrino–antineutrino flux ratios translate to a constraint in the χ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2$$\end{document} analysis that results in a significant improvement in the precision with which neutrino oscillation parameters such as sin2θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{23}$$\end{document} can be determined. Such an effect is unique to all magnetisable detectors and constitutes a great advantage in determining neutrino oscillation parameters using such detectors. Such a study has been performed for the first time here. Along with an increase in the kinematic range compared to earlier analyses, this results in sensitivities to oscillation parameters in the 2–3 sector that are comparable to or better than those from accelerator experiments where the fluxes are significantly higher. For example, the 1σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\sigma $$\end{document} precisions on sin2θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{23}$$\end{document} and |Δm32(31)2|\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\Delta {m^2_{32(31)}}|$$\end{document} achievable for 500 kton year exposure of ICAL are ∼9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }9$$\end{document} and ∼2.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }2.5$$\end{document}%, respectively, for both normal and inverted hierarchies. The mass hierarchy sensitivity achievable with this combination when the true hierarchy is normal (inverted) for the same exposure is Δχ2≈8.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \chi ^2\approx 8.5$$\end{document} (Δχ2≈9.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \chi ^2\approx 9.5$$\end{document}).


Introduction
One of the open questions in neutrino physics is the mass ordering of the neutrinos; whether they are ordered normally or inverted. Many experiments intend to determine the mass ordering, of which the 50 kton Iron Calorimeter (ICAL) detector at the proposed India-based Neutrino Observatory (INO) is one ambitious experiment [1]. ICAL will be a magnetised iron calorimeter mainly sensitive to muons a e-mail: slakshmi@prl.res.in produced in the charged-current (CC) interactions of atmospheric muon neutrinos (and antineutrinos) with the iron target in the detector. It can distinguish CC muon-neutrinoinduced events from antineutrino-induced ones since the former interaction produces μ − while the latter produces μ + in the detector and ICAL has excellent muon charge identification (cid) capability. This is also crucial to determine precisely the momentum of the muons through bending in the magnetic field. Since matter effects are different between neutrino and antineutrino propagation in the Earth, this feature can help resolve the neutrino mass ordering by determining the sign of the 2-3 mass-squared difference m 2 32 ≡ m 2 3 − m 2 2 , m i , i = 1, 2, 3, being the neutrino mass eigenstates [2]. In addition, the matter effects improve the sensitivity to the magnitude | m 2 32 | of this mass-squared difference as well as to the 2-3 mixing angle, θ 23 , provided the across-generation mixing angle θ 13 is rather well-known, which is indeed the case [3][4][5][6][7].
Many previous analyses have been reported, projecting the sensitivities of ICAL detector to oscillation parameters in the 2-3 sector [1,[9][10][11] as also the mass ordering. The sensitivity to the mass hierarchy is directly proportional to the value of sin 2 2θ 13 , which is quite precisely known [12][13][14][15][16][17]. It also depends on the ability of ICAL to separate neutrino and antineutrino events which is possible since ICAL is magnetised. While there is an uncertainty of about 20% on the atmospheric neutrino fluxes themselves, the uncertainty on their ratios is much smaller, about 5%, and it was ignored in earlier analyses [9][10][11]. In this paper, we show that this smaller uncertainty on the ratio acts as a constraint that in turn significantly shrinks the allowed parameter space, especially for sin 2 θ 23 . For instance, we will see that the precision on sin 2 θ 23 decreases from 13 to 9% in a certain analysis mode when this constraint is included. This is generally true for all magnetised detectors. To the best of our knowledge such an effect has not been discussed in the literature earlier.
The paper is structured as follows. All the results in this paper are with detailed simulation studies of the physics processes at the ICAL detector. The main steps involved in this are neutrino event generation, inclusion of detector responses and efficiencies, inclusion of oscillations, binning in observables and χ 2 analysis. The procedure of neutrino event generation with the NUANCE neutrino generator and the implementation of oscillations are discussed in detail in Sect. 2. The choice of observables and kinematic regions used in the analysis, along with the inclusion of detector responses are discussed in detail in Sect. 3. The effect of increasing the energy range of observed muons is also explained in this section. The detailed χ 2 analysis and a discussion of the systematic errors that have been considered are presented in Sect. 4. The results of precision measurements and hierarchy sensitivity studies are shown in Sect. 5. The impact of the additional pull in the ν/ν flux ratio implemented in this analysis is discussed in detail in Sect. 6. The summary and conclusions are given in Sect. 7.

Neutrino events generation
The interactions of interest in ICAL are the CC interactions of ν μ and ν μ with the iron target in ICAL. These ν μ (ν μ ) in ICAL come from both ν μ and ν e atmospheric fluxes via ν μ → ν μ and ν e → ν μ oscillations. The first channel gives the number of ν μ events which have survived and the second, subdominant, one gives the number from oscillations of ν e to ν μ . The number of events ICAL sees will be a sum of these events. Thus, where n d is the number of target nucleons in the detector, σ μ is the differential neutrino interaction cross section in terms of the energy and direction of the CC lepton produced, μ and e are the ν μ and ν e fluxes and P αβ is the oscillation probability of ν α → ν β . The number of unoscillated events over an exposure time t in a bin of (E μ , cos θ μ ) is obtained from the NUANCE (version 3.504) neutrino generator [18] using the Honda 3D atmospheric neutrino fluxes [19], neutrino-nucleus cross sections, and a simplified ICAL detector geometry. While NUANCE lists details of all the final state particles including the muon and all hadrons, ICAL will be optimised to determine accurately the energy and direction of the muons (seen as a clean track in the detector) and the summed energy of all the hadrons in the final state (since it cannot distinguish individual hadrons).
Even though the analyses are done for a smaller number of years (say 10), a huge sample of NUANCE events for a very large number of years (here 1000 years) is generated and scaled down to the required number of years during the analysis. This is mainly done to reduce the effect of statistical (Monte Carlo) fluctuations on sensitivity studies, which may alter the results. A detailed discussion about the effect of fluctuations on oscillation sensitivity studies will be discussed in Appendix.
A sample of 1000 years of unoscillated events was generated using NUANCE (version 3.504). Two sets were generated: 1. CC muon events using the μ flux and 2. CC muon events obtained by swapping e ↔ μ fluxes.
This generates the so-called muon-and swapped-muon events that correspond to the two terms in Eq. (1).

Choice of event generator
The NUANCE neutrino generator incorporates the so-called Honda-3D atmospheric neutrino fluxes [19]. Recently, the GENIE [20] neutrino generator was modified to include atmospheric neutrino fluxes and some preliminary analysis has been done by the INO collaboration. In both cases, the Honda-3D fluxes have been used. While the fluxes used by both codes are therefore identical, there are differences in cross sections.
In both generators, the deep inelastic scattering (DIS) cross sections are the same (with parton distribution functions from the standard CERN libraries) and so also the quasi-elastic (QE) cross section. The main difference is that NUANCE includes nearly a hundred resonance (88) (RES) cross sections and so has the most detailed consideration of even highly suppressed resonance channels. The two differ in the implementation of the relativistic Fermi gas model to include final state effects as well as in the transition region between resonance and deep inelastic scattering. Hence, the hadron spectra may differ; however, there should be no difference in the muon spectrum. One limitation of GENIE is that the total exposure time has not been implemented so far and so the total number of events in ICAL are normalised using NUANCE. However, a movement towards using GENIE for ICAL simulations has begun since it is being regularly updated and improved. However, the present analysis used the NUANCE neutrino generator alone. For this reason, a systematic error of 5% on neutrino cross section was included to take care of these effects.

Oscillation probabilities
These events are oscillated depending on the neutrino oscillation parameters being used. The oscillation probabilities are  [21] has been used to model the varying Earth matter densities encountered by the neutrinos during their travel through the Earth. The Runge-Kutta solver method is used to calculate the oscillation probabilities [22] for various energies E ν and distances L, or equivalently, cos θ ν (θ ν being the zenith angle) of the neutrino. Further discussion of the oscillation probabilities and plots of a few sample curves are presented in the next section after listing the kinematical range of interest. The oscillation is applied event by event (for both muon and swapped muon events) as discussed in detail in Ref. [11] and it is a time consuming process since the actual sample contains 1000 years of events. The central values of the oscillations parameters that we have used are from Ref. [11] and are given in Table 1 along with their known 3σ range. The value of θ 13 alone has been updated from Ref. [8]. Note that δ C P is currently unknown and its true value has been assumed to be 0 • for the purposes of this calculation. Furthermore, since ICAL is insensitive [10] to this parameter, it has been kept fixed in the calculation, along with the values of the 1-2 oscillation parameters m 2 21 and sin 2 θ 12 , which also do not affect the results.
It is convenient to define the effective mass-squared difference m 2 eff , which is the measured quantity whose value is related to m 2 31 and m 2 21 as [23,24]: 2 21 (cos 2 θ 12 − cos δ CP sin θ 13 sin 2θ 12 tan θ 23 ). ( When m 2 eff is varied within its 3σ range, the mass-squared differences are determined according to m 2 31 = m 2 eff + m 2 21 (cos 2 θ 12 − cos δ CP sin θ 13 sin 2θ 12 tan θ 23 ); for the normal hierarchy when m 2 eff > 0, with m 2 31 ↔ − m 2 32 for the inverted hierarchy when m 2 eff < 0. A neater definition of the mass ordering can be obtained by defining the quantity Then, switching the ordering from normal to inverted is exactly equivalent to the interchange m 2 ↔ − m 2 , with no change in its magnitude. However, since the marginalisation is to be done on the observed quantity m 2 eff , we use this quantity, but need to keep in mind that m 2 32 ↔ − m 2 31 when the ordering is flipped between NH and IH in this case.

Choice of observables and kinematic regions
The expression in Eq. (1) is for the ideal case when the detector has perfect resolutions and 100% efficiencies. In this analysis, realistic resolutions and efficiencies obtained from GEANT-4-based simulation studies of ICAL [25][26][27][28][29] have been incorporated; this not only reduces the overall events due to the reconstruction efficiency factor but also smears out the final state (or observed) muon energy and direction and that of the hadron energy as well.

ICAL detection efficiencies
Detailed simulations analyses of the reconstruction efficiency, direction and energy resolution of muons in ICAL have been presented in Refs. [25,26]. In addition, the relative cid (charge identification) efficiency of muons (ability of ICAL to distinguish μ − from μ + ) has also been presented here. The detector has good direction reconstruction capability (better than about 1 • for few-GeV muons) and excellent cid efficiency (better than 99% for few-GeV muons) also for muons. The detailed simulation studies of the response of ICAL to hadrons have been presented in Refs. [28,29]. Hadron hits are identified and calibrated to reconstruct the energy of hadrons in neutrino-induced interactions in ICAL. The present analysis has used these results to simulate the observed events in ICAL. Note that the efficiency of an event is taken to be the ability to see a muon, that is, to be able to reconstruct it. Hence when hadron energy is added as the third observable, the efficiency in detecting an event remains the same. At the time this calculation was begun, the responses of both muons and hadrons in the peripheral parts of ICAL was not completely understood. Hence, instead of propagating the NUANCE events through the simulated ICAL detector in GEANT and obtaining a more realistic set of "observed" values of the energy and momentum of the final state particles, the true values of these variables were smeared according to  [25,28].
Studies in the peripheral region Ref. [26] have shown that the direction reconstruction is similar to that in the central region while the muon momentum resolution is somewhat worse with a reduced reconstruction efficiency of about ∼65%. Since the peripheral region comprises about ∼50% of the detector, this leads to an overall decrease in the reconstruction efficiency from ∼85 to ∼75%, with only a slight worsening of the energy resolution. We have not accounted for this difference; this would lead to a corresponding increase in exposure time by a factor 85/75 = 1.13 so that a result with say 10 years of running the experiment will only be achieved in 11.3 years instead.
It should be noted that instead of reconstructing the neutrino energy and direction using the muon and hadron information and then binning in neutrino energy, the analyses have been done by taking all the observables separately. This is because of the poor energy and direction resolution of neutrinos in ICAL detector owing to the fact that they are driven by the responses of the detector to hadrons, which are worse compared to those of muons. Still, the addition of the extra information regarding hadrons improves the sensitivity of ICAL to oscillation parameters, as shown in Ref. [11].

Effect of extending the energy range of observed muons
The first highlight of this paper is widening the energy range of the observables, especially that of observed muons. Since ICAL is optimised for muon detection, it is desirable to make use of all the events available to perform the oscillation analysis. As opposed to all the earlier studies in ICAL [1,[9][10][11] which restricted themselves to performing the analyses in the energy range of only 1-11 GeV of the observed muon energy, the analysis we present here uses the region of the observed muon energy E obs μ = 0.5-25 GeV. It will be seen in Sect. 5 that the inclusion of the higher energy bins beyond the upper limit of 11 GeV, used in earlier studies, improves the results.
The motivation to use the extended range of the observed muon energy is seen in Fig. 1 where the dominant oscillation probability, P μμ , is shown as a function of the zenith angle cos θ for different values of the neutrino energy, E ν ≥ 5 GeV. The large matter effects are clearly visible at lower energies. With increase in energy, the curve smooths out (matter effects become small so that P μμ ∼ P μμ ) and correspondingly P eμ becomes vanishingly small. The features observed around θ ≤ 33 • or cos θ 0.84 occur because the neutrinos just graze the core-mantle boundary at this angle and the matter effects due to the large density change are consequently large. (The figure shows the survival probability P μμ for the normal hierarchy. Similar effects will be seen in the antineutrino survival probability, P μμ , in the case of the inverted hierarchy). Note also the vanishing of P μμ for high energies, E ν 20 GeV, in the upward direction (cos θ → 1).
The sensitivity to m 2 eff is shown for two different energies, E ν = 10, 22 GeV (representing the last energy bins of the previous analysis and the present one, respectively) in Fig. 1 as well. The minimum moves to the left with increasing m 2 eff so that the solid (dashed) line corresponds to m 2 eff = 2.1 (2.6) × 10 −3 eV 2 , which is the presently allowed 3σ range. It can be seen that the position of the minimum of P μμ is more sensitive to the value of m 2 eff at the larger value of energy, although the probability itself is not sensitive to the sign of this quantity at this energy. Hence the inclusion of the higher energy bins improves the sensitivity to these oscillation parameters, as we shall see.
In addition, there is a small fraction of up-going events arising due to ν τ CC interactions in the detector. While P eτ is small, P μτ is driven by θ 23 and can be large, for instance, whenever P μμ is small, as shown in Fig. 1. However, the cross section for ν τ N CC interactions is severely suppressed due to the large tau mass and has a threshold E ν τ 3.5 GeV, where the atmospheric neutrino fluxes are small. In a separate study [30] it has been shown that roughly 230 (90) τ − (τ + ) events are expected in 10 years running at ICAL (assuming 100% detector efficiency), of which 17% decay into muons and thus add to the direct muon events. This works out to about 40 (15) additional events in 10 years, and this will be further decreased when detector effects are included; this is to be compared to the 6600 (2500) up-going direct muon μ − (μ + ) events estimated over the same time period. Furthermore, these events will be distributed roughly equally over the angular range 0.2 ≤ cos θ ≤ 1.0 in the upward direction. This region has been divided into 12 bins in our analysis, resulting in a negligible contamination arising from tau CC events with the subsequent decay of τ → μ in the detector. Hence we are not concerned here with these events.
We note in passing that in contrast, about 66% of these events decay into hadrons and these add to the NC background at high energies. Since their contribution is significant (and only in the upward direction), they can be separated out and studied [30]. These events are not studied in this paper; they are mentioned here for completeness.

The binning scheme
This is similar to, and an extension of, the one used in the earlier analyses [11]. The observables in the analysis are the observed (i.e., smeared) muon energy E obs μ , observed muon direction cos θ obs μ and observed hadron energy E obs had , where the true total hadron energy is defined as [11] E had ≡ E ν − E μ . There are two different analysis sets, one in which only the muon energy and direction, (E obs μ , cos θ obs μ ), are used, called the 2D (mu only) binning scheme and the other in which all the three observables (E obs μ , cos θ obs μ , E obs had ) are used, which is also known as the 3D (or with-hadron) binning scheme. The details of the two binning schemes are shown in Table 2.
It should be noted that in the current analysis the direction cos θ obs μ = +1 is taken as the up direction. Since atmospheric neutrino oscillations are mainly in the up direction, more bins are assigned in this region than in the down direction. The E obs μ bins up to 1-11 GeV are taken to be same as those used in Ref. [11]. A bin of width 0.5 GeV is added in the lower energy range. In the higher range of E obs μ , four bins are added, two of them with bin width each of 1.5 and 2.5 GeV, respectively, and the last two bins of width 5 GeV, thus making the total number of 15 bins in E obs μ . The cos θ obs μ bins and the E obs had are kept the same as in Ref. [11]. Thus same bins are used in the overlapping energy range, and the additional bin sizes were optimised to obtain reasonable event rate as well as sensitivity to oscillation parameters. This gives an optimised result on the whole.
For the bins from 11-20 GeV, we have used the resolutions and efficiencies obtained in Ref. [25]. Although Ref. [25] showed results up to 20 GeV only, results had been obtained in that study up to 25 GeV. We note that the results are very similar between 20 and 25 GeV, as can be seen from the trend of the results in this paper.
For the low energy muons below 1 GeV, a new analysis, again in the central region, was performed [27] and the results from that analysis were appended to the existing look-up table to obtain a look-up table from 0.5-25 GeV which was used in the analysis presented in our paper. For example, the muon momentum resolution varies from 24-16% (20-10%) for input muon momenta from 0.5-1.0 GeV and cos θ = 0.35 (0.95), respectively. The reconstruction efficiency is 2-23% (78-90%) for these two zenith angles while the relative cid efficiency is roughly constant at 62% (72%) respectively.
As mentioned above, the number of hadron bins was retained as before, as also the energy range. No extension of hadron energies beyond E obs had = 15 GeV was used, since this gave only a marginal improvement in χ 2 while ICAL's sensitivity to hadrons at higher energies in terms of the number of hits in the detector tends to saturate [28].

Number of events
The true number of oscillated events is given by where P αβ is the oscillation probability of a flavour α to a flavour β, N 0 μ ± and N 0 e ± refer to unoscillated muon and swapped-muon events generated by NUANCE arising from survived atmospheric ν μ → ν μ and oscillated atmospheric ν e → ν μ interactions in ICAL corresponding to the two terms in Eq. (1).
The number of events per bin including the charge misidentified ones is given as where N tot μ − (N tot μ + ) is the total number of oscillated ν μ (ν μ ) CC muon-neutrino events observed in the bin (E obs μ , cos θ obs μ ). The quantity rec is the reconstruction efficiency of muons with a given energy and direction and cid is the relative charge identification efficiency of the same. The reconstruction and charge identification efficiencies for μ − and μ + have been taken to be the same; studies show [25] that they are only marginally different in a few energy-cos θ bins. Finally, the events in a bin are considered non-zero if there is at least one event in that bin. Now this 1000 years sample, oscillated according to the central values of the oscillation parameters listed in Table 1, is scaled to the required number of years to generate the "data". The current precision analysis is done for 10 years of exposure of 50 kton ICAL (500 kton year). In order to generate the "theory" for comparison with "data" for the χ 2 analysis, the oscillation parameters are changed within their 3σ ranges and the aforementioned processes are repeated. Different theories are generated by changing the oscillation parameters.

χ 2 analysis
Systematic uncertainties play a very important role in determining the sensitivity to oscillation parameters in any experiment. The inclusion of these uncertainties always gives a worse χ 2 than the one obtained when we have no uncertainties at all.
In this new analysis, an extra systematic uncertainty compared to the older analyses is included: the uncertainty on the neutrino-antineutrino flux ratio. This has been considered for the first time in such an analysis and will be seen to have a great impact because ICAL is a magnetised detector that can separate μ − and μ + events. With the inclusion of this uncertainty, the χ 2 can no longer be expressed as a sum of the separate contributions of neutrino and antineutrino events. When the systematic errors are implemented using the usual method of pulls [31][32][33][34][35] we have where i, j, k sum over muon energy, muon angle, and hadron energy bins. The number of theory (expected) events in each bin, with systematic errors, is given by where T 0± i j (k) is the corresponding number of events without systematic errors, D ± i j (k) is the number of "data" (observed) events in each bin, and ξ l ± are the pulls corresponding to the same five systematic uncertainties, l = 1, . . . , 5, for each of neutrino and antineutrino contributions, as considered in the earlier analyses by the INO collaboration [1,[9][10][11].
The five systematic uncertainties include the flux normalisation, shape (spectral or energy dependence) uncertainty or "tilt" and zenith angle uncertainties, the cross-section uncertainties, and an overall systematic uncertainty due to the detector response (see Ref. [11] for details).
The cross-section uncertainty is assumed to be process independent. At high energies, the cross section is dominated by deep inelastic scattering (DIS) where the uncertainties are smaller; however, in the energy range of interest here, all processes [quasi-elastic (QE), resonance (RES) and DIS] have significant contributions and so a common cross-section uncertainty is used.
In Eq. 6 ξ 6 is the 11th (additional) pull and π 6 is taken to be 2.5%. The effect of the new pull can be understood by considering its contribution alone on the ratio of neutrino to antineutrino events: This pull therefore accounts for the uncertainty in the flux ratio; 2π 6 corresponds to the 1σ error (when ξ 6 = 1); this gives the 1σ error on the ratio to be 5%, consistent with Refs. [19,36]. In the earlier analysis with 10 pulls only, the pulls for N − and N + were independent so that they could be in the same or opposite directions. The introduction of the 11th pull constrains the ratio and results in a (negative) correlation between the normalisations of the T + and T − events, as will become clear from the discussions presented in Sects. 5 and 6. Without this pull, the total χ 2 can be simply expressed as a sum over the μ − and μ + contributions: Note that the observed muon events have contributions from both the μ and the e fluxes; here we have assumed the same systematic error on both the μ / μ and the e / e ratios (although, in principle, this can be included separately).
We have also ignored the small differences due to additional possible uncertainty in the μ / e flux ratios since the contribution from the second term in Eq. (1), i.e., from ν e → ν μ oscillations, is subdominant/small. An 8% prior at 1σ is also added on sin 2 2θ 13 , since this quantity is known to this accuracy [4,5]. No prior is imposed on θ 23 and m 2 32 , since the precision measurements of these parameters are to be carried out with ICAL. The contribution to χ 2 due to prior is defined as where σ (sin 2 2θ 13 ) = 0.08 × sin 2 θ true 13 . Thus the total χ 2 is defined as: where χ 2 corresponds suitably to χ 2 11 (new analysis) or χ 2 10 (repeat of the older analysis with extended energy range). During χ 2 minimisation, χ 2 ICAL is first minimised with respect to the pull variables ξ l for a given set of oscillation parameters, then marginalised over the ranges of the oscillation parameters sin 2 θ 23 , m 2 eff and sin 2 2θ 13 given in Table 1. The third column of the table shows the 3σ range over which the parameter values are varied. These along with the best-fit values of θ 12 and m 2 21 are obtained from the global fits in Refs. [37][38][39][40][41]. As mentioned earlier, the parameter δ C P is kept fixed at zero throughout this analysis. The relative precision achieved on a parameter λ (here λ being sin 2 θ 23 or | m 2 eff |) at 1σ is expressed as where λ max-2σ and λ min-2σ are the maximum and minimum allowed values of λ at 2σ ; λ true is the true choice. The statistical significance of the obtained result is denoted by nσ , where n = χ 2 , which is given by χ 2 0 being the minimum value of χ 2 ICAL in the allowed parameter range. With no statistical fluctuations, χ 2 0 = 0.

Results: precision measurement of parameters
The precision measurement of oscillation parameters in the atmospheric sector in the energy range 0.5-25 GeV is probed in the current analysis using 500 kton year exposure of ICAL detector. Comparisons with previous analyses in the [1][2][3][4][5][6][7][8][9][10][11] GeV energy range are also done to illustrate how much better is the sensitivity with our new analysis. Sensitivities to the oscillation parameters θ 23 and | m 2 eff | are found separately, when the other parameter and θ 13 are marginalised over their 3σ ranges. Marginalisation is also done over the two possible mass hierarchies; since atmospheric neutrino events in ICAL are sensitive to the mass hierarchy (also discussed below), the wrong hierarchy always gives a worse value of χ 2 during marginalisation. Typically, the normal hierarchy (NH) is taken to be the true hierarchy (a couple of results with inverted hierarchy (IH) are shown for completeness) and 500 kton years of exposure is used (10 years of running the experiment).

Precision measurement of sin 2 θ 23
The relative 1σ precision on sin 2 θ 23 obtained from different analyses, with the normal hierarchy as the true hierarchy, is shown in Fig. 2 for different cases which are the combinations of energy ranges, binning schemes and number of pulls.
The other parameters | m 2 eff | and θ 13 have been marginalised over their 3σ ranges as given in Table 1. Percentage precisions on sin 2 θ 23 at 1σ obtained with different analyses are shown in Table 3.
2D Analysis of sin 2 θ 23 : The extension of the E obs μ range to 0.5-25 GeV improves the precision to 13% from 14% obtained in the earlier analysis with 1-11 GeV [9,11]. A very large enhancement in precision to 9% is obtained when the 11th pull is included. This is a very significant observation for  all magnetised detectors and the reason for this improvement is discussed in Sect. 6. 3D Analysis of sin 2 θ 23 : A similar result of 9% precision is obtained with the 3D analysis when the 11th pull is included. This is much better than the earlier result of 13% [11]. The 2D muon-only analysis with 11 pulls gives a comparable result; this is significant for the following reason: A ν μ CC event typically has a muon and several hadrons produced at the common interaction vertex, propagating in the detector. In the 2D analysis only the muon momentum and direction information are used. This is obtained from the Kalman filter algorithm which uses the bending of the track in the presence of the magnetic field to determine both these quantities. The technique is efficient even if a few hadron hits are misidentified as muon hits, or even if there are no hits at all in a layer, since it is able to extrapolate the expected track across such "empty" layers. The clean muon track typically extends beyond the hadron hits (which rarely traverse more than 6-7 layers in ICAL) and this enables precise reconstruction of the muon momentum and direction. Hence, the reconstruction is more robust with respect to mis-id of muon hits as hadrons and vice versa.
On the other hand, the only way hadron energy can be calibrated in ICAL is by counting the number of hits in the shower. Here, any mis-id of hadron hits as muon hits or vice versa will alter the estimate of the hadron energy. The hadron energy reconstruction is thus more susceptible to such misidentification. Hence the most robust results are obtained from 2D analysis and therefore it is advantageous to have as good a 2D analysis as possible. Of course an improved measurement of hadron energy can further improve this result.
The earlier best result (with 1-11 GeV in muon energy, hadron bins as listed in Table 2, and including only 10 pulls) [11] is shown in comparison with the best results of the  Fig. 3 (left) for NH. Indeed, it can be seen that the earlier best result to 3σ precision is worse than that already known from other experiments as listed in Table 1. The precision obtainable with IH as the true hierarchy was also studied with our latest analysis. The relative 1σ precision obtained with IH is 9.25% which is comparable to that obtained from NH (8.9%), as can be seen in Fig. 3 (right). This is comparable to or slightly better than the precision on sin 2 θ 23 obtained by T2K [42,43]. That is, the 11th pull acts as a constraint on the relative μ + and μ − events so that a magnetised atmospheric neutrino detector can achieve the precision obtained by an accelerator experiment. We believe that this fact has been pointed out for the first time in this paper.

Precision on
Since ICAL is a magnetised iron calorimeter, it can measure | m 2 32 | with very good precision. As in the case of sin 2 θ 23 , there are six different analyses which give the results as shown in Fig. 4. The percentage precisions at 1σ obtained for the magnitude of m 2 32 ( m 2 31 ) are shown in Table 4 when the true hierarchy is normal (inverted). It can be seen  Left comparison of χ 2 ICAL obtained from 2D and 3D analyses in the energy range E obs μ = 0.5-25 GeV and with 11 pulls with the normal hierarchy as the true one. The best result is for the 3D anal-ysis (which is similar to either 10 or 11 pulls. Right best sensitivity to | m 2 32 | in the current analysis (3D, 11 pulls) for both normal and inverted hierarchies that ICAL will be able to determine | m 2 32 | (| m 2 31 |) with a greater precision than sin 2 θ 23 , in all energy ranges.
The current best results with 2D and 3D analyses are shown in Fig. 5 (left). The new 2D and 3D analyses in the range E obs μ = 0.5-25 GeV constitute an improvement over the older 2D and 3D analyses [11]. However, as the precision on | m 2 32 | achievable by ICAL is already quite good, the improvement does not seem to be as pronounced as in the case of sin 2 θ 23 although the 3D analysis is itself an improvement over the current bounds.
2D Analysis of | m 2 32 |: The muon-only (2D) analysis with 10 pulls gives a precision of 4%, which is a ∼22% improvement over the old value. With the additional constraint of the 11th pull in the E obs μ = 0.5-25 GeV case, the precision achievable is similar, which shows that the new pull does not improve the precision further. This is in contrast to the precision measurement of sin 2 θ 23 , which is impacted mainly by the constraint on the ν μ -ν μ flux ratio. The reason for this is discussed in Sect. 6.
3D Analysis of | m 2 32 |: The precision obtained with 3D binning and 10 pulls in 0.5-25 GeV improves to 2.5% from the older value of 3%, which corresponds to a ∼17% increase. The addition of the 11th pull again does not improve the precision further. Precision measurement with IH as the true hierarchy gives a precision which is comparable to that obtained for NH. This is shown in the right panel of with NH. This is also the best-fit point since the data are scaled down from 1000 to 10 years and so fluctuations are suppressed. Right comparison of the ICAL projected 10 years sensitivity with the current T2K result at 90% CL. The ICAL contour is shifted to match the best-fit value of | m 2 32 | from T2K [44] | m 2 32 | (| m 2 31 |). We now discuss the parameter space allowed by our latest analysis. The analysis was done for the 11 pull case with E obs μ = 0.5-25 GeV and with hadrons, for 500 kton years of ICAL exposure. The normal hierarchy is assumed to be the true hierarchy.
The 90 and 99% confidence contours for 500 kton year exposure of ICAL for NH with the true choices of (sin 2 θ 23 , | m 2 32 |) = (0.5, 2.366 × 10 −3 eV 2 )) are shown in the left panel of Fig. 6. Similar results hold true for the inverted hierarchy as the true hierarchy. It can be seen that the extension of the energy range for analysis and constraining the ν μ -ν μ flux ratio in ICAL result in an improved sensitivity to the precision on both parameters. The projected ICAL precision on sin 2 θ 23 is better than or comparable to the current T2K precision as can be seen from the right panel of Fig. 6, where the 90% CL contour is compared to the results from T2K [44]. A comparison of these projected results for ICAL (NH) with the current results from MINOS [45] and T2K [44] are shown in Fig. 7. It must be remembered, though, that these experiments are already taking data while ICAL is yet to be constructed! 5.4 Sensitivity to neutrino mass ordering ICAL with its magnetisability is an exclusive mass hierarchy machine. Most importantly, its ability to discriminate the normal and inverted mass ordering is independent of the CP phase [1]. The ability of ICAL to distinguish the correct mass ordering in the 2-3 sector is given by where χ 2 false-MO is the minimum χ 2 calculated using the false ordering, while allowing θ 23 , θ 13 and | m 2 eff | to vary over the 3σ range given in Table 1  The addition of hadron energy as the third observable increases the χ 2 MO-ICAL to ∼8.5, for an exposure time of 10 years, an improvement over the earlier value of 7.7. Again the addition of the 11th pull has no effect on the hierarchy sensitivity when normal ordering is taken to be the true ordering. Inverted ordering: While trends are similar to that with normal order as the true order, the inclusion of the 11th pull has significant impact on the mass hierarchy sensitivity when the true ordering is inverted. In fact, the best sensitivity (E obs μ = 0.5-25 GeV, with 11 pulls and 3D binning) is better with the inverted than with the normal hierarchy (by ∼16%). The reason for this effect with the 11th pull is discussed in the next section.
It should be noted that the values of χ 2 MO-ICAL are lower than those reported in Ref. [11], for the same exposure time. This is due to the fact that the earlier analysis used the input value θ 13 ∼ 9.217 • while our analysis uses the current 1 bestfit value [8] θ 13 ∼ 8.729 • . Given that the best fit of θ 13 has reduced further to ∼8.5 • [37, 46] it is even more important 1 The best fit at the time when this analysis was begun.
to perform the analysis in as wide a kinematic range (energy and direction) as possible to get the best possible hierarchy discrimination.

Dependence on θ 23
The sensitivity to mass ordering is also known to depend on the true value of θ 23 ; it is higher for larger θ 23 as can be seen from Fig. 9. While ICAL has better sensitivity to the inverted ordering when θ 23 is in the first octant, the reverse is true when it is in the second octant. This is due to the dominance of the P μμ term in Eq. (1) arising from the survived ν μ s compared to the P eμ term from oscillated ν e s, and the nature of its dependence on θ 23 . In fact, the results would also be vastly improved by "switching off" the P eμ term since their dependence on the oscillation parameters (especially θ 23 ) are practically the opposite of each other [22]. This is not possible with atmospheric neutrinos where Nature provides both flavours, but is possible with neutrino beams. In the latter case, however, the fact that there is a single base-line which results in a significant dependence on and correlation with the CP violating phase δ CP complicates the analysis, as with MINOS [45], T2K [42], LBNE [47,48] or NOνA [10,49].
To summarise, the sensitivity to the mass ordering in the 2-3 sector improves with the addition of higher energy bins in the analysis, while constraining the ν μ -ν μ flux ratio improves the sensitivity only when the true ordering is inverted. ICAL's ability to determine this mass ordering is significant owing to its magnetisability and its 50 kton mass. Improvement in energy resolutions will further improve the  Fig. 9 The best case values of χ 2 MO-ICAL as a function of exposure time in years. It can be seen that the sensitivity to the mass hierarchy increases linearly with the true value of θ 23 when it is varied over the 3σ range of the parameter detector's sensitivity to this parameter. Also, in ICAL, the mass ordering can be determined independent of the CP violating phase δ CP because of the range of baselines involved in atmospheric neutrinos [1]; this is in contrast to beam/shortbase-line experiments where there is a non-trivial sensitivity to the 2-3 mass ordering depending on the true value of the CP phase. Hence ICAL will be important in the determination of this mass ordering, and any amount of improvement in determining this parameter is noteworthy.

Impact of the 11th pull on determination of the oscillation parameters
The 11th pull accounts for the fact that the ratios of the ν μ and ν μ fluxes are better known (to within 5%) than the absolute fluxes themselves [19]. This is implemented by using a pull π 6 = 0.025, which contributes with the opposite sign for neutrino and antineutrino events, as can be seen from Eq. (9). It is seen that the inclusion of the 11th pull is most visible in the determination of θ 23 , which becomes more constrained when this pull is included. One way to understand this is to re-bin the events in a single variable L/E (of the final state muon) and consider the effect of just this pull. Figure 10 shows the effect of θ 23 on both the neutrino and the antineutrino events. The thick solid line is the "data" corresponding to θ 23 = 45 • while the thin solid "theory" line corresponds to a fit with θ 23 = 40 • and without any pull. In both cases, reducing θ 23 from the true value increases the event rate in every bin (the opposite will hold with the inverted hierarchy; here the normal hierarchy is shown). Note that the downgoing events are not shown here. The dot-dashed and dashed curves corresponding to the labels 11-and 10-pulls show the effect of changing the (single) normalisation of the theory with and without the 11th pull. The overall normalisation of the events in the 10-pull case can be independently varied for neutrinos and antineutrinos (decreased by 4 and 3% in figure) to improve the agreement of the 40 • theory line to the data, resulting in smaller overall χ 2 in this fit. On the other hand, a 2.5% decrease in the normalisation of neutrino-induced events in the 11-pull case is accompanied by a 2.5% increase in the antineutrino case, so that the agreement with the neutrino data becomes better, but that with the antineutrino data becomes worse.
(Of course, it can be applied vice versa, but the smaller χ 2 is obtained with this choice since there are about twice as many neutrino events as antineutrino ones due to the smaller cross section of the latter.) Hence it is not possible to improve the agreement of the 40 • theory with data by tuning the normalisation in the analysis when including this pull; this results in a larger χ 2 compared to the 10-pull case where the normalisations of neutrino-and antineutrino-induced events can be independently varied. This gives rise to the tighter constraints on θ 23 when the 11th pull is added. We note that only a detector like ICAL that is capable of charge identification can successfully implement this pull as a constraint.
It can also be seen from Fig. 10 that there is greater sensitivity to θ 23 in the neutrino rather than in the antineutrino sector. The reverse is true with inverted mass ordering. In the determination of the sensitivity to the mass ordering, the minimum χ 2 with the false ordering is found. That means, when the true ordering is inverted, that the "theory" is obtained using the normal ordering, where there is greater sensitivity to θ 23 as just mentioned. When the 11th pull is included, therefore, the discrepancy between theory and data cannot be achieved by changing θ 23 and hence there is more sensitivity to determination of the mass ordering when the true ordering is inverted, provided θ 23 is in the first octant. The reverse is true when θ 23 is in the second octant, as can be seen from Fig. 9.
In addition, it can be seen from Fig. 10 that there is sensitivity to oscillation parameters near and beyond log 10 (L(km)/E(GeV)) ∼ 4. This is precisely the region that is included when the range of E obs μ is extended from 1 GeV down to 0.5 GeV. Similarly, although smaller, sensitivity to the parameters is also seen for log 10 (L(km)/E(GeV)) ∼ 2-3, which corresponds to the extension in the higher energy end from 11 to 25 GeV.

Summary and discussion
This paper contains a simulation study of the physics potential of the proposed 50 kton magnetised Iron Calorimeter detector (ICAL) at INO which aims to probe neutrino oscillation parameters by observing atmospheric neutrino oscillations and studying their Earth matter effects as they propagate through the Earth. This will be done by detecting (mainly) the charged-current interactions of ν μ and ν μ in the detector by means of the final state muons. The detector, which is optimised for the detection of muons in the GeV energy range, will have a magnetic field which will enable the distinction of ν μ and ν μ events by identifying the charge of the muon in the final state, thus making ICAL an excellent detector to determine the neutrino mass ordering. Not only this, the magnetic field helps to improve the precision measurement of the mixing angle θ 23 (and | m 2 32 | and the mass ordering as well).
The main themes of our study were the effects of extending the observed muon energy range to 0.5-25 GeV (from 1 to 11 GeV used in earlier studies) and that of a constraint on the ν μν μ flux ratio on the sensitivity of a 50 kton ICAL to neutrino oscillation parameters in the 2-3 sector. The second-and the aspect which is found to have the biggest impact so far on oscillation sensitivity-arises from two facts; one that ICAL detects atmospheric neutrinos and the other that this massive detector will be magnetised.
In particular, we show that the relatively small uncertainties on the atmospheric neutrino-antineutrino flux ratios act as a constraint in analyses where the neutrino and antineutrino events can be separated. This is true in a magnetised detector such as ICAL where the magnetic field distinguishes muons and antimuons produced in charged-current interactions of neutrinos and antineutrinos respectively with the detector material. The presence of the magnetic field enables ICAL to distinguish between neutrino and antineutrino events; hence, inclusion of the uncertainty on the ν μ /ν μ flux ratio as an additional pull translates to a constraint on this ratio which in turn significantly improves the precision with which sin 2 θ 23 can be determined. Such a constraint is applicable for all magnetised detectors which have good charge identification capabilities and the impact of this constraint has been shown for the first time in this paper. As a consequence, our simulation studies show that, with 10 years of data taking, ICAL will not only be able to determine | m 2 32 | with good precision (as expected) but can also pin down sin 2 θ 23 to a precision better than the current limits set by T2K.
The precision that can be achieved on these parameters in about 10 years' running is about 2.5 and 9% for | m 2 32 | and sin 2 θ 23 , respectively; see details in Tables 3 and 4.
The studies presented here assume that there is perfect separation between different types of charged-current (CC) and neutral current (NC) events. The event separation efficiency will affect the results of the analysis since they will determine the actual number of events in each bin apart from the contamination from NC oscillation-independent events. However, the inclusion of this consideration is beyond the scope of this paper. Preliminary studies [50] show that certain selection criteria can be applied so that an event sample which comprises more than 95% CC muon events can be obtained. The criterion results in cutting out events with E ν 0.5 GeV; this has determined largely the range (lower limit) of muon energy analysed in this paper. Neutral current (NC) events with a hard hadron leaving a track in the detector can mimic muon CC events. According to a preliminary study, the ratio of the neutral current events to the charged-current muon events satisfying the cuts for track reconstruction is about ∼1-2% in the bins of reconstructed track momentum from 1-25 GeV/c. In the bin 0.5-1 GeV/c, this is about 5%. This can be further reduced by including the information on hit multiplicity per layer, but this study is not yet complete. Thus the NC contamination to the CC ν μ sample is around 1-2% in most of the bins. The contamination in the lower bin is larger owing to the fact that the muons satisfying the track criteria are lower in number themselves and hence the contamination is larger in that bin.
Also, improvements in the reconstruction efficiencies and resolutions (which have been studied with GEANT-4 simulations of the detector) as well as possible changes in detector geometry can all alter the results of this analysis. Note that the current studies were all done with the atmospheric neutrino fluxes computed at the Super Kamiokande site [19]. The fluxes at Theni where ICAL is proposed to be built are slightly different and are smaller at energies less than 10 GeV [19,51]; this will also impact the analysis.
Even within these limitations, it appears that the physics results of ICAL will have the capability to impact global fits to neutrino data and thus any new analysis will open a window to understanding the neutrino oscillation parameters better on the whole and the atmospheric neutrino fluxes themselves, in particular. Fig. 11 The value of χ 2 ICAL vs. the number of years of the sample theory to be scaled to 10 years for θ 23 = 39 • (left) and | m 2 32 | = 2.466×10 −3 eV 2 (right). It can be seen that the scaling of only 100 years of the sample size gives too high a χ 2 , which is more prominent in the case of m 2 32 . As the number of years of the sample to be scaled increases, the value of χ 2 also decreases as well as saturates for large enough sample sizes the theory itself. The χ 2 stabilises at a point when the sample is fairly large and this is the number of years of exposure to be taken and scaled down for the analysis. It can be seen that use of 1000 years sample yields fairly stable results and thus the analyses have made use of such a sample of charged-current muon-neutrino events. We note, however, that the results are much more sensitive to the sample size in the case of | m 2 32 | than for θ 23 .