A simulation study of tau neutrino events at the ICAL detector in INO

We present the first detailed simulation study of tau neutrino-induced charged current (CC) events from atmospheric neutrino interactions in the Iron Calorimeter (ICAL) detector at the proposed India-based Neutrino Observatory (INO) laboratory. Since the intrinsic atmospheric neutrino flux at few to 10s of GeV energy comprises only electron and muon neutrinos (and anti-neutrinos) with negligible tau neutrino component, any signature of atmospheric tau neutrinos is a signal for neutrino oscillations. We study the tau leptons produced through these CC interactions via their hadronic decay. These events appear as an excess over the neutral current (NC) background where hadrons are the only observable component. We find that the presence of tau neutrinos in the atmospheric neutrino flux can be demonstrated to nearly $4\sigma$ confidence with 10 years data; in addition, these events are sensitive to the neutrino oscillation parameters, $\sin^2\theta_{23}$ and $\vert \Delta m_{31}^2 \vert$ (or $\vert \Delta m_{32}^2 \vert$), in the 2--3 sector. Finally, we show that combining these events with the standard muon analysis which is the core goal of ICAL further improves the precision with which these parameters, especially the octant of $\theta_{23}$, can be measured.


Introduction
The proposed deep underground facility for carrying out research in basic sciences, the Indiabased Neutrino Observatory (INO), will house the magnetised Iron Calorimeter (ICAL) detector for carrying out atmospheric neutrino experiments. The ICAL detector will consist of 51 ktons of magnetised iron arranged in 151 layers, interspersed with Resistive Plate Chambers (RPCs) as the active detector. The experiment aims to precisely determine some of the important neutrino oscillation parameters using atmospheric electron and muon neutrinos (and anti-neutrinos) as sources [1].
The central goal of the proposed magnetised ICAL detector is the study of the neutrino mass ordering/hierarchy through the separate detection of muons and anti-muons (so-called standard muon events) produced in charged-current (CC) interactions of atmospheric muon neutrinos and anti-neutrinos, taking advantage of the magnetic field in ICAL. The detector will be optimised to detect muons (and associated hadrons) produced in the quasi-elastic, resonant, or deeply inelastic interactions of neutrinos with few to 10s of GeV energy. Such a detector has been shown, through detailed simulations, to be capable of making precision measurements of the 2-3 oscillation parameters sin 2 θ 23 and ∆m 2 32 , while being completely insensitive to the CP phase, δ CP [1]. It is also sensitive, due to Earth matter effects and the ability to identify the charge of muons, to the important open question of the neutrino mass ordering. The detector also has (reduced) sensitivity to these oscillation parameters through measurement of the electrons (and positrons) produced in CC interactions of electron neutrinos with the detector [2].
In this paper, we explore the other conventional physics that can be addressed by such a detector and study tau lepton production through atmospheric neutrinos in particular. Studies of sensitivity to electrons, as well as other exotic possibilities such as CPT/Lorentz invariance violation, neutrino decay, etc., have been studied by various members of the collaboration and the main results can be found in Ref. [1]. Here we present the sensitivity to the presence of tau neutrinos in the atmospheric neutrino flux that arise purely from oscillations of ν e and ν µ (and their anti-particles) in the atmospheric neutrino fluxes.
The Super Kamiokande collaboration has analysed their data for signatures of tau neutrinos and found that a no-tau hypothesis is rejected at the 3.8σ level [3]. They also determined the overall normalisation of the tau neutrino flux. Recently, the IceCube Collaboration has shown the presence of tau neutrinos arising from neutrino oscillations in atmospheric neutrinos to 3σ confidence and has also analysed the data to determine some of the neutrino oscillation parameters [4]. In fact, it was pointed out in Ref. [5] that the presence of tau neutrinos can be established at IceCube without knowing the oscillation parameters, and also without any assumption about the unitarity of the PMNS mixing matrix. Tau neutrino signatures are also being explored by the DUNE [6] and KM3NeT/ORCA experiments [7]. We show here that the presence of tau neutrinos in the atmospheric neutrino flux can be established to nearly 4σ confidence with 10 years data at ICAL. Moreover, we find that the tau sample (although contaminated with the neutral current (NC) sample) is indeed sensitive to the 2-3 parameters sin 2 θ 23 and ∆m 2 32 ; in addition, we show that a combined study of tau neutrino events with the standard muon events will significantly improve the precision with which these parameters can be measured at ICAL.
2 Interaction of tau neutrinos with the detector Atmospheric neutrinos are produced from the primary and secondary interactions of cosmic rays with Earth's atmosphere and comprise electron and muon neutrinos and anti-neutrinos. The majority of these interaction processes occur within a height of approximately 15 km from the Earth's surface. Since tau neutrinos arise from the production and decays of D ± S mesons in contrast to electron and muon neutrinos which arise from the decays of the lighter π ± pions, there is a negligible intrinsic tau neutrino flux in the atmosphere at energies less than 100 GeV. In fact, it has been shown that the intrinsic neutrino fluxes produced in the Earth's atmosphere via p p interactions are present in the ratio ν e : ν µ : ν τ :: 1 : 2 : 3 × 10 −5 [8]. Since muon neutrinos arise both from the decay of the pion component of cosmic rays as well as from the subsequent decays of the muons produced in these decays, the ratio of ν e : ν µ in atmospheric neutrinos is approximately 1 : 2 in the few GeV energy range. When these atmospheric neutrinos pass through the Earth to reach the detector, neutrino flavour oscillations occur, and tau neutrinos can arise via oscillations of both electron-and muon-type neutrinos. Consequently, the atmospheric neutrino species get re-distributed in the ratio 1 : 1 : 1, provided the ratio of the path length traversed to their energy, L(km)/E ν (GeV) 330 [9,10,11]. For neutrinos in the 10s of GeV range, this holds for the neutrinos entering the detector from below, or the so-called up-going neutrinos. Hence the atmospheric neutrino fluxes contain a significant fraction of ν τ in the upward direction (and practically none in the downward direction).
In principle, these tau neutrinos can be detected via their charged current (CC) interaction with the nucleons in the detector that produce charged τ leptons and hadrons via where X are hadrons (containing at least one nucleon to conserve baryon number). These τ leptons will decay promptly, primarily into hadrons. The branching fractions in the leptonic and hadronic channels are with similar fractions for τ + as well. It can be seen from the last expression in Eq. 2 that the dominant branching in the hadronic mode gives rise to additional hadrons, H ′ . The total energy in hadrons in such charged current interactions of tau neutrinos therefore arises from the contributions from X and H ′ , that is, from both primary production as well as subsequent decay. Hence it is expected that CC-tau events will present as events with high hadronic energy and no observed final state lepton. Due to the kinematics of the tau production and decay process, these hadrons are also peaked in the direction of the incident neutrino. In addition to CC interactions, Neutral Current (NC) interactions from all flavours of neutrinos also produce hadrons (X ′ ) via These cannot be distinguished from the CC-tau events, i.e., hadrons (X + H ′ ) due to CC interactions of ν τ and hadrons (X ′ ) due to NC interactions of all neutrino flavours, cannot be separated. Therefore, the NC events act as an inseparable background to the CC-tau events. Hence we will study the combined sample of all NC and CC-tau events in what follows. Note that the muon events arising from CC interactions of muon neutrinos (and antineutrinos) in ICAL give rise to events with a characteristic muon track associated with the hadron shower and hence can be distinguished from the NC or CC tau events which are present only as hadronic showers. Electron events from CC interactions of electron neutrinos (and anti-neutrinos) also present as showers since the electromagnetic shower of the electron and the hadronic shower from the associated hadrons cannot be separated in ICAL. Such showers can be distinguished from purely hadronic showers due to the larger number of hits per layer [2]; hence in the present analysis we assume that the NC and tau CC events can be (fully) separated from the CC electron and muon events although not from each other.
Relevant neutrino energies in our study are E ν > 3.5 GeV (which is the threshold for CC-tau production due to the large mass of the τ ). The processes involved are thus predominantly in the deep inelastic scattering (DIS) region. In this work, we have used the NUANCE neutrino generator to generate events for CC-tau and NC events in a simulated ICAL detector [12]. We use a GEANT4-based simulation toolkit to study the response of hadrons in ICAL with respect to energy and direction reconstruction [13]. Using this, we study the sensitivity of ν τ events (this henceforth refers to the combined CC-tau + NC sample) to the neutrino oscillation parameters. While these events are indeed sensitive to these parameters, a vast improvement is found on combining the tau sample with the standard CC muon sample. That will be discussed in the following sections.

