From oscillation dip to oscillation valley in atmospheric neutrino experiments

Atmospheric neutrino experiments can show the"oscillation dip"feature in data, due to their sensitivity over a large $L/E$ range. In experiments that can distinguish between neutrinos and antineutrinos, like INO, oscillation dips can be observed in both these channels separately. We present the dip-identification algorithm employing a data-driven approach -- one that uses the asymmetry in the up and down events, binned in the reconstructed $L/E$ of muons -- to demonstrate the dip, which would confirm the oscillation hypothesis. We further propose, for the first time, the identification of an"oscillation valley"in the reconstructed ($E_\mu$,$\,\cos\theta_\mu$) plane, feasible for detectors like ICAL having excellent muon energy and direction resolutions. We illustrate how this two-dimensional valley would offer a clear visual representation and test of the $L/E$ dependence, the alignment of the valley quantifying the atmospheric mass-squared difference. Taking into account the statistical fluctuations in the 10-year simulated data at ICAL, we estimate the precision to which atmospheric neutrino oscillation parameters would be determined using our procedure. Our approach can be adopted by other atmospheric neutrino experiments to establish neutrino oscillations in terms of ($L$,$\,E$) rather than $L/E$.


Introduction and Motivation
We have witnessed remarkable developments in neutrino physics over the last two decades, driven by the astonishing discovery of "neutrino mass-induced flavor oscillation" suggesting that neutrinos have non-degenerate mass and they mix with each other [1,2]. Due to these two unique features, neutrinos change their flavor as they move in space and time implying that leptonic flavors are not symmetries of Nature [3]. Such a path-breaking discovery of fundamental significance, recently recognized with the Nobel Prize [4], has created enormous interest in the global neutrino community to measure the fundamental oscillation parameters with much better precision [5]. Such measurements of neutrino oscillations parameters are essential to provide a stringent test of the standard three-flavor neutrino oscillation framework, which has six fundamental parameters: (i) the solar masssquared difference, ∆m 2 21 (≡ m 2 2 − m 2 1 ≈ 7.4 × 10 −5 eV 2 ), (ii) the solar mixing angle, θ 12 ≈ 34 • , (iii) the atmospheric mass-squared difference, |∆m 2 32 | (≡ |m 2 3 − m 2 2 | ≈ 2.5 × 10 −3 eV 2 ), (iv) the atmospheric mixing angle, θ 23 ≈ 48 • , (v) the reactor mixing angle, θ 13 ≈ 8.6 • , and (vi) the Dirac CP phase, δ CP , which at present lies in the third quadrant with large uncertainty.
It is quite remarkable to see that starting from almost no knowledge of the neutrino masses or lepton mixing parameters twenty years ago, we have been able to construct a robust, simple, three-flavor oscillation framework which successfully explains most of the data [6][7][8][9]. Atmospheric neutrino experiments have contributed significantly to achieve this milestone [10,11] by providing an avenue to study neutrino oscillations over a wide range of energies (E ν in the range of ∼100 MeV to a few hundreds of GeV) and baselines (L ν in the range of a few km to more than 12,000 km) in the presence of Earth's matter with a density varying in the range of 0 to 10 g/cm 3 .
An important breakthrough in the saga of atmospheric neutrinos came in 1998 when the pioneering Super-Kamiokande (Super-K) experiment reported convincing evidence for neutrino oscillations in atmospheric neutrinos by observing the zenith angle (this angle is zero for vertically downward-going events) dependence of µ-like and e-like events [12]. The data accumulated by the Super-K experiment, based on an exposure of 33 kt·yr, showed a clear deficit of upward-going events in the zenith angle distributions of µ-like events with a statistical significance of more than 6σ, while the zenith angle distribution of e-like events did not show any significant updown asymmetry. These crucial observations by the Super-K experiment were successfully interpreted in a two-flavor scenario assuming oscillation between ν µ and ν τ , leading to the disappearance of ν µ . In this scenario, the survival probability of ν µ can be expressed in the following simple fashion: P (ν µ → ν µ ) = 1 − sin 2 2θ 23 · sin 2 1.27 · |∆m 2 32 | eV 2 · L ν (km) E ν (GeV) . (1.1) Due to the hierarchies in neutrino mass pattern (∆m 2 21 |∆m 2 32 |) and in mixing angles (θ 13 θ 12 , θ 23 ), the above simple two-flavor ν µ survival probability expression was sufficient to explain the following broad features of the Super-K atmospheric data, providing a solution of the long-standing atmospheric neutrino anomaly in terms of "neutrino mass-induced flavor oscillations".
• Around 50% of upward-going muon neutrinos disappear after traveling long distances.
For very long baselines, the above ν µ survival probability expression (see Eq. 1.1) approaches 1 − (1/2) sin 2 2θ 23 and the observed survival probability becomes close to 0.5, suggesting that the mixing angles is close to the maximal value of 45 • .
• The disappearance of ν µ events starts for the zenith angle close to the horizon indicating that the oscillation length should be around ∼ 400 km for neutrinos having energies close to 1 GeV. This suggests that the atmospheric mass-squared difference, |∆m 2 32 | should be around 10 −3 eV 2 .
• There is no visible excess or deficit of electron neutrinos. It suggests that the oscillations of muon neutrinos mainly occur due to ν µ → ν τ transitions.
• The Super-K experiment was also the first experiment to observe the sinusoidal L/E dependence of the ν µ survival probability [13]. A special analysis was performed by the Super-K collaboration, where they selected events with good resolution in L/E, largely excluding low-energy and near-horizon events. Using this high-resolution event sample, they took the ratio of the observed and expected event rates and the oscillation dip started to appear on the canvas around L/E ∼ 100 km/GeV, with the deepest location of the dip around L/E ∼ 500 km/GeV (see Fig. 4 in Ref. [13]). This study was quite useful to disfavor the other alternative models which showed different L/E behaviors.
The proposed atmospheric neutrino detector Iron Calorimeter (ICAL) at the Indiabased Neutrino Observatory project aims to detect atmospheric neutrinos and antineutrinos separately over a wide range of energies and baselines. The performance of ICAL is optimized for the reconstructed muon energy (E rec µ ) range of 1 GeV to 25 GeV, and reconstructed baseline (L rec µ ) from 15 km to 12000 km except around the horizon. Thus, the range of L rec µ /E rec µ of muon to which ICAL is sensitive, is from 1 km/GeV to around 10 4 km/GeV. Therefore, we expect to see an oscillation dip (or two) as in Super-K. In this paper, we investigate the capability of the ICAL detector to reconstruct the oscillation dip in the reconstructed L/E distribution of observed muon events. Moreover, it would be possible to identify the dip in µ − and µ + events separately, due to the muon charge identification capability of ICAL.
In contrast to the study in Ref. [13] where the ratio between data and unoscillated Monte Carlo (MC) as a function of L/E was used, in this paper we use the ratio between upward-going and downward-going events for the analysis -a data-driven approach. We exploit the fact that the downward-going atmospheric neutrinos and antineutrinos in multi-GeV energy range do not oscillate, and up-down asymmetry reflects the features mainly due to oscillations in upward-going events in these energies. The magnetic field in the ICAL detector enables us to study these distributions in neutrino and antineutrino oscillation channels separately. In this paper, we also discuss the procedure to measure the atmospheric mass-squared difference |∆m 2 32 | and the mixing angle θ 23 using the information from the oscillation dip, and provide the results.
In this paper, we point out for the first time that there is an "oscillation valley" feature in the two-dimensional plane of (E rec µ , cos θ rec µ ) at ICAL. This feature would be recognizable at ICAL due to its sensitivity to muons having energies in the wide range of 1 GeV to 25 GeV and excellent angular resolution at these energies (the muon energy resolution at ICAL is around 10% to 15%, whereas the direction resolution is less than 1 • for multi-GeV muons with directions away from the horizon [14]). We show that, due to this excellent energy and direction resolution, the oscillation valley in two dimensional (E rec µ , cos θ rec µ ) plane would appear prominently at ICAL. We also demonstrate for the first time that the mass-squared difference may be determined using the alignment of the oscillation valley.
There are studies to measure the sensitivity of the ICAL detector for measuring the atmospheric oscillation parameters |∆m 2 32 | and θ 23 , see Ref. [15][16][17]. The difference between these analyses and our study is that Ref. [15][16][17] employ the sophisticated χ 2 method and also take into account possible systematic uncertainties. On the contrary, the focus of our study is to identify the oscillation dip and the oscillation valley, independently in µ − and µ + events, and determine the oscillation parameters from them. This is possible due to the wide range of energy and baselines reconstructed at ICAL with excellent detector properties. We also include statistical fluctuation with 10-year data at ICAL. This paper is organized in the following fashion. In Sec. 2, we present the survival probabilities of neutrinos and antineutrinos as one-dimensional functions of L ν /E ν and as two-dimensional oscillograms in (E ν , cos θ ν ) plane, to pinpoint the oscillation dips and the oscillation valley, respectively. For exploring these features in atmospheric neutrino and antineutrino data separately, we simulate ν µ andν µ interactions at ICAL as discussed in Sec. 3. Next, in Sec. 4, we formulate a data-driven variable, the ratio between upwardgoing and downward-going events (U/D), to use it as an observable for all the analyses in this paper. We propose a novel algorithm to identify the oscillation dip and to measure its "location" in L/E distributions. We present the 90% C.L. range for |∆m 2 32 | expected with 500 kt·yr exposure of ICAL using the calibration curve between the location of the dip and |∆m 2 32 |. Sec. 5 is devoted to discuss our analysis of U/D distributions in (E rec µ , cos θ rec µ ) plane. For identifying the oscillation valley in (E rec µ , cos θ rec µ ) plane of these distributions and to measure the "alignment" of the valley, we suggest a unique algorithm and demonstrate its use. We present the 90% C.L. range of |∆m 2 32 | using 500 kt·yr exposure of ICAL with the help of the calibration curve between the alignment of the valley and ∆m 2 32 . In Sec. 6, we summarize our findings and point out the utility of our novel approach in establishing the oscillation hypothesis in atmospheric neutrino experiments.

Neutrino Oscillation Probabilities
In this section, we discuss the survival probabilities of atmospheric ν µ andν µ , reaching the detector after their propagation through the atmosphere and possibly the Earth. The oscillation probabilities are functions of energy (E ν ) and zenith angle (θ ν ) of neutrino. The net distance L ν traversed by a neutrino (its "baseline") is related to its zenith angle via where R, h, and d denote the radius of Earth, the average height from the surface at which neutrinos are produced, and the depth of the detector underground, respectively. In this study, we take R = 6371 km, h = 15 km, and d = 0 km. Since R h d, the changes in the survival probabilities with small variations in h and d are negligible.
The survival probabilities may be represented in two ways: (i) as one-dimensional functions of L ν /E ν , and (ii) as two-dimensional oscillograms in the plane of energy and zenith angle. Below we discuss these two representations, and the major physics insights obtained via them.
2.1 L ν /E ν dependence of ν µ andν µ survival probabilities Figure 1 shows the survival probabilities for ν µ andν µ as functions of log 10 [L ν /E ν ] for different zenith angles for a benchmark set of oscillation parameters as shown in Table 1. It may be observed that at log 10 [L ν /E ν ] ≈ 2.7, the first oscillation minimum appears for We use the oscillation parameters given in Table 1.   [6][7][8][9]. Normal mass ordering corresponds to m 1 < m 2 < m 3 .
ν µ as well as forν µ . For L ν /E ν values smaller than the location of the first oscillation minimum, the survival probabilities for ν µ andν µ are almost equal, and the overlapping nature of different cos θ ν curves shows that there is no θ-dependence beyond that coming from L ν /E ν . This is an indication that the survival probabilities are dominated by vacuum oscillations in this region. However for higher L ν /E ν values, Earth matter effects introduce an additional θ-dependence beyond that via L ν /E ν , which can be prominently seen for ν µ survival probability in the case of NO. Since matter effects on antineutrinos are not significant in NO, theν µ survival probability does not show this feature.
2.2 Oscillograms in (E ν , cos θ ν ) plane Figure 2 shows the survival probabilities of ν µ andν µ in the (E ν , cos θ ν ) plane, with NO as the neutrino mass ordering. In both the panels of Fig. 2, we can see the prominent dark diagonal band, which represents the minimum survival probability. This oscillation minimum region corresponds to the same "oscillation dip" as that earlier observed in Fig. 1 at log 10 [L ν /E ν ] ≈ 2.7. Its nature in neutrino and antineutrino survival probabilities is almost identical, indicating that this region is dominated by vacuum oscillations. In this study, we will explore whether we can reconstruct this band of oscillation minima (an "oscillation valley") from the atmospheric neutrino data. The part of the oscillogram below plane, using three-flavor neutrino oscillation framework in the presence of Earth matter (PREM profile [18]).
We use the oscillation parameters given in Table 1.
the blue band corresponds to lower energies and longer baselines, where the matter effects are significant for neutrinos (in the NO scenario considered here). For neutrinos, the MSW resonance [19][20][21] can be observed for 6 GeV < E ν < 8 GeV and −0.7 < cos θ ν < −0.5. The parametric resonance [22,23] or oscillation length resonance [24,25] can be seen for 3 GeV < E ν < 6 GeV and cos θ ν < −0.8. The region above the blue band corresponds to vacuum oscillations, and is almost identical for ν µ andν µ survival probabilities.
In the next sections, we demonstrate how these multi-GeV range features in the probability may be reflected in the expected event spectra of neutrinos and antineutrinos in an atmospheric neutrino experiment like ICAL.