Generation of CC-tau and NC events in ICAL
We use the NUANCE neutrino generator to generate the unoscillated events arising from both ν e and ν µ fluxes in the atmosphere. We then oscillate them in a 3-flavour model using the Pontecorvo-Maki-Nakagawa-Sakata mixing matrix (U). Detector-dependent resolutions and characteristics are then folded into the distributions. The events are analysed for their sensitivity to the neutrino oscillation parameters. We begin by decribing the generation of unoscillated events.

Generation of unoscillated events using the NUANCE neutrino generator
The NUANCE neutrino generator was used to generate CC-tau events with atmospheric muon and electron neutrinos [12]. The HONDA3d fluxes were used to generate events for 1000 years of exposure at the ICAL detector [14,15,16]. The distribution of tau events produced in the detector due to the interaction of these tau neutrinos with the material (mostly iron) of the detector is given by (4) A similar expression holds for tau anti-neutrinos as well. Here, T is the exposure time in seconds, N D the total number of targets in the detector, Φ e,µ the electron-and muon-type atmospheric neutrino fluxes, P iτ the relevant oscillation probabilities from flavours i → τ , and σ CC τ the CC cross section for ν τ interactions. Note that the oscillation probabilities are independent of the azimuthal angle and depend on (E ν , cos θ ν ) alone. In addition, the differential cross section, given by is used to obtain events in terms of the final lepton energy and angle. The NUANCE generator further fragments the final baryon into hadrons and lists the events in terms of the neutrino parameters (E ν , cos θ ν , φ ν ) as well as the individual parameters for each particle in the final state including the lepton: (E f , cos θ f , φ f ). When the final lepton is tau, it also completes the decay of the tau in either the hadronic or semi-leptonic modes and lists these as well. "Unoscillated" CC-tau events were generated assuming the ν µ and ν e atmospheric fluxes to be ν τ fluxes and using the CC-tau cross sections (for quasi-elastic, resonance and deep-inelastic processes) coded into NUANCE, that is, the events corresponding to the two terms in Eq. 4 were generated, excluding the oscillation probability. These unoscillated events are weighted with the appropriate oscillation probability during "data" generation and analysis. A similar procedure was used to generate the neutral current (NC) events using NU-ANCE. Since NC events are insensitive to oscillations, the NC events are simply generated from the ν e and ν µ atmospheric neutrino fluxes, using the NC cross sections instead of the CC cross sections shown in Eq. 4 and without including any oscillation probabilities. Both the CC-tau 1 and NC events give rise to events in the ICAL detector with hadron showers alone and no charged lepton track. We now discuss the implementation of neutrino oscillations.

Neutrino oscillation parameters
In the model with 3-flavour neutrino oscillation, the Pontecorvo-Maki-Nakagawa-Sakata mixing matrix (U) is commonly parametrised with three mixing angles θ 12 , θ 13 , θ 23 and one charge parity violating phase (δ CP ). If the neutrinos are Majorana particles there exist two additional Majorana CP phases; these are not visible in neutrino flavour oscillations. The neutrino flavour oscillation probabilities depend in general on the neutrino energy (E ν ), propagation distance between the source and the detector (L), parameters of the oscillation matrix (U) and the mass squared differences (∆m 2 ij = m 2 i − m 2 j , i = j) [17]. There are only two independent mass square differences to be considered for the case of three flavour neutrino oscillations.
The central values of the neutrino oscillation parameters and their 3σ ranges used in this analysis are given in Table 1. While the 1-2 parameters are kept fixed, the CP phase is poorly known [18]; moreover, the ICAL experiment is not sensitive to this phase. Hence we fix δ CP = 0 throughout this analysis. Since the octant of θ 23 is still unknown, its central value has been taken to be θ 23 = 45 • (sin 2 θ 23 = 0.5) in this analysis, unless otherwise specified.
The as-yet unknown neutrino mass hierarchy depends on the neutrino mass ordering: ∆m 2 31 > 0 for Normal Ordering , ∆m 2 31 < 0 for Inverted Ordering .
1 In the entire analysis, we always refer to the CC-tau events that decay hadronically. The 17% each of CC tau events that decay semi-leptonically producing muons or electrons will add (insignificantly) to the CC muon or electron events arising from the direct interaction of muon and electron neutrinos with the detector.
Since ∆m 2 21 > 0 and |∆m 2 21 | ≪ |∆m 2 31 |, this means that ∆m 2 32 and ∆m 2 31 have the same sign. For convenience, we define [19] Note that ∆m 2 flips sign without changing its magnitude when the hierarchy/ordering changes and hence is a convenient parameter compared to ∆m 2 31 or ∆m 2 32 , which change in both sign and magnitude depending on the mass ordering. We shall use ∆m 2 throughout in the analysis. Depending on the mass ordering, and using the value of ∆m 2 21 given in Table 1, the values of ∆m 2 31 and ∆m 2 32 can be found from ∆m 2 . We shall assume the normal ordering throughout this analysis, unless otherwise specified. More details and the current values and ranges of these parameters for both the normal and inverted mass ordering can be found in Refs. [18,20].