Event Generation at ICAL detector
While the analysis in this paper can in principle be performed with any atmospheric neutrino detector, we use the proposed ICAL detector at the India-based Neutrino Observatory [26,27] for our simulations. There are two main reasons for this. (i) ICAL would be able to distinguish neutrinos from antineutrinos, and hence independent analyses for these two channels may be carried out. (ii) ICAL would be able to reconstruct the energy as well as the direction of µ − and µ + (produced in the charged-current interactions of ν µ and ν µ , respectively) with a good precision, which will be crucial for the analysis.
ICAL will be composed of 50 kt magnetized iron as the target, and glass resistive plate chambers (RPCs) as the sensitive detector elements. Iron plates of 5.6 cm thickness, and around 30000 glass RPCs, will be stacked alternately in three modules, each of dimension 16 m (L) × 16 m (W) × 14.5 m (H). The magnetic field in the detector volume will be around 1.5 T. Charged particles passing through the RPCs will ionize the gas inside them, thereby inducing an electrical signal in the X-and Y-pickup strips and providing the X-Y coordinates of the charged particle. The layer number provides the Z coordinate. In the charged-current (CC) interactions of neutrinos, the muon, and hadrons in the final state will be distinguished from their track-like and shower-like features respectively. The time resolution of less than 1 ns for each RPC layer [28][29][30] would enable ICAL to distinguish between upward-going and downward-going events with a high confidence. The direction of bending of the charged muon track will enable the identification of its charge, and hence the identification of whether the original interacting particle was a neutrino or antineutrino.
To simulate the neutrino interactions in the ICAL detector, we use the MC neutrino generator NUANCE [31] and the neutrino flux calculated for the INO site [32,33]. In this context, it is worthwhile mentioning that the studies of atmospheric neutrino oscillations [27] for the ICAL detector have been performed using the neutrino fluxes calculated for the Kamioka site. The neutrino fluxes vary from place to place due to the different geomagnetic field of the Earth, and there are important differences in the fluxes at Kamioka and at the proposed site of INO at Theni district of Tamil Nadu, India. Theni is close to the region having the largest horizontal component (≈ 40 µT ) of the Earth's magnetic field. In comparison to this, the horizontal component of the geomagnetic field at Kamioka is ≈ 30 µT . As a result, the neutrino flux at the INO site is smaller by approximately a factor of 3 at low energies. However, for neutrino energies around 10 GeV, the fluxes at INO and Kamioka sites are almost the same [33]. In this study, we use the neutrino fluxes given in [32,33], taking into account the 1 km rock coverage of the mountain. Also, to take into account the effect of solar modulation on the neutrino fluxes, we use the fluxes with high solar activity (solar maximum) for half the exposure and the fluxes with low solar activity (solar minimum) for the remaining half. These fluxes are given in the website [34].
In this work, we focus on the CC interactions of muon neutrinos and antineutrinos at the ICAL detector. These may be a result of the unoscillated original ν µ /ν µ , or original ν e /ν e which have oscillated to ν µ /ν µ . The three-flavor neutrino oscillations in the presence of matter (PREM profile [18]) are incorporated using the re-weighting algorithm as in [15][16][17]. The events are then folded with the detector properties as obtained from the GEANT4 simulation for muons at ICAL [14]. The reconstructed energy and zenith angle of the charged muon (E rec µ and θ rec µ ) are taken to be the observable quantities, and no attempt is made at reconstructing the neutrino energy or incoming direction. The energy deposited at the detector by hadrons in the final state of neutrino interaction can be reconstructed at ICAL [35], but we do not use the hadron energy information in this study. Indeed, we use a data-driven approach in this paper, where we only use the ratios of upward-going and downward-going charged muons in various energy and direction bins, thus reducing the impact of uncertainties in the calculations of absolute values of the original neutrino flux and its directional dependence. Note that the upward-going vs. downward-going events may lead to some ambiguity when the events are near-horizon, i.e., −0.2 < cos θ rec µ < 0.2 (or 73 km < L rec µ < 2621 km). For these events, there is a significant change in the neutrino path length for a small variation in the estimated neutrino arrival direction. However, the detector response of ICAL to such events is anyway very poor, owing to the horizontally stacked structure, and it is observed that these events (corresponding to log 10 [L rec µ /E rec µ ] in the range of 1.5 -2.0) do not affect the analysis 1 . Note that the Super-K collaboration selected only the "high-resolution" events in their L/E analysis [36,37]. In particular, their analysis rejected neutrino events near the horizon, as well as low-energy events where the large scattering angles would have led to large errors in the reconstruction of neutrino direction. The low efficiency of ICAL for the near-horizontal event, and our analysis threshold of 1 GeV for muons, automatically incorporates both these filters. However, a difference between our analysis and the L/E analysis of Super-K [36] also needs to be pointed out. While the analysis in [36] is in terms of the inferred L ν and E ν of neutrinos, our analysis is in terms of the L rec µ and E rec µ of muons. In Sections 4 and 5, we discuss the event distributions as the one-dimensional functions of L rec µ /E rec µ , and in the two-dimensional (E rec µ , cos θ rec µ ) plane, respectively.
In this section, we analyze the expected event distributions of µ − and µ + as functions of L rec µ /E rec µ , with an MC event sample corresponding to 1000 years, and simulated event sample corresponding to 10 years, for ICAL. The 1000-year MC sample will correspond to the expected observations if an unlimited amount of data were available, and analysis using this sample will guide our data-driven algorithms. It will also lead to the calibrations of oscillation parameters like |∆m 2 32 | and sin 2 θ 23 . On the other hand, the analysis with the 10-year sample would give an idea of the effects due to statistical fluctuations, an important consideration for low-statistics experiments like those with atmospheric neutrinos. To this end, we will analyze 100 independent sets of the 10-year samples.  Table 2: Total number of upward-going (U) and downward-going (D) µ − and µ + events expected at the 50 kt ICAL detector in 10 years (total exposure of 500 kt·yr). We use the oscillation parameters given in Table 1.
The number of events expected in 10 years at the 50 kt ICAL are shown in Table 2. Here, U is the number of upward-going events (cos θ rec µ < 0), while D is the number of downward-going events (cos θ rec µ > 0). The U/D ratio is defined as where N (E rec µ , cos θ rec µ ) is the number of events with energy E rec µ and the zenith angle θ rec µ . We associate the U/D ratio as a function of cos θ rec µ with the upward-going events, i.e., with cos θ rec µ < 0. We can also associate this U/D asymmetry with the L rec µ corresponding to these upward-going events, which is related to cos θ rec µ by Eq. (2.1). The ratio U/D may be taken to be a proxy for the probability P (ν µ → ν µ ) since in the ICAL detector, around 98% events are contributed from the oscillation channel ν µ → ν µ , and only ∼ 2% of total events come from the ν e → ν µ channel for the benchmark values of the oscillation parameters shown in Table 1. With the use of this ratio, the analysis becomes less dependent on the uncertainties in the absolute neutrino flux values. Below, we shall analyze the event distributions and U/D ratios for µ − and µ + events separately, as the functions of L rec µ /E rec µ .  Table 1.