Neutrino oscillation probabilities
We use these values of neutrino oscillation parameters to generate the CC-tau events in ICAL. The threshold for this production is E th = 3.5 GeV. Since the atmospheric neutrino fluxes fall steeply with energy, hence the dominant contributions to CC tau events are at neutrino energies around 5-10 GeV. The tau oscillation probabilities were computed using the Preliminary Reference Earth Model (PREM) profile for the density distribution in the Earth using a 3-flavour model [21]. In the simulations, neutrinos of different flavours were propagated through the atmosphere and through Earth matter using a fast Runge-Kutta solver [19,22]. The oscillation probabilities P ατ , α = e, µ, are shown in Fig. 1 as a function of the zenith angle cos θ for a neutrino energy of E ν = 5, 10 GeV for both neutrinos and anti-neutrinos, using the central values of the oscillation parameters given in Table 1.
It can be seen that P µτ is much larger than P eτ so that the contribution from intrinsic muon atmospheric neutrinos and anti-neutrinos will dominate over the contribution from electron neutrinos.
Note also the discontinuity and features in both P eτ and P µτ at cos θ ≃ 0.86, which corresponds to the core-mantle boundary occurring in the neutrino oscillation probabilities due to the choice of the normal mass ordering. Such features will instead occur in P ατ if the inverted ordering was assumed.
While horizontal events are hard to measure due to the geometry of ICAL, it can be seen that the maximum contribution to tau events arises from the region cos θ ∼ 0.5 with  Figure 1: Oscillation probabilities (a) P ατ for neutrinos with E ν = 5 GeV, (b) P ατ for antineutrinos with E ν = 5 GeV, α = e, µ, as a function of the zenith angle cos θ (cos θ = 1 corresponds to UP neutrinos), with (c) and (d) corresponding to the same plots for energy E ν = 10 GeV respectively. a smaller contribution at cos θ ∼ 1.0 (due to the larger cross sections, neutrino events are about three times larger than anti-neutrino events).
The tau oscillation probabilities P ατ , α = e, µ, are shown in Fig. 2 as a function of energy for a zenith angle of cos θ = 0.5 (θ = 60 • ) for both neutrinos and anti-neutrinos, again using the central values of the oscillation parameters given in Table 1.  Figure 2: Oscillation probabilities (a) P ατ for neutrinos and (b) P ατ for anti-neutrinos, α = e, µ, as a function of the energy, E ν , for a zenith angle of cos θ = 0.5.
Note the presence of a broad maximum at E ν ∼ 10-15 GeV. We will see in the next section that this gives rise to CC tau events with hadron energies of 5-10 GeV, to which ICAL has reasonably good sensitivity. Before going on to the analysis of the sensitivity of CC-tau + NC events to the neutrino oscillation parameters, we briefly describe the simulation of the ICAL detector and its response to the hadrons of interest here.

Hadrons in ICAL
The proposed ICAL detector comprises 151 layers of 56 mm thick iron plates placed 40 mm apart, interleaved with the resistive plate chambers (RPCs) which are the active detector elements. This 51 kton detector will be magnetised upto about 1.4 T although this does not affect the hadron response [23]. When charged particles pass through the RPCs, they leave electrical signals called "hits" that are picked up which are localised to 3 cm × 3 cm × 0.2 cm in the x, y, z directions. The localisation in the vertical z direction is due to the small 2 mm gas gap in the RPCs. While the minimum ionising muons leave long tracks in the detector, the hadrons shower and hence can be calibrated in energy and direction only from the shape and number of hits in the shower [24]. A detailed study of the hadron energy response from an analysis of these hits can be found in Ref. [23] and has been used as-is in this analysis. A brief study of the dependence of the sensitivity to neutrino oscillation parameters on the hadron energy response is presented later in the paper since the analysis relies heavily on this observable. Here we will only highlight the angular response of the detector since it plays an important role in the analysis, due to the upward nature of tau events.  [26]; it can be seen that ǫ reco increases as cos θ increases.
As already mentioned, we are interested here in the events where the tau decays hadronically. There are two sets of hadrons in such events: the set of hadrons comprising the final state labelled X, which are the hadrons produced in the original interaction, and those labelled H ′ arising from tau decay.
Since the sets of hadrons cannot be distinguished and the tau decays rapidly (τ τ = 2.9 × 10 −13 s), all hadrons in each event are detected as a single shower in the detector. Due to the kinematics of both tau production and decay, the hadrons are peaked in the direction of the incident neutrino [25]. Earlier simulations studies of hadrons arising from NUANCE events showed that the hadron shower can be identified and their total energy calibrated by the hits in the RPCs, when the hadrons pass through them [23]. In addition, the direction of these hadrons can also be determined, although quite crudely compared to the direction of muons [26]. Due to the geometry of ICAL (with horizontal iron plates), the reconstruction ability in the vertical direction is expected to be better than for more horizontal events. While it was shown that the zenith angle in vertical directions can be reconstructed to within 10 • for 10 GeV hadrons, we are here interested only in the efficiency with which an upward/downward going hadron is reconstructed as an upward/downward going hadron, that is, reconstructed in the correct quadrant, and not in the actual zenith angle itself. The result of the simulation study is shown in Fig. 3 where the angular reconstruction efficiency of hadrons is shown as a function of both the true hadron energy and angle (cos θ). Here the reconstruction efficiency is defined as ǫ reco = Number of events reconstructed in the correct quadrant Total number of events reconstructed .
It is seen that the best direction reconstruction occurs at higher energies and for vertical angles (as expected), and is just about greater than 50% at more horizontal angles. We use this information in the analysis to separate the τ events into two direction bins, viz., up and down.
The neutrino energy and angle are used for calculating the oscillation probabilities, while the parameters used in the analysis are the summed reconstructed hadron energy in the event (denoted as where E τ is the final lepton energy), and the direction of the hadron shower (in two bins of UP/DOWN alone). From the direction reconstruction efficiency plot shown in Fig. 3, we see that a fraction (1 − ǫ reco ) of UP events will be reconstructed as DOWN events and vice versa. Details of the hadron energy resolution used in this analysis can be found in Ref. [23].

Generation of oscillated events at ICAL
The hadron energy and angle in the oscillated CC-tau events are smeared into the observed hadron energy and angle event-by-event, as per the detector response described above, and binned appropriately. Similarly, the NC events generated from NUANCE are also binned in the same bins. The dependence on the obtained sensitivity to the neutrino oscillation parameters due to errors in the hadron energy response is discussed in Section 4.8. Fig. 4 shows the energy distribution of a 10 year sample of CC-tau and NC events separately from neutrino and antineutrino sources (although they are not separable in the detector). The CC-tau events have been generated using the central values of the oscillation parameters given in Table 1. Fig. 4a shows the events for which the hadron shower direction is reconstructed as being in the UP direction (cos θ > 0) while Fig. 4b is for the DOWN events (cos θ < 0). It can be seen that there are a small number of CC-tau events reconstructed in the DOWN bin due to the error in direction reconstruction, as seen already from Fig. 3. However, the increase over the NC background is clearly visible for events in the UP bin, where the error bars shown are statistical. We are now ready to analyse these events for their sensitivity to neutrino oscillation parameters.

Numerical analysis and physics reach
We now go on to a numerical study to examine the sensitivity of NC + CC tau events to various neutrino oscillation parameters of interest.

Generation of data samples
Statistical fluctuations are significant in the analysis of neutrino events for their sensitivity to the oscillation parameters. This is in fact what distinguishes the sensitivity for different exposure times. The statistical limitations of data contribute in two differing ways. One is the error or precision with which a neutrino oscillation parameter (say θ 23 ) can be measured. It is determined by the amount of data available for the analysis. The other limitation is the fluctuations present in the data sample which may yield a best fit value for the parameter which may not coincide with its true value. In simulations analysis, therefore, different sets of "data" samples, all randomly generated for the same n number of years, will yield different best values. Hence using a single arbitrary "data" sample in such simulations analysis raises the risk of over-or under-estimation of sensitivity for a given input value. This can be avoided by repeating the analysis N times with different "data" sets of generated data samples for the given input value. The average sensitivity in these repeated analyses (for N 60), will approach a median sensitivity for the given input value of oscillation parameter; in fact, the various results obtained with each sample will cluster around this median sensitivity.
Alternatively, using a procedure which enormously saves computational time, a large data sample (1000 years) is generated and scaled to the required n number of years (typically 10 years) during the analysis. This procedure will yield the median sensitivity for a given input value of parameters in the χ 2 analysis for which the sensitivity is being measured. It was shown in several analyses [27,28] that such a procedure of generating large data samples and scaling them to the required number of years correctly determines the precision with which parameters can be determined. We will use this technique for our analysis in what follows and examine this in more detail with an example in Section 4.7.

Best fit approach
A large data sample (1000 years) of CC-tau events was generated as per the process described above. A similar 1000 year sample of NC events was generated as well. Notice that the source of atmospheric neutrino fluxes is the same in both cases. The simulated "data" is generated by applying neutrino flavour oscillations using a set of input values of oscillation parameters (typically the central values of the oscillation parameters given in Table 1), and scaling this sample to 10 years, as required. In order to test the sensitivity of this sample to neutrino oscillations, the original unoscillated sample is oscillated using a different value of one or more oscillation parameters, scaled to the same number of years, and labelled "theory". The χ 2 for this set of "data" and "theory" is defined as where, is the total number of CC-tau and NC events arising from both antineutrino (+) and neutrino (−) fluxes, generated with the set of input values of oscillation parameters, in the i th energy and j th cos θ bin, defined to be the "Data" in the simulations study, T ij is the corresponding number of predicted theory events in the same bins, generated using a different set of oscillation parameters, where we have included the systematic uncertainties via the pulls technique; see Refs. [29,30], so that the number of theory events includes the systematic uncertainty from five sources which are described in detail in Section 4.3: Here T ij,0 ± are the number of antineutrino/neutrino theory events, without systematic errors in the corresponding bins, and ξ 2 k ≡ ξ 2 k,− +ξ 2 k,+ includes the penalty from each pull parameter for neutrino and anti-neutrino events respectively.
A measure of the sensitivity of the data to the input value of any parameter is given by ∆χ 2 , defined as where χ 2 (input) corresponds to the minimum χ 2 obtained when the theory events are calculated with the same value of the parameter as its input value, and χ 2 (par) is the minimum value obtained when a different theory value of the parameter is used. Note that when the scaling procedure is used to generate the "data", the value of χ 2 (input) = 0. The definition can be appropriately extended to the case when more than one parameter is varied from its input value when calculating the "theory" events. Finally, a prior on the well-known parameter sin 2 2θ 13 is included via Note that sin 2 2θ 13 is very well constrained; also, while the sensitivity of ICAL to the neutrino mass ordering depends on this parameter, it is not very sensitive to the value of this parameter itself, and hence any changes in the central values used for this parameter will not affect our results. While including the systematics, the minimisation with respect to the systematic nuisance parameters is done by analytically solving for the set of ξ i for a given set of oscillation parameters. The minimum χ 2 is then found by varying the oscillation parameters for the "theory" set of events and marginalising over irrelevant parameters as required.