Events and U/
The upper panels of Fig. 3 present the distributions of µ − (left panel) and µ + (right panel) events for the 1000-year MC sample. We take events with the reconstructed muon energy in the range (E rec µ ) min = 1 GeV to (E rec µ ) max = 25 GeV. Since L rec µ is in the range of 15 km to 12757 km, the minimum and maximum values of log 10 respectively. Note that for upward-going events, the minimum value of log 10 [L rec µ /E rec µ ] is 1.2, since L rec µ (min) for upward-going events is 437 km. We have binned the data in 100 bins, uniform in log 10 [L rec µ /E rec µ ]. The number of upward-going events is clearly less than that of the downward-going events, since the upward-going ν µ have traveled larger distances and have had a larger chance of oscillating to the other neutrinos flavors. The number of µ + events is less than that of µ − events mainly because, at these energies, the cross-section of antineutrinos is smaller than that of neutrinos by a factor of approximately 2.
For the lower panels of Fig. 3, we divide the range of upward-going neutrinos (cos θ rec µ < 0), i.e. log 10 [L rec µ /E rec µ ] = 1.2 -4.1, in 66 uniform bins. A downward-going event with a given cos θ rec µ > 0 value is then assigned to the bin corresponding to − cos θ rec µ . The U/D ratio in each bin is then calculated by dividing the number of upward-going events by the number of downward-going events, and plotted with the corresponding value of log 10 [L rec µ /E rec µ ]. Two major dips in the U/D ratio may be observed in the figure, in both the µ − and µ + channels. These dips occur at log 10 [L rec µ /E rec µ ] ∼ 2.75 and ∼ 3.2, respectively. In this study, we focus on the observation of the first dip.
Note that, while the two dips discussed above correspond to the first two minima in the survival probability as shown in Fig. 1, they are not exactly at the same location. This is because of the crucial difference that, while Fig. 1 uses the actual values of energy and zenith angle of the atmospheric neutrino, Fig. 3 uses the reconstructed values of energy and zenith angle of the muon produced from the CC interaction of that neutrino. In CC deep inelastic scattering of neutrino, hadrons in the final state take away some energy of the incoming neutrino, resulting in E rec µ < E ν . As a result, the dips move towards higher values of L rec µ /E rec µ . Since the average inelasticities 2 in the antineutrino events ( yν ≈ 0.3) are smaller than that in the neutrino events ( y ν ≈ 0.45) in the multi-GeV energy range [17], the shift of the dip is smaller in the case of µ + than µ − .
Note that the dips in Fig. 3 are shallower and broader than those in Fig. 1, where the dips reach all the way to the bottom, that is P (ν µ → ν µ ) ≈ 0, with sin 2 θ 23 = 0.5. The reason for this smearing lies in the difference between the momenta of the neutrino and muon, as well as the limitations on the muon momentum reconstruction in the detector. It may also be observed that the smearing is more in µ − events as compared to the µ + events in Fig. 3. This difference is due to the broader spread in the inelasticities of neutrino interactions than that in antineutrino interactions [17]. As a result of this spread, the dips in µ − events get smeared more, and hence become shallower, as compared to the dips in µ + events.  . The height of shaded boxes represents expected uncertainties in the U/D ratio. The exposure is taken to be 10 years of 50 kt ICAL, i.e. 500 kt·yr. We use the oscillation parameters given in Table 1.

Events and U/D ratio using 10-year simulated data
In this section, we discuss the expected event distributions as functions of L rec µ /E rec µ at the ICAL detector, with an exposure of 500 kt·yr. In order to take into account the statistical fluctuations, which are expected to have a significant impact on the accuracy of the results, we take 100 independent sets, and calculate the mean as well as root-mean-square (rms) deviations of relevant quantities.
We distribute the events in bins corresponding to the log 10 [L rec µ /E rec µ ] values, as described in Sec. 4.1. However, since the number of events available here is much smaller than that for 1000 years, we choose a non-uniform binning scheme shown in Table 3 such that typically, all the bins will have at least 10 down-going events. We have a total of 34 bins of log 10 [L rec µ /E rec µ ] in the range of 0 to 4. The distribution of the number of events in these bins is shown in the top panels of Fig. 4. The uncertainties due to statistical fluctuations, calculated as the rms deviation obtained from 100 independent simulated data sets, are also shown.
The lower panels of Fig. 4 show the U/D ratio in log 10 [L rec µ /E rec µ ] bins, using the same procedure outlined in Sec. 4.1. We choose the range of log 10 [L rec µ /E rec µ ] to be in the range 2.0 -4.0, which is the most interesting range. The fluctuations shown in the figure are the rms deviations obtained from the distributions of the U/D ratio in 100 independent simulated data sets. It is observed that, while the first dip is quite prominent, the second dip is lost in the statistical fluctuations, and due to the broad binning scheme that we have to employ owing to a small number of events. While this does not rule out the identification of the second dip using more efficient algorithms, in this work, we shall focus on the first dip. Henceforth, when we refer to the dip position, we refer to the position of the first dip. It may be observed that the dip in the µ + channel is deeper than that in the µ − channel similar to 1000-year MC sample as described in Sec. 4.1.

Identifying the dip with 10-year simulated data
The left panel of Fig. 5 shows the U/D ratios as functions of log 10 [L rec µ /E rec µ ] for five independent simulated data sets with 500 kt·yr exposure of ICAL. Clearly, the statistical fluctuations are significant, and may lead to the misidentification of the dip position. An identification of the dip position should not correspond simply to the lowest value of the U/D ratio, but also be guided by the values of the U/D ratios in surrounding bins. In order to achieve this, we use the dip-identification algorithm, which is a data-driven approach as described below.
To start with, we consider the region corresponding to the range log 10 [L rec µ /E rec µ ] = 3.2 -4.0 as a single bin, since in this region, we expect the oscillations to be quite rapid, leading to the U/D ratio averaging out to a constant value. Let the measured value of the U/D ratio in this bin be R 0 . From the simulations done with 100 independent sets, we have found that the statistical fluctuations in data sets give rise to a rms deviation of ∆R 0 ≈ 0.02 and 0.03 in this ratio for µ − and µ + respectively. We take this fluctuation into account, and start with an initial ratio threshold of R th ≡ R 0 − 2∆R 0 , shown by the black dashed line in the right panel of Fig. 5. All the bins with measured U/D ratio less than R th form the initial candidates for the dip position. However, all these bins need not be a part of the actual dip, due to fluctuations. In order to identify the actual bins surrounding the dip, we try to find the cluster of consecutive bins, all of which have U/D ratio less than that in all the other bins. This is achieved by lowering the value of R th till all the bins with U/D ratio less than R th are contiguous. This final value of R th is denoted by the blue dashed line in the right panel of Fig. 5.
Once the dip-identification algorithm has thus identified a cluster of contiguous bins as the dip region, we fit the U/D ratios in these bins with a parabola as shown in the right panel of Fig. 5. The value of log 10 [L rec µ /E rec µ ] corresponding to the lowest U/D ratio obtained from this fit, denoted by x min , would be identified with the "location" of the dip. 2 32 | and θ 23 from the dip and the U/D ratio The first dip in the survival probability P (ν µ → ν µ ), in the approximation of two-flavor oscillations in vacuum, appears at L ν (km)/E ν (GeV) = π/(2.54 |∆m 2 32 |(eV 2 )) (see Eq. 1.1). However as discussed in Sec. 4.1, the location of the dip x min in the L rec µ /E rec µ distribution would have a complicated dependence on the three flavor neutrino oscillation, matter effect, neutrino fluxes, CC cross-sections, inelasticities of events, and detector response. However, we can obtain a calibration of |∆m 2 32 | with the location of the dip, using the 1000-year MC sample. The binning scheme used is the same as in Sec. 4.3. The blue points in the upper panels of Fig. 6   different values of |∆m 2 32 |. These points are observed to lie close to a straight line, and we can draw a calibration curve that would enable us to infer the actual value of |∆m 2 32 |, given the x min value determined from the data. Note that this calibration curve also has a dependence on the analysis procedure, including the binning scheme and bin-identification algorithm.