Systematic uncertainties in the analysis
The inherent systematic uncertainties associated with the prediction of the events and their rates affect the sensitivity of these events to the oscillation parameters. There are five different types of systematic unertainties which are considered in our analysis, viz., the relevant sum in Eqs. 9 and 10 runs over k = 1, · · · , N k = 5 [1]. These values are standardly used by the INO Collaboration in all its analyses [1].
1. We calculate the energy dependent flux tilt error by considering a deviation of δ = 5% from the standard behaviour, viz., E −2.7 ν . Hence the systematic error π tilt from this uncertainty can be calculated for each energy bin through We choose the standard value of E 0 = 2 GeV [15].
2. We consider the flux angular uncertainty to be π zenith = 5% cos θ, in the given zenith angle bin.
3. We take the overall flux normalisation uncertainty to be π norm = 20%.
4. The systematic error due to uncertainty in computing the cross section is taken to be π σ = 10%.
5. Finally, an overall uncertainty of π D = 5% is included to take care of any uncertainty in characterising the detector response.
We have summarised the systematic uncertainties in Table 2. Note that we have taken the uncertainties to be the same for the neutrino and antineutrino samples, and these events are combined to compute the χ 2 since they are not distinguishable in the detector. Note also that the last three overall bin-independent errors can be replaced numerically by a single one: π = π 2 norm + π 2 σ + π 2 D ; they are separately retained to allow for later refinement in the inclusion of systematic errors. We discuss in Section 4.8 the effect of an energydependent uncertainty while implementing the hadron energy response of the detector, which is a preliminary study on the bin dependence of the pull corresponding to the detector response, π D .

The binning scheme
We have optimised the number of hadron energy bins by computing the χ 2 values for various sets of data and theory. The final set of bins used in the analysis is listed in Table 3.

Sensitivity to the presence of tau events
We begin by asking whether ICAL will be sensitive to the presence of the CC-tau events in the sample. That is, we consider the situation where the "theory" does not account for the presence of tau neutrinos in the sample. Hence the theory events arise purely from the NC events and are independent of neutrino oscillations. We find that, without including sytematic uncertainties, the significance to the presence of tau is at a confidence level of about ∆χ 2 ∼ 6σ for 10 years exposure; when the systematic uncertainties are included as described above, this decreases to 3.6σ; see Table 4. The major impact is due to the tilt and the zenith angle uncertainties, while the impact of the remaining uncertainties is small or negligible. For instance, if the zenith angle and tilt uncertainties reduce to 70% or 50% of the values given in Table 2 by the time ICAL is operational, the sensitivity will increase to 3.9σ and 4.1σ respectively. Hence it is possible that sensitivity to the presence of tau neutrinos can be achieved to nearly 4σ level with ICAL.  Table 4: Sensitivity to the presence of tau neutrino events in various scenarios including the number of angular bins, and whether systematic uncertainties were included or not. The effects of a reduction in tilt and zenith angle uncertainties to 70% and 50% of their values given in Table 2 are also shown.
Also note the effect of separating the events in UP/DOWN angle bins. If no angular information is taken into account, the significance falls to just 2.8σ. This highlights the importance of angular information in this analysis, due to the fact that almost all tau neutrinos are produced in the upward direction. An improvement in the angular reconstruction of hadron showers will therefore vastly improve this result, which is comparable to that obtained by the SuperKamiokande Collaboration also with atmospheric neutrino data [31]. This is discussed further in the next two sections.

Sensitivity of tau neutrino events to the oscillation parameters
We now proceed to study the sensitivity of the tau events (with NC background) to the neutrino oscillation parameters, in particular, to the 2-3 parameters sin 2 θ 23 and ∆m 2 . We begin with the mixing angle, sin 2 θ 23 . Now the data and theory include both the NC and CC-tau events, with the "theory" events being generated with a different value of sin 2 θ 23 than the input value of sin 2 θ in 23 used to generate the data. The resulting ∆χ 2 , defined in Eq. 11, is then a measure of the sensitivity of the data to the input value of the parameter. Fig. 5a shows the resulting ∆χ 2 as a function of the value of sin 2 θ 23 used to generate the theory events, for an input value of sin 2 θ 23 = 0.5 (θ in 23 = 45 • ). Here sin 2 θ 23 is kept fixed while ∆m 2 and sin 2 2θ 13 are marginalised over their 3σ ranges as listed in Table 1; the systematic uncertainties mentioned above are also included in the analysis. It can be seen that including the zenith angle dependence (using the information on UP/DOWN bins shown in Fig. 3) improves the sensitivity even though the hadron direction is poorly determined. The same trend is seen for ∆m 2 , as seen in Fig. 5b. Henceforth we always consider the data binned in both hadron energy and angle. It may be argued that including even more zenith angle bins will further improve the sensitivity. Indeed, when the up-going events are binned in 2 zenith angle bins (0 • < θ < 45 • and 45 • < θ < 90 • ) while retaining the down-going events in a single bin, there is a small improvement in the sensitivity. However, due to the limited reconstruction capability of the zenith angle of hadrons, there will be large correlations between these bins; such an analysis requires a deeper study of the angular reconstruction of hadrons in ICAL than is available at present. In addition, there is also a modest improvement when the hadron energy resolution is improved, for instance, on using the energy resolution that would be obtained if ICAL used 2 cm iron plates rather than the design value of 5.6 cm.
We find that the sensitivity is dominated by the systematic uncertainties as can be seen in Fig. 6a where the results of an analysis with no consideration of systematic uncertainties and either keeping all parameters fixed at their central values (other than sin 2 θ 23 , of course) or marginalised over their 3σ ranges listed in Table 1 are shown. The analysis is insensitive to the value of the CP phase and this as well as the 1-2 oscillation parameters have been kept fixed to their central values listed in Table 1. In addition, the neutrino ordering is assumed to be normal, unless stated otherwise. While there is practically no effect on including marginalisation, the curve labelled "Marg., with pulls" indicates clearly that the maximum loss of sensitivity occurs on including systematic errors.
In Fig. 6b, we see the sensitivity to the oscillation parameter sin 2 θ 23 as a function of the exposure time. The results correspond to inclusion of systematic errors and marginalisation over the remaining parameters. While the sensitivity to the oscillation parameter precision measurement improves as the number of years increases, this is not exactly linear in the number of years due to systematic effects.
In Fig. 7, we have shown the sensitivity study to the oscillation parameter |∆m 2 |. It is seen that the tau events are able to discriminate better against |∆m 2 | values that are lower than the input value since the dependence on this parameter is via the dominant P µτ oscillation probability and arises as sin 2 [∆m 2 L/E ν ]. Finally, we remark that the analysis indicates negligible sensitivity to the neutrino mass ordering and is therefore insensitive to the sign of ∆m 2 (that is, to the sign of ∆m 2 31 or ∆m 2 32 ). However, we believe that this is the first study to explore the possibility of extracting neutrino oscillation information from