|∆m
While the calibration curve allows a one-to-one correspondence between the actual |∆m 2 32 | value and the location of the dip for sufficiently large data (like 1000 years of exposure), practically speaking the data available is going to be limited. The statistical fluctuations introduced due to this will lead to uncertainties in the determination of |∆m 2 32 | via this method. To estimate these uncertainties, we fix a value of |∆m 2 32 | (true) = 2.46 × 10 −3 eV 2 to generate 100 independent simulated data sets with an exposure of 500 kt·yr each, and determine the location of the dip in each of them. The distribution of the dip locations thus obtained allows us to determine 90% C.L. allowed regions for the dip location, and hence the 90% C.L. allowed regions for the calibrated |∆m 2 32 | values. These allowed regions are represented in the upper panels of Fig. 6 by gray bands. The figure indicates that the 90% C.L. allowed range for |∆m 2 32 | from µ − events is (2.15 -2.65)×10 −3 eV 2 , while that from µ + data is (2.22 -2.80)×10 −3 eV 2 .
We can follow the same procedure to determine the value of the mixing angle θ 23 from the data. In principle, this is related to the depth of the dip, and such a calibration could be obtained using the 1000-year MC events. However, the statistical fluctuations in the 500 kt·yr simulated data are observed to give rise to large uncertainties. On the other hand, it is found that the total U/D ratio of all events with L rec µ /E rec µ in the range of 2.2 -4.1 gives a much better estimate of θ 23 . In the lower panels of Fig. 6, we show in blue the calibration for sin 2 θ 23 with the measured total U/D ratio, obtained using the 1000-year MC sample for several sin 2 θ 23 (true) values. We also show the 90% C.L. allowed values of sin 2 θ 23 , inferred using the 100 independent simulated data sets with 500 kt·yr exposure each, at sin 2 θ 23 (true) = 0.5. The expected allowed range for sin 2 θ 23 at 90% C.L. is (0.44 -0.64) from µ − events, and (0.37 -0.63) from µ + events.
Note that we do not claim or demand that our analysis gives better results than those obtained with the complete χ 2 analysis done for ICAL in Ref. [16,17]. The χ 2 analyses take into account the complete event spectra, while our focus is on the identification of the dip, a feature that would reconfirm the oscillation paradigm for neutrinos. The observation that the inferred |∆m 2 32 | range is comparable to the range expected from the complete χ 2 analysis points to the conclusion that most of the information about the |∆m 2 32 | is concentrated in the location of the U/D dip. In our method, that uses the U/D ratio, the systematic uncertainties like flux normalization, cross section, overall detection efficiency, and energy dependent tilt error get canceled out automatically.

Oscillation valley in the
The oscillation dip discussed in the last section was a dip in the U/D ratio as a function of L rec µ /E rec µ . The dependence on L rec µ /E rec µ was motivated from the approximate form of the neutrino oscillation probability (with two flavors in vacuum). However, if the detector can reconstruct both L µ (hence cos θ µ ) and E µ accurately, then one can go a step ahead and look at the distribution of the U/D ratio in the (E µ , cos θ µ ) plane. We perform such an analysis in this section, and find that such a distribution can have many interesting features with physical significance, which may be identifiable with sufficient data. In particular, we point out an "oscillation valley" corresponding to the dark diagonal band in Fig. 2, whose nature can provide stronger tests for the oscillation hypothesis, albeit only with a large amount of data. On the other hand, we show that the identification of the oscillation valley, and the determination of |∆m 2 32 | based on its "alignment", is possible with 500 kt·yr of simulated data at ICAL.  Note that the U/D(E rec µ , cos θ rec µ ) is defined only for cos θ rec µ < 0 (see Eq. 4.1).

Events and U/D ratio using 1000-year Monte Carlo simulation
In this section, we discuss the distributions of events and the U/D ratio, in the plane of reconstructed energy E rec µ and zenith angle cos θ rec µ of the µ − and µ + events, for a 1000year MC sample. The main reason for considering such a huge exposure here is not to miss those features, which survive in spite of our not using E ν and cos θ ν directly, but which may disappear due to the limitation of statistics.
The quantities shown in Fig. 7 are binned based on the reconstructed values of E rec µ and cos θ rec µ . The cos θ rec µ range of -1.0 to 1.0 has been divided into 80 bins of equal width. We have a total of 84 bins in E rec µ in the range of 1 GeV to 25 GeV. The bin width is not uniform in E rec µ , since at higher energies the number of events decreases rapidly. The first µ bins, in the range of 1.0 -5.5 GeV, have a width of 0.1 GeV each, and the last 39 bins, in the range of 5.5 -25 GeV, have a width of 0.5 GeV each. The upper panels of Fig. 7 present the distributions of µ − (left panel) and µ + (right panel) event density. The abrupt lowering of the event density for −0.2 < cos θ rec µ < 0.2 is due to the poor reconstruction efficiency of the ICAL detector for the events coming near the horizontal direction. Note that this feature was absent in Fig. 3 since it was distributed over a large L rec µ /E rec µ range. In lower panels of Fig. 7, we show the U/D ratio in each bin, in the (E rec µ , cos θ rec µ ) plane. Note that the U/D(E rec µ , cos θ rec µ ) is defined only for cos θ rec µ < 0 (See Eq. 4.1). This also makes the features in Fig. 7 resemble those in Fig. 2, thus enabling their understanding in terms of the survival probabilities shown in the oscillograms. The lower panels may be observed to have many interesting features. The most prominent one is of course the light/ dark blue band, corresponding to U/D 0.25, extending diagonally from (E rec µ , cos θ rec µ ) ≈ (5 GeV, −0.3) to (E rec µ , cos θ rec µ ) ≈ (25 GeV, −1.0). This is the band with the lowest values of U/D ratio in this plane, and we shall henceforth refer to it as the "oscillation valley". Clearly, this valley is deeper for µ + , since the number of dark blue bins with U/D < 0.15 is seen to be larger for µ + . This is consistent with the observation of a deeper oscillation dip in µ + , as discussed in Sec. 4.1. The main reason for this is the smaller inelasticity of antineutrino CC events. The differences observed between the left and right panel may be attributed to the Earth matter effects, which affect the survival probabilities of neutrinos and antineutrinos differently for a given neutrino mass ordering.