A discussion on scaling the data sample
Here, we show with an example how the procedure of generating the "data" sample using the 1000 year generated events scaled to the required number of n years yields the correct precision with which the parameters are determined. In this example, we calculate the ∆χ 2 values for a theory value of sin 2 θ par 23 = 0.25 (θ 23 = 30 • ) using 10 different "data" sets of 10 years data generated without scaling. Here the other oscillation parameters are kept fixed and systematic uncertainties have been ignored, for clarity.
As we see from Fig. 8a, the calculated ∆χ 2 values for these samples are clustered around the value of ∆χ 2 = 1.833 which is the value obtained using the alternate method when the 1000 year set is scaled to 10 years and used as "data". It clearly indicates the risk of over-or under-estimation of sensitivity to sin 2 θ 23 if we only use a single sample which was randomly generated for 10 years.
We also further see in Fig. 8b the consistency between the procedure of scaling the events to the required number of years and the procedure when actual data for n years are used with no scaling. It shows the ∆χ 2 obtained when all parameters are kept fixed (for convenience) and the "theory" is generated with sin 2 θ input 23 = 0.25, as a function of the number of years of exposure in ICAL. The smooth red line corresponds to the sensitivity obtained when 1000 years of events are taken and scaled to the required number of years (1, 2, · · · , 20) to generate the "data". The green histogram, in contrast, corresponds to the case when "data" are generated for the exact number of years required and then compared to the (scaled) "theory". It can be seen that the trend of the two lines is the same, and the ∆χ 2 value in the second case fluctuates about the median red line obtained with scaling. This validates our use of the scaling procedure in our analysis.

Effect of errors in hadron energy reconstruction
It can be seen that the primary sensitivity to tau events over the NC background occurs because the former correspond to higher hadron energies. Hence the results obtained in the previous sections are dependent on the correct hadron energy reconstruction. We examine briefly here the effect of errors in the hadron energy reconstruction.
One of the important issues in the analysis is the hadron energy reconstruction. ICAL has rather poor sensitivity to the hadron energy, compared to its ability to reconstruct muons. In our analysis, the hadron energy of the events generated by the NUANCE generator is smeared event-by-event and binned appropriately in the observed hadron energy bins. In order to understand the dependence of the sensitivity on the hadron energy, we have therefore simulated the following: the "data" is generated according to the "true" hadron energy reconstruction determined by the simulations group of the collaboration [23]; the "theory" however, is generated using a different width, while retaining the correct central value.
Three sets of widths were used, with σ/E of the reconstructed hadron energy distribution being 5%, 10%, and 50% worse than the true value 2 . As a consequence, some of the events which would have been binned in a given energy bin may now be binned in any of the adjacent bins. Since there are more events at low enery, especially below 6 GeV, as can be seen from Fig. 4, the application of a worse reconstruction width for the hadron energy causes more of the low energy events to smear into even lower energies (and hence are lost to the analysis if the reconstructed energy is less than 1 GeV), or into the higher energy, E H > 6 GeV, bins. This results in the higher energy bins having a larger number of events than the "data" set, even when the oscillation parameters are not changed. The dominant events at low energies are NC events, which are independent of oscillations; the mis-match caused by this error in hadron energy reconstruction can only be compensated by changing the CC-τ events, using a different set of oscillation parameters; this change has to be quite large due to the smaller number of CC-τ events compared to the NC sample.
In Fig. 9, we compare the sensitivity to sin 2 θ 23 when the "true" hadron energy response is used for both the "data" and "theory" with that when the latter has a 10% larger width than the former. Here the sensitivity is defined as the difference in χ 2 : where χ 2 0 is the minimum value of χ 2 when both "data" and "theory" use the same hadron energy reconstruction. While the sensitivity worsens a little, and the minimum ∆χ 2 is no longer at the true value (see below for more details on this), there is still consistency of the results at 1σ, i.e., ∆χ 2 = 1.
True Eh 1.1Eh Figure 9: The dependence on the hadron energy response. The figure shows the change in the sensitivity to sin 2 θ 23 when the "theory" events are smeared by a hadron reconstruction code that has 10% larger width than "data". See text for details. Fig. 10 shows the two dimensional contour plot of allowed values in sin 2 θ 23 and ∆m 2 at the 70% CL along with the best fit value (shown as a plus). The best fit values when the "theory" events are reconstructed according to a hadron energy response with a width 5%, 10% and 50% larger than that used for the "data", are also shown (as a cross, a filled square and a filled circle respectively).
We see that the effect of mismatch in true and fitted hadron energy response is higher in sin 2 θ 23 compared to ∆m 2 as the deviation in the best fit value is rather small for the latter. Both 5% and 10% worse widths give acceptable central values, with an acceptable the worse results. Since the Vavilov has a longer tail than the corresponding Gaussian distribution, we have therefore ignored higher energy tails that would actually improve the result we obtain. True Eh Best, 1.05Eh Best, 1.10Eh Best, 1.50Eh Figure 10: 70CL contour plot of allowed parameter range in the sin 2 θ 23 -∆m 2 plane with the best fit value marked with a 'plus'. The best fit values when the "theory" events are reconstructed with hadron energy resolution having 5%, 10% and 50% larger widths are also shown as cross, square and circled points. See text for details. minimum ∆χ 2 (< 0.4). When the width used for "theory" is as much as 50% worse than for the "data", the best fit point lies well outside the contour. In addition, the minimum χ 2 in this case is ∆χ 2 min = 12.4, which is nearly 3.5σ worse than when the correct hadron energy reconstruction is used for both "data" and "theory". This is an important consideration in the tau neutrino analysis. The mini-ICAL prototype has been used [32] to validate the energy and momentum response of the detector to muons. However, there is no data as yet to validate the hadron simulations results that have been obtained so far, although the simulations themselves have been validated against data from the MONOLITH experiment [33]. A more detailed analysis using bin-to-bin correlations will allow for a detailed study of the effect of systematic errors in the hadron energy response and will improve the quality of the tau neutrino analysis. At this point, the short summary is that the analysis will tolerate about 10% deviations from the expected hadron energy response and such a detector will continue to be sensitive to these events.