Events and U/D ratio using 10-year simulated data
In the upper panels of Fig. 8, we show the expected event distributions of µ − (left panel) and µ + (right panel) events, for 500 kt·yr exposure of ICAL, in the plane of reconstructed muon energy (E rec µ ) and reconstructed muon zenith angle (cos θ rec µ ). The cos θ rec µ range of -1.0 to 1.0 is divided into 20 uniform bins of width 0.1 each. The binning in E rec µ is non-uniform, and is such that there are at least a few number of events in each bin (except in the largest energy bins, where this may not be possible due to the smaller neutrino flux at higher energies). The E rec µ range of 1 -25 GeV has been divided into 16 non-uniform bins, as given in Table 4 Table 4: The binning scheme considered for E rec µ and cos θ rec µ for µ − and µ + events of 10year simulated data while we show event distribution and U/D plots in the (E rec µ , cos θ rec µ ) plane. This binning scheme is used in analysis of oscillation valley as well.  for µ − (left panels) and µ + (right panels) in the plane of cos θ rec µ and E rec µ . The mean is calculated from 100 independent data sets of exposure 500 kt·yr each of ICAL. The number of events written in white are rounded to the nearest integer. We use the oscillation parameters given in Table 1. Note that the U/D(E rec µ , cos θ rec µ ) is defined only for cos θ rec µ < 0 (see Eq. 4.1).
In order to take care of the statistical fluctuations that will clearly be significant with 500 kt·yr data, we generate 100 independent data sets, each with this exposure. In the upper panels of Fig. 8, we show the mean number of events obtained from these 100 data sets. Similarly, in the lower panels of Fig. 8, we calculate the U/D ratios in each bin for the above 100 sets, and present their mean values in the figure, using appropriate colors. Note that the U/D(E rec µ , cos θ rec µ ) is defined only for cos θ rec µ < 0 (see Eq. 4.1).
Clearly, while many of the features corresponding to matter effects seem to get diluted due to the coarser bins and statistical fluctuations, the oscillation valley along the diagonal survives. In the following sections, we shall provide a procedure for identifying this oscillation valley and extracting information about the oscillation parameters from it.

Identifying the valley with 10-year simulated data
We now describe the methodology adopted for the identification of the oscillation valley in the (E rec µ , cos θ rec µ ) plane, and the determination of its alignment. To this end, we use the binning described in sec. 5.2 for the 500 kt·yr data set.
• Since the events with log 10 [L rec µ /E rec µ ] > 3.1 are susceptible to significant matter effect and rapid oscillations, the information from these bins cannot be directly connected to the dip corresponding to the vacuum oscillations. As far as the region log 10 [L rec µ /E rec µ ] < 2.2 is concerned, the oscillation are not developed yet. Therefore, we only choose the region 2.2 < log 10 [L rec µ /E rec µ ] < 3.1 for our analysis, where vacuum oscillations are expected to dominate, and to be resolvable.
• We then analyze each horizontal "energy strip" in this region, i.e. the U/D ratio for all cos θ rec µ , as shown in the left panel of Fig. 9. We perform a parabolic fit to the points in this strip that are in the allowed region as specified above. The equation of the parabola to be fitted is y = y E + a E (x − x E ) 2 , where y is the U/D ratio while x is the value of cos θ rec µ . The parameters (x E , y E ) correspond to coordinates (cos θ rec µ , U/D) at the bottom of the parabola, while the parameter a E describes the "width" of the parabola. Higher the value of a E , more prominent the dip in that energy strip.
• Once the parabolic fits in all energy strips are obtained, all the points (E rec µ = E, cos θ rec µ = x E ) are plotted on the (E rec µ , cos θ rec µ ) plane, as shown in the right panel of Fig. 9. A straight line, passing through the origin (E rec µ = 0, cos θ rec µ = 0), is fitted to these points. For this fitting, each of the points is given a weight of √ a E N E , where N E is the number of downward-going events in the analysis region (between the solid and dashed black lines in the right panel Fig. 9) of the energy strip. The points lying outside the analysis region are not taken into account in the fit. The best-fit line is thus obtained of the form y = m · x, where y represents E rec µ , and x represents cos θ rec µ . The slope m of this line quantifies the "alignment" of the valley.
Note that the above procedure is motivated by our expectation that the valley would correspond to the value of L rec µ /E rec µ which would correspond to the dip in the onedimensional L/E analysis described in Sec. 4.3. Our analysis is rather simple, focused only on the identification of the valley and on determining its alignment. With more data and a refined analysis, one maybe able to use the goodness of the fit to the line, or the intercepts of the line on either axes, for testing the oscillation framework in more detail. We show here that 500 kt·yr of exposure would be enough to locate the valley and quantify its alignment even with our simplified treatment. In the following section, we shall also show that the value of |∆m 2 32 | may be calibrated using the alignment of the valley.  For the same reasons described in Sec. 4.4, we expect the alignment of the valley, as determined through the procedure outlined in the previous section, to lead to the value of |∆m 2 32 |. Indeed, for upward-going events, L rec µ is approximately proportional to cos θ rec µ , so that the alignment of the valley, i.e. E rec µ (GeV)/ cos θ rec µ , is approximately equal to (5.08/π) |∆m 2 32 |(eV 2 ) R(km). We calibrate the actual value of |∆m 2 32 | with the alignment of the valley using the same procedure as in Sec. 4.4. The 1000-year MC samples, with different true values of |∆m 2 32 |, are used to locate the blue calibration points, and hence the blue calibration curve, in Fig. 10, the statistical fluctuations with 500 kt·yr data are taken into account by generating 100 independent data sets with a fixed value of |∆m 2 32 | (true) = 2.46×10 −3 eV 2 , and showing the 90% C.L. allowed values for the alignment of the valley, and hence the 90% C.L. allowed values of the calibrated |∆m 2 32 |. These allowed regions are represented in Fig. 10 by gray bands. The figure indicates that the 90% C.L. allowed range for |∆m 2 32 | from µ − events is (2.31 -2.92)×10 −3 eV 2 , while that from µ + data is (2.19 -2.81)×10 −3 eV 2 .