Combined study of tau and "standard" muon events
So far we have examined the sensitivity of pure hadron events (including NC and CC-tau events) in ICAL to the neutrino oscillation parameters. Notice that these arise from the same atmospheric neutrino fluxes as the standard muon events in ICAL that are the main goal of this detector [1]. While the unoscillated muon neutrino fluxes give rise to the dominant CCmuon events in the detector via P µµ , the electron neutrino fluxes also contribute to this signal via P eµ . Hence, uncertainties such as overall flux normalisation, zenith angle dependence, and energy tilt error are the same for both sets of analyses. Since tau is massive, the CC-tau cross section is highly suppressed at smaller energies, E ν 5 GeV. However, due to the larger threshold, tau production events dominantly arise from deep inelastic scattering and it is reasonable to assume that the cross section uncertainties (for CC-tau, NC, and CC-mu) are roughly the same. It then becomes obvious that the systematic uncertainties for the dominant CC-mu processes are the same as for the CC-tau + NC processes currently being studied. Since the systematic uncertainties are the dominant factor limiting the sensitivity in the tau analysis, there is expected to be a significant improvement on combining the analyses from the two data sets. That is, given that the two data sets have common systematic uncertainties, the sensitivity (∆χ 2 ) of the combined analysis is expected to be better than just the sum of the two individual values.
Previous simulations studies of the INO collaboration have demonstrated the capability of ICAL with respect to precision measurement of the 2-3 neutrino oscillation parameters: sin 2 θ 23 and its (as-yet unknown) octant, ∆m 2 (including its sign), while being insensitive to the CP phase δ CP [1,22]. This was done using the CC muon events generated when muon neutrinos interact with the ICAL detector, producing charged muons whose momentum, direction and sign of charge can be accurately reconstructed from the long tracks they leave in the (magnetised) detector using a Kalman filter-based alogrithm [1]. Since the tau events indicate (albeit admittedly limited) sensitivity to these parameters, we now go on to a combined analysis of these data sets.
The sensitivity to a given neutrino oscillation parameter is again defined through the χ 2 : where the individual contributions are given by and χ 2 (prior) is given by the second term on the right hand side of Eq. 12, with ξ k summing over all nuisance parameters as described earlier. Note that due to the ability to reconstruct the sign of the muon charge, the contribution from the individual neutrino and anti-neutrino events are considered for the standard muon analysis while the summed contribution is analysed for the tau events. Here the "tau" contribution is understood to include both CCtau and NC events. The terms corresponding to the "data" and "theory" events for the standard muons are analogous to those for the tau events given in Eq. 10, and binned in the three variables, E µ , cos θ µ and E Hµ , corresponding to the muon momentum magnitude, muon zenith angle, and total hadron energy in the muon CC event.
Note on coding the systematics : The µ + and µ − events are separately analysed and their individual contributions to the χ 2 determined. The theory contributions to each of these involves either the set ξ − or ξ + of systematic uncertainties; however, due to a small charge identification inefficiency (less than 2% for muon momenta beyond about 1 GeV/c), there is a small fraction of "wrong-sign" events in each events sample. In the case of the CC-tau+NC analysis, of course, the neutrino-and anti-neutrino-induced events are added together. Hence, while solving for the set of ξ min k , a certain simplification is applied: since the anti-neutrino events are about three times smaller than the neutrino events (due to their relatively smaller cross sections), terms of the order of (N + /N − ) 2 are dropped from the expressions. While this error is very small for the case of standard muons, it is about 2-3% for the case of the tau events. However, with this approximation, it was found possible to implement a fast analytical invertor for the corresponding 10 × 10 pulls matrix, which speeded up the analysis extensively.
Effect on sensitivity to sin 2 θ 23 : Fig. 11 shows the effect of combining the tau events with the standard muon events on the 10 year sensitivity to the oscillation parameter sin 2 θ 23 , for the input value sin 2 θ in 23 = 0.5 (θ in 23 = 45 • ). While the sensitivity of the tau events alone is marginal in the range of sin 2 θ 23 shown in the figure, it significantly improves the precision reach for this parameter. This is defined as where ∆V p n is the allowed range of the values of the parameter p at nσ, when the remaining parameters are marginalised over their 3σ ranges, and V p 0 is its central value. At 2σ, we see that the precision P 2σ (sin 2 θ 23 ) reduces from 11% to 9.5% on the inclusion of the tau events, and improves even more significantly from 15% to 12% at 3σ. This is because the source of atmospheric neutrinos is the same in both cases, and this helps reduce the effects of including systematic uncertainties, which are common to both analyses.  Figure 11: 10 year sensitivity of ICAL to the oscillation parameter sin 2 θ 23 with NC + CC tau events alone (this analysis), from standard muons alone (old analysis [1,22]), and a combined analysis of the two data sets for sin 2 θ in 23 = 0.5. It can be seen that while the sensitivity of tau events to the parameter is very small in the range of values shown, it causes a dramatic improvement in the precision measurement of this parameter when combined with the standard muon events.
Sensitivity to the octant of θ 23 : This question is important for model builders. The dominant contribution to the survival probabilities such as P µµ is from the octant-insensitive term sin 2 2θ 23 and the dependence on the octant-sensitive sin 2 θ 23 is proportional to sin 2 2θ 13 and hence is small. There is a larger octant dependence in the oscillation probabilities (P ij , i = j), but all such dependences are also modulated by sin 2 2θ 13 , thus making it challenging to measure.   Table 1.
As seen from Fig. 12a, there is a significant improvement in the case of the input value of this parameter being in the lower octant, viz., θ in 23 = 40 • , with the combined analysis being able to discriminate against the maximal mixing value of θ in 23 = 45 • (sin 2 θ in 23 = 0.5) at 2σ. Hence, while the standard muon analysis could not distinguish the octant (or even deviation from maximality) for the input value of θ in 23 = 40 • , the combined analysis can do both. As expected, the improvement is more modest for the case when the input value is in the upper octant, with θ in 23 = 50 • (see Fig. 12b), due to the nature of the dependence of P µτ on sin 2 θ 23 ; however, as seen earlier, the inclusion of the tau events improves the precision to which sin 2 θ 23 can be determined in every case.
Effect on sensitivity to ∆m 2 : Fig. 13 shows the effect of combining the tau events with the standard muon events on the 10 year sensitivity to the oscillation parameter ∆m 2 , for the input value ∆m 2 in = 2.47 × 10 −3 eV 2 . It is seen that the sensitivity on including the tau events is marginal. Note that the results from standard muons alone are slightly different from those shown in Ref. [22] due to the slightly different central values of parameters used in the analysis.
Sensitivity to the mass ordering or the sign of ∆m 2 : The combined standard muon and tau+NC events have a marginally better sensitivity to the sign of ∆m 2 . For instance, the value of ∆χ 2 from the combined analysis is 0.5 better than that obtained with standard   Figure 13: Sensitivity to oscillation parameter ∆m 2 from standard muons alone, tau events alone, and from combining the two sets in a simultaneous analysis.
muons alone for the central values listed in Table 1, assuming the normal ordering and 10 years' running of ICAL [22]. There is hardly any improvement if the true hierarchy is assumed to be inverted.