|∆m
So far, our analysis has been based upon 500 kt·yr of simulated data. However, ICAL should be able to identify the oscillation dip/oscillation valley even with half the exposure, and may provide some early measurement of |∆m 2 32 | with data in the multi-GeV range. We show some results with 250 kt·yr exposure in Appendix A.

Summary and Concluding Remarks
Atmospheric neutrino experiments access a large range of energies and baselines for neutrino oscillations, hence they are capable of testing basic features of the neutrino flavor conversion framework over a large parameter space. In this work, we study the potential of the ICAL detector to visualize the L/E dependence of survival probabilities in neutrino and antineutrino channels separately, by distinguishing µ − and µ + events. For this, we adopt a data-driven approach, and use the reconstructed muon energy and zenith angle, as opposed to the inferred neutrino energy and zenith angle. Further, we consider the ratio of upward-going (U) and downward-going (D) events, instead of the ratio of observed oscillated events and simulated unoscillated events. The reconstructed L rec µ /E rec µ distributions of the U/D ratio for 1000-yr MC indicate that the feature of the first oscillation dip in this ratio in muon neutrinos (and antineutrinos) is still preserved in the detected µ − and µ + . We demonstrate that this first oscillation dip can be clearly identified at ICAL with an exposure of 500 kt·yr, i.e. 10 years running of the 50 kt ICAL.
We develop a novel dip-identification algorithm to identify the oscillation dip and measure the value of |∆m 2 32 |, using the information from the first oscillation dip in the L rec µ /E rec µ distributions of the U/D ratio, in µ − and µ + events separately. This algorithm finds a contiguous set of L rec µ /E rec µ bins with the lowest U/D values, and uses the combined information in all of them to determine the dip location precisely. The value of |∆m 2 32 | is then calibrated against the location of the dip using the 1000-yr MC, and the statistical uncertainties expected with the 10-year data are estimated using simulations of 100 independent data sets. We summarize the expected 90% C.L. ranges of |∆m 2 32 | in Table 5. We also calibrate θ 23 , however, here simply the ratio of the total number of U and D events is found to be a good choice to calibrate against. We find that, when the true value of sin 2 θ 23 is 0.50, the U/D ratio would determine the 90% C.L. range of its value to be (0.44 -0.64) with µ − events, and (0.39 -0.61) with µ + events.  Table 5: The 90% C.L. allowed ranges for expected for |∆m 2 32 |, from the analyses using L rec µ /E rec µ and (E rec µ , cos θ rec µ ) distributions of the U/D ratio, with 500 kt·yr exposure. The true value of |∆m 2 32 | is taken to be 2.46 × 10 −3 eV 2 , and the other oscillation parameters are kept fixed at the benchmark values given in Table 1.
In this paper, we point out for the first time that the identification of an "oscillation valley" feature is possible, in the distribution of the U/D ratio in the plane of (E rec µ , cos θ rec µ ). The oscillation dip observed in the one-dimensional L rec µ /E rec µ analysis above now appears as an oscillation valley in the two-dimensional (E rec µ , cos θ rec µ ) plane. We go on to develop an algorithm for finding the alignment of the valley, i.e. the slope of the best-fit line to the valley in the (E rec µ , cos θ rec µ ) plane, and show that it serves as a good proxy for |∆m 2 32 |. As in the oscillation dip analysis, we calibrate |∆m 2 32 | against the alignment of the valley using 1000-yr MC data, and estimate the uncertainties with 10-year data using 100 independent data sets. The expected 90% C. L. ranges of |∆m 2 32 | are summarized in Table 5. We also show the expected results for 90% C.L. ranges for |∆m 2 32 |, obtained with half the data (5 years data, or an exposure of 250 kt·yr), in Appendix A. This shows that preliminary results with the oscillation dip and oscillation valley analyses may already start appearing within the first few years of ICAL. Note that the aim of this exercise is not to compare the precisions achieved by different methods -indeed, this study does not aim to compete with the conventional χ 2 method for precision. Our aim is to present an approach ("oscillation dip") that has proved useful in the past for establishing the oscillation hypothesis and rejecting some alternative hypotheses, and to go one step beyond it, by performing the "oscillation valley" analysis in the two-dimensional (E rec µ , cos θ rec µ ) plane. This two-dimensional analysis would open up avenues for testing the oscillation framework from more angles -for example, by looking for effects of new physics on the shape, alignment, width, or depth of the valley -which would be addressed in future work. An important feature of our analysis is that it is data-driven, i.e. we only use the reconstructed observables in our analysis.
Though we have presented our results for ICAL, the same procedure can be used for any present or upcoming atmospheric neutrino experiment that has access to a large range of energies and baselines. An atmospheric neutrino experiment can perform the analysis based on the oscillation dip as well as the oscillation valley. Of course, if the detector cannot identify the muon charge, the data from neutrino and antineutrino channels will have to be combined, which may dilute some of the observable features. However, the principle of performing the oscillation valley analysis in the (E rec µ , cos θ rec µ ) plane -going one step ahead of the L rec µ /E rec µ analysis -may be adopted in any atmospheric neutrino experiment. Note that while the fixed-baseline experiments may in principle be able to identify the L ν /E ν oscillation dip, they do not have access to the (E ν , cos θ ν ) plane, and hence they cannot perform the oscillation valley analysis.
Even then, certain features of ICAL make it uniquely suited for such an analysis. It will likely be the first detector to identify the oscillation dip separately for the neutrino and antineutrino channels, which will also be a test for CPT violation in neutrino oscillations. Further, as far as the oscillation valley analysis goes, a crucial requirement for the capability to identify the alignment of the valley is enough number of points to determine the best-fit line as shown in the right panel of Fig. 9. This needs a large energy range to which the detector is sensitive (large number of y-bins), and an excellent angular resolution (large number of x-bins). ICAL is thus uniquely suited for the oscillation valley analysis, for providing an orthogonal approach to establish the nature of neutrino oscillations, and hence for making the neutrino oscillation picture more robust.