Discussion and Conclusions
Tau neutrinos do not naturally occur in the atmospheric neutrino spectrum, which comprises electron and muon neutrinos (and anti-neutrinos) primarily arising from pion and subsequent muon decays of the primary cosmic ray spectrum. Due to neutrino oscillations, those atmospheric neutrinos that traverse significant path lengths can oscillate into one another as well as into tau neutrinos. Hence it is expected that a significant fraction of upward-going atmospheric neutrinos are of the tau flavour. Such tau neutrinos can provide an independent test of neutrino oscillations through their direct detection via charged current (CC) interactions of these neutrinos with the material of the detector. In this paper we have made a detailed simulations study of such a process at the proposed magnetised ICAL detector at the Indiabased Neutrino Observatory, INO. While ICAL is being optimised to detect charged muons from the CC interactions of atmospheric muon neutrinos (and anti-neutrinos), the so-called standard muon sample, it is also sensitive to hadrons that are produced along with muons in the CC interaction.
In particular, we have analysed the events where the charged taus produced in the CC interaction decay hadronically so that the event is rich in hadrons (produced at the interaction vertex as well as during tau decay). Due to the characteristics of the ICAL detector, such events are indistinguishable from neutral current (NC) events where the final state neutrino escapes and only a hadron residue is observable. Hence both CC-tau events and NC events are analysed together in this work. Since the tau production threshold is high, E ν > 3.5 GeV, due to the large mass of the tau lepton, the tau events are few in number but make a significant addition over the high energy NC events which are also limited due to the steeply falling neutrino spectrum (∝ E −2.7 ν ).
Another signature of tau events is the fact that they are exclusively produced by tau neutrinos moving in the upward direction. Since the neutrino energies involved are sufficiently large, the final state tau as well as its decay products are also peaked in the same direction as the initial neutrino. Hence good angular discrimination should enable extraction of these events over the uniformly distributed NC events. Due to the limitations of detector reconstruction, the direction resolution of hadrons is rather poor in ICAL. However, the addition of just two angle bins (corresponding to UP going and DOWN going events) leads to a significant sensitivity of these events to the presence of tau neutrinos. We have found that ignoring the tau component in the theory fits to the simulated data (which contains both NC and CC-tau events) leads to a mismatch with ∆χ 2 = 38 which drops to 15 when various systematic errors are included in the analysis, thus indicating that the pure hadron events sample can unambiguously signal the presence of tau neutrinos in the atmospheric neutrino flux.
For the first time, we have also analysed the CC-tau + NC sample for its sensitivity to the neutrino oscillation parameters themselves. A modest sensitivity to both the 2-3 parameters sin 2 θ 23 and |∆m 2 | (although the sample had no sensitivity to the sign of ∆m 2 ) was found, mainly limited by the systematic uncertainties. Although small, the result was encouraging since this data sample originates from the same atmospheric neutrino fluxes as the standard muon sample that has been extensively studied by the INO collaboration [1]. Hence systematic unertainties pertaining to the fluxes, such as overall normalisation, zenith angle, and energy dependent tilt uncertainties are the same for both samples. In addition, at higher energies dominated by deep inelastic scattering processes, the cross section uncertainties can also be considered to be the same for the two sets. Hence it is possible to perform a combined analysis of the CC-tau+NC and the standard muon sample thereby increasing the statistics without worsening the systematic uncertainties. Such a combined analysis gave a sensitivity (∆χ 2 ) that was significantly better than just the sum of the individual contributions. In particular, while there was not much improvement in the precision of |∆m 2 |, there was a significant improvement in the precision of sin 2 θ 23 . In addition, it was found that the sensitivity to the octant of θ 23 significantly improved. For instance, it is possible to determine both the octant as well as establish deviation from maximality at 2σ from a 10 year combined sample when the input value of θ 23 was in the first octant, θ in 23 = 40 • . A somewhat more modest improvement was seen when the input value of θ 23 was assumed to be in the second octant, as expected. A small improvement in the determination of the neutrino mass ordering was found when the true ordering was assumed to be normal; no such improvement was seen when the true ordering was assumed to be inverted.
Neutrino experiments are low counting experiments. Hence it is important to analyse every possible channel to yield more information on the neutrino oscillation parameters. Combining tau events with standard muon events opens up a way of improving the precision and possible measurement of parameters such as the as-yet unknown octant of the 2-3 oscillation parameter, θ 23 . Many current and up-coming experiments are focussing on this relatively unknown sector. In addition, information from the tau sector can also probe the 3-flavour structure of the PMNS matrix and unitarity violation [34]. This should give even more impetus to the study of such tau events.