Invited review: Physics potential of the ICAL detector at the India-based Neutrino Observatory (INO)

Abstract

The upcoming 50 kt magnetized iron calorimeter (ICAL) detector at the India-based Neutrino Observatory (INO) is designed to study the atmospheric neutrinos and antineutrinos separately over a wide range of energies and path lengths. The primary focus of this experiment is to explore the Earth matter effects by observing the energy and zenith angle dependence of the atmospheric neutrinos in the multi-GeV range. This study will be crucial to address some of the outstanding issues in neutrino oscillation physics, including the fundamental issue of neutrino mass hierarchy. In this document, we present the physics potential of the detector as obtained from realistic detector simulations. We describe the simulation framework, the neutrino interactions in the detector, and the expected response of the detector to particles traversing it. The ICAL detector can determine the energy and direction of the muons to a high precision, and in addition, its sensitivity to multi-GeV hadrons increases its physics reach substantially. Its charge identification capability, and hence its ability to distinguish neutrinos from antineutrinos, makes it an efficient detector for determining the neutrino mass hierarchy. In this report, we outline the analyses carried out for the determination of neutrino mass hierarchy and precision measurements of atmospheric neutrino mixing parameters at ICAL, and give the expected physics reach of the detector with 10 years of runtime. We also explore the potential of ICAL for probing new physics scenarios like CPT violation and the presence of magnetic monopoles.

Appendices

Appendix A. Neutrino fluxes at the Kamioka and INO sites

In this section, we briefly discuss the preliminary atmospheric neutrino flux calculations corresponding to the INO site, and compare them with those for the Kamioka site. This work is based on refs [5] and [100]. The primary cosmic ray flux model based on AMS [169] and BESS [170,171] data has been used, and the hadronic interactions have been implemented with DPMJET-III [172] above 32 GeV, and JAM below 32 GeV. For the propagation of cosmic rays in the atmosphere, the model NRLMSISE-00 [173], which takes into account the temperatures and densities of the atmosphere’s components and takes care of the position dependence and the time variation in a year, is used. The geomagnetic field model used is IGRF2010 [174]. The atmospheric neutrino flux is obtained using a three-dimensional scheme below 32 GeV, and a one-dimensional scheme above that. The flux calculated in both the schemes agree with each other at 32 GeV [5].

Figure A.1 shows the calculated atmospheric neutrino flux averaged over one year by folding over all directions and summing over all types of neutrinos ($\nu _{e},\bar \nu _{e},\nu _{\mu },\bar \nu _{\mu }$) for the Super-Kamiokande, INO and South Pole sites. These results, obtained with the NRLMSISE-00 model, agree well with those obtained from the US-standard76 [175] model. However, the fluxes with the NRLMSISE-00, as used here, have an advantage in the study of seasonal variations, which can be appreciable at the INO site.

Figure A.1

Atmospheric neutrino flux averaged over all directions and summed over $\nu _{e}\,{+}\,\bar \nu _{e}\,{+}\,\nu _{\mu }\,{+}\,\bar \nu _{\mu }$, as a function of neutrino energy for Super-Kamiokande, INO, and South Pole sites. This is for all the energy ranges of the calculation. The fluxes calculated with US-standard76 are also plotted in dashed line for the Kamioka and the INO sites.

It is observed that the total flux at INO is slightly smaller than that at Kamioka at low energies (E ≲ 3 GeV), but the difference becomes small with the increase in neutrino energy. It may be noted that this is true only for the angle-integrated fluxes. Figure A.2 shows the zenith angle dependence (integrating over all azimuthal angles) of fluxes at Kamioka and at the INO site at two values of energy. It is found that at 1 GeV, there are large up–down asymmetries in the fluxes at the INO site; the upward-going flux is larger than the downard-going one. These asymmetries decrease with the increase in neutrino energy and almost disappear at 10 GeV.

Figure A.2

The zenith angle dependence of atmospheric neutrino flux at E = 1.0 GeV (left) and E = 3.2 GeV (right), averaged over all azimuthal angles calculated for the Super-Kamiokande (top) and the INO (bottom) sites. Here 𝜃 is the arrival direction of the neutrino, with $\cos \theta \,=\,1$ for vertically downward-going neutrinos and $\cos \theta \,=\,-1$ for vertically upward-going neutrinos.

The results for the atmospheric neutrino fluxes as functions of the azimuthal angle ϕ, in five zenith angle bins, are shown in figure A.3 for Kamioka and in figure A.4 for the INO site. These results are presented for ν _e , $\bar \nu _{e}$, ν _μ , and $\bar \nu _{\mu }$, for the (anti)neutrino energy of 3.2 GeV. It is observed that even for such a high energy, the variation of the atmospheric neutrino flux has a complex structure. This is a result of the rigidity cut-off and muon bending in the geomagnetic field. At the INO site, the horizontal component of the geomagnetic field is ∼40 μ T, larger than that at Kamioka, where it is ∼30 μ T, and so the azimuthal angle dependence is also more complex. This complex azimuthal angle dependence continues even above 10 GeV for the near horizontal directions.

Figure A.3

The azimuthal angle dependence of atmospheric neutrino flux, averaged over zenith angle bins calculated for the Super-Kamiokande site for ν _e at E = 3.2 GeV.

Figure A.4

The azimuthal angle dependence of atmospheric neutrino flux, averaged over zenith angle bins calculated for the INO site for ν _e at E = 3.2 GeV.

The atmospheric electron neutrino flux also shows a rapid variation in the azimuth angles, but the statistical errors in the simulations of production of atmospheric neutrinos are still large, and more statistics is needed for a better understanding. As the INO site is close to the equator, the seasonal variation is also expected, the calculation for which is also in progress. These updates are expected to be reported after the accumulation of sufficient statistics.

Detailed physics re-analyses with the fluxes at the INO site need to be carried out in order to determine the final physics potential of ICAL. However, the effect of the change of flux may be estimated by comparing the number of events calculated using the fluxes at Kamioka and at the INO site. The comparison of the number of μ ⁻ and μ ⁺ events at these two sites, as a function of the muon energy, is shown in figure A.5. For the sake of this sample comparison, we have considered charged-current muon events in the energy range 1–11 GeV, with no oscillations, 100% efficiency for detection and charge identification of muons and extremely accurate energy measurement.

Figure A.5

Comparisons of energy distributions of μ ⁻ and μ ⁺ events in 500 kt-yr of ICAL, with no oscillations, 100% efficiency for detection and charge identification of muons and extremely accurate energy measurement.

The total number of muon events will be less with the fluxes at the INO site. As a result, the performance may be expected to be slightly worse than that calculated from the Kamioka fluxes. The extent to which the performance will be affected will depend on the quantity of interest, though. For example, the accuracy in the measurement of sin22𝜃 ₂₃ typically is controlled by the total number of events. As the total number of events with the fluxes at the INO site are about 14% smaller, we expect to take about 14% additional exposure to obtain the same level of accuracy as described in this review. On the other hand, as has been pointed out in figure 5 of [14], the hierarchy sensitivity comes mainly from the events with muon energy greater than 4 GeV. As figure A.5 indicates, the numbers of such muon events calculated using the two fluxes are nearly the same within statistical uncertainties. The results for the mass hierarchy determination are thus expected to be unaffected.

Appendix B. Neutrino oscillation probabilities in matter

Neutrino oscillation probabilities in matter are obtained by solving the propagation equation, which may be written in the flavour basis as

$$ i\frac{\mathrm{d}|\nu_{\alpha}(x)\rangle}{\mathrm{d}x} = H{\mid}{\nu_{\alpha}(x)}\rangle, $$

(B.1)

where |ν _α (x)〉 = (ν _e (x),ν _μ (x),ν _τ (x))^T. H is the effective Hamiltonian, given as

$$\begin{array}{@{}rcl@{}} H &=& \frac{1}{2E} ~U ~ \text{diag}(0,{\Delta} m^{2}_{21}, {\Delta} m^{2}_{31})U^{\dagger}\\ && +\, \text{diag}(V(x),0,0). \end{array} $$

(B.2)

Here E is the energy of the neutrino and ${\Delta } m^{2}_{ij} = {m_{i}^{2}} - {m_{j}^{2}}$ is the mass-squared difference between the neutrino mass eigenstates. The PMNS mixing matrix U relates the neutrino flavour eigenstates and mass eigenstates. V (x) is the matter potential arising due to the charged-current interaction of ν _e with electrons and is given in terms of the electron density n _e by $V(x) = \sqrt {2} G_{F} n_{e}(x)$. For antineutrinos, U → U ^∗ and V →−V. In general, for an arbitrary density profile one needs to solve the above equation numerically to obtain the probabilities. However, simplified analytic expressions can be obtained by assuming constant matter density. In such cases, one can diagonalize the above Hamiltonian to obtain

$$ H = \frac{1}{2E}~ U^{\mathrm{m}}~ \text{diag}((m_{1}^{\mathrm{m}})^{2}, (m_{2}^{\mathrm{m}})^{2}, (m_{3}^{\mathrm{m}})^{2})U^{\mathrm{m}\dagger}, $$

(B.3)

where $m_{i}^{\mathrm {m}}$ and U ^m denote the mass eigenvalues and mixing matrix in matter respectively. For a neutrino travelling a distance L, the flavour conversion probability in matter of constant density has an analogous expression as in the case of vacuum (see eq. (1.2)), and can be expressed as

$$\begin{array}{@{}rcl@{}} P_{\alpha\beta} (L)\! &=&\! \delta_{\alpha\beta} - 4 \sum\limits_{j>i} \text{Re} (U_{\alpha i}^{\mathrm{m}} U_{\beta i}^{\mathrm{m}^{\star}} U_{\alpha j}^{\mathrm{m}^{\star}} U_{\beta j}^{\mathrm{m}}) \sin^{2} {\Delta}_{ij}^{\mathrm{m}} \\ &&+ \,2 \sum\limits_{j>i} \text{Im}(U_{\alpha i}^{\mathrm{m}} U_{\beta i}^{\mathrm{m}^{\star}} U_{\alpha j}^{\mathrm{m}^{\star}} U_{\beta j}^{\mathrm{m}}) \sin (2 {\Delta}_{ij}^{\mathrm{m}}),\\ \end{array} $$

(B.4)

where the quantity ${\Delta }_{ij}^{\mathrm {m}}$ in the presence of matter is defined as

$${\Delta}_{ij}^{\text{m}} = \frac{ 1.27 ~ ({\Delta} m_{ij}^{2})^{\mathrm{m}} (\text{eV}^{2}) ~L (\text{km}) } { ~E (\text{GeV}) }, $$

with $({\Delta } m_{ij}^{2})^{\mathrm {m}}$ = $({m_{i}^{\mathrm {m}}})^{2}$ −$({m_{j}^{\mathrm {m}}})^{2}$ the difference between the squares of the mass eigenvalues $m_{i}^{\mathrm {m}}$ and $m_{j}^{\mathrm {m}}$ in matter.

To obtain tractable expressions, further assumptions need to be made. Many approximate analytic expressions for probabilities exist in the literature. However, different assumptions that lead to different approximate forms have different regimes for validity. To understand the results presented in this review, the probability expressions obtained under the following two approximations are mostly relevant:

• the one mass scale dominance (OMSD) approximation which assumes ${\Delta } m^{2}_{21}=0$,

• the double expansion in terms of small parameters $\alpha = {\Delta } m^{2}_{21}/{\Delta } m^{2}_{31}$ and sin𝜃 ₁₃ [176,177].

The condition on the neutrino energy and baseline for the validity of both approximations can be expressed as ${\Delta } m^{2}_{21} L /E \ll 1$. This translates to L/E ≪ 10⁴ km /GeV for typical values of the solar mass-squared difference ${\Delta } m^{2}_{21}$, and hence to L ≤ 10⁴ km for neutrinos of energy ${\mathcal {O}}(\text {GeV})$. Thus, these approximations are valid for most of the energy and path length ranges considered here. The OMSD approximation is exact in 𝜃 ₁₃ and works better near the resonance region. Below we give the probabilities relevant for the study presented in this Review in both OMSD and double expansion approximations, and discuss in which L and E regimes these are appropriate. Note that we give the expressions only for neutrino propagation through a constant matter density. This approximation is not applicable for neutrinos passing through the Earth’s core. However, it is enough for an analytic understanding of our arguments. All our numerical calculations take the variation of Earth’s density into account through the preliminary reference Earth model (PREM).

2.1 Appendix B.1: One mass scale dominance approximation

In this approximation, the Hamiltonian in eq. (B.2) can be exactly diagonalized analytically. Below we give the expressions for the muon neutrino survival probability P(ν _μ → ν _μ ) ≡ P _{μ μ} and conversion probability of electron neutrinos to muon neutrinos P(ν _e → ν _μ ) ≡ P _{e μ} , which are relevant for the atmospheric neutrinos at ICAL because the detector is sensitive to muon flavour. In the OMSD approximation, these can be expressed as

$$\begin{array}{@{}rcl@{}} P_{\mu \mu} &=& 1 - \cos^{2} \theta_{13}^{\mathrm{m}} \ \sin^{2} 2 \theta_{23} \\ &&\times\sin^{2}\big[1.27 ({\Delta} m^{2}_{31} + A + ({\Delta} m_{31}^{2})^{\mathrm{m}}) L/2E\big] \\ &&- \sin^{2} \theta_{13}^{\mathrm{m}} \ \sin^{2} 2 \theta_{23} \ \\ &&\times\sin^{2}\big[1.27 ({\Delta} m^{2}_{31} + A - ({\Delta} m_{31}^{2})^{\mathrm{m}}) L/2E\big] \\ &&- \sin^{4} \theta_{23} \ \sin^{2} 2\theta_{13}^{\mathrm{m}} \ \sin^{2} \big[1.27 {\Delta} m^{2}_{31} L/ E\big],\\ \end{array} $$

(B.5)

and

$$\begin{array}{@{}rcl@{}} P_{e \mu} = \sin^2 \theta_{23} \ \sin^2 2 \theta^{\mathrm{m}}_{13} \ \sin^2 \big[1.27 ({\Delta} m_{31}^2)^{\mathrm{m}} L/E\big].\\ \end{array} $$

(B.6)

In the above expressions, $({\Delta } m_{31}^{2})^{\mathrm {m}}$ and $\sin 2 \theta _{13}^{\mathrm {m}}$, the mass-squared difference and mixing angle in matter, respectively, are given by

$$\begin{array}{@{}rcl@{}}({\Delta} m_{31}^2)^{\mathrm{m}}\! &=&\! \sqrt{({\Delta} m^2_{31} \cos 2 \theta_{13}\! -\! A)^{2}\! +\! ({\Delta} m^{2}_{31} \sin 2 \theta_{13})^{2}},\\ &&\sin 2 \theta_{13}^{\mathrm{m}} = \sin 2 \theta_{13} \frac{\Delta m^2_{31}}{({\Delta} m_{31}^2)^{\mathrm{m}}}, \end{array} $$

(B.7)

where

$$\begin{array}{@{}rcl@{}} A (\text{eV}^{2}) &=& 2EV= 2\sqrt{2} G_{F} n_{e} E\\ &=& 0.76 \times 10^{-4} ~\rho(\text{g}/\text{cc}) ~E (\text{GeV}). \end{array} $$

When $A = {\Delta } m^{2}_{31} \cos 2 \theta _{13}$, we see a resonance. The resonance energy is given by

$$\begin{array}{@{}rcl@{}} E_{\text{res}} = \frac{\Delta m^2_{31} \cos 2\theta_{13}} {2 \sqrt{2} G_F n_e}. \end{array} $$

(B.8)

In table B.1, we give the average resonance energies for neutrinos travelling a given distance L through the Earth, for baselines ranging from 1000 to 10,000 km.

Table B.1 Values of E _res at various baselines using the line-averaged PREM [178] density ρ _avg. We have used $\frac {}{}{\Delta } m^{2}_{31}=2.5 \times 10^{-3} \text {eV}^{2}$ and $\frac {}{}\sin ^{2} 2\theta _{13}=0.1$.

It is seen from the table that the resonance energy is in the range 6–10 GeV for path lengths in the range 1000–10,000 km. These ranges are relevant for atmospheric neutrinos passing through Earth and hence provide an excellent avenue to probe resonant Earth matter effects. The importance of this can be understood by noting that the resonance condition depends on the sign of ${\Delta } m^{2}_{31}$. For ${\Delta } m^{2}_{31} \!\!>\!\! 0$ there is a matter enhancement in $\theta _{13}^{\mathrm {m}}$ for neutrinos, and a matter suppression in $\theta _{13}^{\mathrm {m}}$ for antineutrinos (as A → − A). The situation is reversed for ${\Delta } m^{2}_{31} \!\!<\!\! 0$. Thus, matter effects can differentiate between the two hierarchies and detectors with charge sensitivity (like ICAL) are very suitable for probing this.

For atmospheric neutrinos in ICAL, the most relevant probability is P _{μ μ} . The significance of the P _{e μ} channel is less than that of P _{μ μ} for two reasons: the number of electron neutrinos produced in the atmosphere is smaller, and more importantly, the probability of their conversion to muon neutrinos is also usually smaller than P _{μ μ} , so that their contribution to the total number of events in ICAL is small. It is not completely negligible though, because the value of 𝜃 ₁₃ is moderately large.

Note that P _{e μ} does not attain its maximum value at E = E _res even though $\sin 2 \theta _{13}^{\mathrm {m}}$ achieves its maximum value of unity at this energy, because the mass-squared difference $({\Delta } m_{31}^{2})^{\mathrm {m}}$ hits a minimum [179]. The values of $({\Delta } m_{31}^{2})^{\mathrm {m}} \sin 2 \theta _{13}^{\mathrm {m}}$ and P _{μ e} remain small for path lengths of L ≲ 1000 km. If L is chosen suitably large so as to satisfy $(1.27 {\Delta } m^{2}_{31} \sin 2 \theta _{13} L/E) \geq \pi /4$, then P _{e μ} can reach values ≥0.25 for sin22𝜃 ₂₃ = 1. One needs $L \gtrsim 6000$ km to satisfy the above condition. For such baselines and in the energy range 6–8 GeV, the resonant Earth matter effects lead to P _{e μ} in matter being significantly greater than its vacuum value [180].

Figure B.1

The oscillograms for the muon neutrino (left panel) and antineutrino (right panel) survival probabilities during their passage through Earth in E–$\cos \theta _{z}$ plane. The oscillation parameters used are 𝜃 ₂₃ = 45^∘, δ _CP = 0, ${\Delta } m^{2}_{31} = +2.45 \times 10^{-3}~\text {eV}^{2}$ (NH) and $\sin ^{2} 2\theta _{13} = 0.1$.

The muon neutrino survival probability is a more complicated function and can show both fall and rise above the vacuum value for longer baselines (∼10,000 km). Thus, the energy and angular smearing effects are more for P _{μ μ} . The maximum hierarchy sensitivity is achieved in this channel when resonance occurs close to a vacuum peak or dip, thus maximizing the matter effects because when there is resonant matter effect for one hierarchy, the probability for the other hierarchy closely follows the vacuum value. Figure B.1 shows oscillograms for muon neutrino and antineutrino survival probabilities in the case of normal hierarchy, in the plane of neutrino energy and cosine of the zenith angle 𝜃 _z . The plot in the left panel shows the resonant effect in the muon neutrino probabilities in the region between cos𝜃 _z between −0.6 to −0.8 and energy in the range 6–8 GeV. This feature is not present in the right panel because for the normal hierarchy, the muon antineutrinos do not encounter any resonance effect. The plot also shows the enhanced oscillation features due to the effects of the Earth’s core (cos𝜃 _z between −0.8 and −1.0) for neutrinos. For the inverted hierarchy, the muon antineutrino survival probability will show resonance effects, whereas the neutrino probabilities will not. The ICAL detector being charge sensitive can differentiate between neutrino and antineutrino effects and hence between the two hierarchies.

Equations (B.5) and (B.6) can also help us understand the octant sensitivity of atmospheric neutrinos arising due to resonant matter effects. The leading-order term in P _{e μ} in vacuum depends on sin2𝜃 ₂₃ sin22𝜃 ₁₃. Although this term is sensitive to the octant of 𝜃 ₂₃, the uncertainty in the value of 𝜃 ₁₃ may give rise to octant degeneracies. In matter, the $\sin ^{2} 2\theta _{13}^{\mathrm {m}}$ term gets amplified near resonance, and the combination $\sin ^{2} \theta _{23} \sin ^{2} 2 \theta _{13}^{\mathrm {m}}$ breaks the degeneracy of the octant with 𝜃 ₁₃. Also, the strong octant-sensitive nature of the term $\sin ^{4} \theta _{23} \sin ^{2} 2 \theta _{13}^{\mathrm {m}}$ near resonance can overcome the degeneracy due to the sin22𝜃 ₂₃-dependent terms. Unfortunately, the muon events in ICAL get contribution from both P _{μ μ} channel and P _{e μ} channel, and the matter effect in these two channels act in opposite directions for most of the baselines. This causes a worsening in the octant sensitivity of muon events at atmospheric neutrino experiments.

The OMSD probabilities are in the limit ${\Delta } m^{2}_{21} =0$ and have no dependence on the CP phase. These expressions match well with the numerical probabilities obtained by solving the propagation equation in the resonance region, i.e. for the baseline range 6000–10,000 km. Accelerator-based experiments like T2K and NO ν A have shorter baselines and lower matter effects, and lie far from resonance. For these experiments, the dominant terms in P _{μ μ} are insensitive to the hierarchy and octant. Consequently, the relative change in probability due to the hierarchy /octant-sensitive subdominant terms is small. Therefore, the P _{μ μ} oscillation channel does not contribute much to the hierarchy and octant sensitivity of T2K and NO ν A. However, these experiments get their sensitivity primarily from the ν _μ → ν _e conversion probability. For this case, the double expansion up to second order in α and sin𝜃 ₁₃ works better.

2.2 Appendix B.2: Double expansion in α and sin𝜃 ₁₃

In accelerator experiments, high-energy pions or kaons decay to give muons and muon neutrinos /antineutrinos. One can study the muon neutrino conversion probability P(ν _μ → ν _e ) ≡ P _{μ e} in these with a detector sensitive to electron flavour. In order to study the effect due to ${\Delta } m^{2}_{21}$ and the CP phase δ _CP it is convenient to write down the probabilities as an expansion in terms of the two small parameters, $\alpha = {\Delta } m^{2}_{21}/{\Delta } m^{2}_{31}$ and sin𝜃 ₁₃, to second order (i.e. terms up to α ²,sin2𝜃 ₁₃ and α sin𝜃 ₁₃ are kept) [176,177].

$$\begin{array}{@{}rcl@{}} P_{\mu e }&= & \sin^{2} 2\theta_{13} \sin^{2}\theta_{23} \frac{\sin^{2}{[(1 -\hat A){\Delta}]}}{(1-\hat A)^{2}} \\ && +\, \alpha \sin{2\theta_{13}} \sin{2\theta_{12}} \sin{2\theta_{23}} \cos{\!({\Delta} + \delta_{\text{CP}})}\\ &&\times\,\frac{\sin(\hat A {\Delta})}{\hat A} \frac{\sin{[(1-\hat A){\Delta}]}}{(1-\hat A)} \\ && +\, \alpha^{2} \sin^{2}\theta_{12} \cos^{2}\theta_{23} \frac{\sin^{2} (\hat A {\Delta})}{\hat A^{2}}\\ &&+\, {\mathcal{O}}(\alpha^{3},\alpha^{2} s_{13},\alpha s_{13}^{2},s_{13}^{3}). \end{array} $$

(B.9)

The notations used in writing the probability expressions are: ${\Delta } \equiv {\Delta } m^{2}_{31}L/4E$, s _{i j} (c _{i j} ) ≡ sin𝜃 _{i j} (cos𝜃 _{i j} ), $\hat {A} = 2\sqrt {2} G_{F} n_{e} E / {\Delta } m^{2}_{31}$. For neutrinos, the signs of $\hat {A}$ and Δ are positive (negative) for NH (IH). The sign of $\hat {A}$ as well as δ _CP reverse for antineutrinos. This probability is sensitive to all the three current unknowns in neutrino physics – hierarchy, octant of 𝜃 ₂₃ as well as δ _CP – and is often hailed as the golden channel. However, the dependences are interrelated and extraction of each of these unknowns depends on the knowledge of the others. Specially, the complete lack of knowledge of δ _CP gives rise to the hierarchy- δ _CP degeneracy as well as the octant- δ _CP degeneracy in these experiments, through the second term in eq. (B.9).

The above expressions reduce to the vacuum expressions for shorter baselines for which A → 0. For such cases there is no hierarchy sensitivity. The hierarchy sensitivity increases with increasing baseline and is maximum in the resonance region. The resonance energy at shorter baselines is > 10 GeV and therefore these experiments cannot probe resonant Earth matter effects.

2.3 Appendix B.3: Probability for reactor neutrinos

A crucial input in the analysis presented in this Review is the value of 𝜃 ₁₃, measured by the reactor neutrino experiments. The probability relevant for reactor neutrinos is the survival probability for electron antineutrinos $P(\bar {\nu }_{e} \to \bar {\nu }_{e}) \equiv P_{\bar {e} \overline {e}}$. Since reactor neutrinos have very low energy (order of MeV) and they travel very short distances (order of km), they experience negligible matter effects. The exact formula for the survival probability in vacuum is given by [180,181]

$$ \begin{array}{cccccc} P_{\overline{e} \overline{e}} & = & 1 - c_{13}^4 \sin^2 2\theta_{12} \sin^2 {\Delta}_{21} - \sin^2 2\theta_{13} \sin^2 {\Delta}_{31} \\ && +\,\sin^2 2\theta_{13} s_{12}^2[\sin^2 {\Delta}_{31} - \sin^2 {\Delta}_{32}]. \end{array} $$

(B.10)

As this is independent of δ _CP and matter effects, the probability is the same for neutrinos and antineutrinos (assuming CPT conservation, i.e. the same mass and mixing parameters describing neutrinos and antineutrinos).

Appendix C. The Vavilov distribution function

The Vavilov probability distribution function is found to be suitable to represent the hit distributions of hadrons of a given energy in the ICAL, as has been observed from figure 4.8. The Vavilov probability density function in the standard form is defined by [182]

$$ P(x; \kappa, \beta^{2}) = \frac{1}{2 \pi i}{\int}_{c-i\infty}^{c+i\infty} \phi(s) \mathrm{e}^{x s} \mathrm{d}s, $$

(C.1)

where

$$ \phi(s) = \mathrm{e}^{C} \mathrm{e}^{\psi(s)},~~~C = \kappa (1+\beta^{2} \gamma) $$

(C.2)

and

$$\begin{array}{@{}rcl@{}} \psi(s) &=& s \ln \kappa + (s+\beta^{2} \kappa) \cdot \left[{\int}_{0}^{1} \frac{1 - \mathrm{e}^{-st/\kappa}}{t} \,\mathrm{d} t - \gamma\right]\!\!\\ &&-\, {\kappa}\mathrm{e}^{-s/\kappa}, \end{array} $$

(C.3)

where γ = 0.577… is the Euler’s constant.

The parameters mean and variance (σ ²) of the distribution in eq. (C.1) are given by

$$ {\text{Mean}}= \gamma -1 - \ln\kappa-\beta^{2};~~~ \sigma^{2} = \frac{2-\beta^{2}}{2\kappa}. $$

(C.4)

For κ ≤ 0.05, the Vavilov distribution may be approximated by the Landau distribution, while for κ ≥ 10, it may be approximated by the Gaussian approximation, with the corresponding mean and variance.

We have used the Vavilov distribution function P(x;κ,β ²) defined above, which is also built into ROOT, as the basic distribution for the fit. However, the hadron hit distribution itself is fitted to the modified distribution (P ₄/P ₃)P((x − P ₂)/P ₃; P ₀, P ₁), to account for the x-scaling (P ₃), normalization P ₄ and the shift of the peak to a non-zero value, P ₂. Clearly, P ₀ = κ and P ₁ = β ². The modified mean and variance are then

$$\begin{array}{@{}rcl@{}} \text{Mean}_{\text{Vavilov}} &=& (\gamma-1-\ln {P}_{0} - {P}_{1}) {P}_{3}+ {P}_{2}, \\ &&\sigma^{2}_{\text{Vavilov}}= \frac{(2- {P}_{1})}{2{P}_{0}} {P}_{3}^{2}. \end{array} $$

(C.5)

These are the quantities used while presenting the energy response of hadrons in the ICAL detector.

Appendix D. Hadron energy resolution as a function of plate thickness

A potentially crucial factor in the determination of hadron energy and direction is the thickness of absorber material, namely iron plate thickness in ICAL. In all the simulation studies reported here, we have assumed that the thickness of the iron plate is 5.6 cm, which is the default value. While not much variation in this thickness is possible due to constraints imposed by total mass, physical size, location of the support structure, and other parameters like the cost factor etc., we look at possible variation of this thickness in view of optimizing the hadron energy resolution [10]. The hadron energy resolution is a crucial limiting factor in reconstructing the neutrino energy in atmospheric neutrino interactions in the ICAL detector. This information is also helpful because ICAL is modular in form and future modules may come in for further improvements using such analyses.

Naively, this can be achieved by simply changing the angle of propagation of the particle in the simulation, because the effective thickness is (t/cos𝜃). In the case of muons this itself may be sufficient to study the effect of plate thickness. However, in an actual detector, the detector geometry – including support structure, orientation as well as the arrangement of detector elements – imposes additional non-trivial dependence on thickness. Therefore, we study hadron energy resolution with the present arrangement of ICAL by varying the plate thickness, while other parameters are fixed. The analysis was done by propagating pions in the simulated ICAL detector at various fixed energies, averaged over all directions in each case.

The hit distribution patterns for 5 GeV pions propagated through sample plate thicknesses in the central region are shown in figure D.1. The methodology is already discussed in §4 and we shall not repeat it here. For comparing the resolutions with different thicknesses we use the mean and rms width (σ) of the hit distributions as functions of energy.

Figure D.1

Hit distribution of 5 GeV pions propagated through sample iron plate thicknesses [10].

The hadron energy resolution is parametrized as

$$ \left( \frac{\sigma}{E}\right)^{2} =\frac{a^{2}}{E}+b^{2}, $$

(D.1)

where a is the stochastic coefficient and b is a constant, both of which depend on the thickness. We divide the relevant energy range 2–15 GeV into two subranges, below 5 GeV and above 5 GeV. Below 5 GeV, the quasielastic, resonance and deep inelastic processes contribute to the production of hadrons in neutrino interactions in comparable proportions, while above 5 GeV the hadron production is dominated by the deep inelastic scattering. The results for the energy resolution as a function of plate thickness are shown in figure D.2. Note that here, we show the square of the resolution instead of the resolution itself.

Figure D.2

Plots of (σ/mean)² as a function of 1/E. The data as well as fits to eq. (D.1) are shown in the energy range 2–4.75 GeV (left) and 5–15 GeV (right). The thickness is varied from 2.5 to 8 cm [10].

The stochastic coefficient a as a function of thickness is obtained from a fit to the hadron energy resolution, and is shown in figure D.3 as a function of plate thickness for the two energy ranges as in figure D.2.

Figure D.3

The stochastic coefficient a obtained from the fit to eq. (D.1) is shown in the two energy ranges as a function of plate thickness [10].

The analysis in the two energy ranges shows that the thickness dependence is stronger than $\sqrt {t}$ which is observed in hadron calorimeters at high energies (tens of GeV) [183]. In fact, at the energies of relevance to us, the thickness dependence is not uniform but dependent on the energy. This is borne out by two independent analyses: in the first one we obtain the thickness dependence of the stochastic coefficient a and in the second analysis we directly parametrize the energy resolution as a function of thickness at each energy. Typically, instead of t ^0.5, we find the power varying from about 0.65 to 0.98 depending on the energy.

Finally, we compare the ICAL simulations with varying thicknesses with MONOLITH and MINOS and their test beam runs. This is a useful comparison because test beam runs with ICAL prototype have not been done till now. The data from the above detectors can however be used for the validation of the ICAL simulations results.

The test beam results for the Baby MONOLITH (BM) detector at CERN with 5 cm thick iron plates [184,185] have been obtained when the beam energy is in the range 2–10 GeV. In order to provide a comparison, we have simulated the ICAL detector response with 5 cm iron plates, for single pions of 2–10 GeV energy incident normally on the detector at a fixed vertex. Also, in order to be consistent with the BM parametrization, the energy resolution σ _E /E is fitted to the function $A/\sqrt {E} + B$. A comparison of the ICAL- simulated results with the BM beam results, along with the respective fits, is shown in figure D.4.

Figure D.4

The energy response of ICAL detector with 5 cm thick iron plates with single pions in the energy range 2–10 GeV, propagated from a fixed vertex in the vertical direction [10] compared with the data from MONOLITH test beam run [184,185].

For BM, an energy resolution of $\sigma _{E}/E = 68\%/ \sqrt {E} + 2\%$ was reported [184]. However, no errors on the parameters A and B were specified. Our fit to the same BM data gives A _BM = (66 ± 5)% and B _BM = (1 ± 2)%, which also gives an estimation of errors on these parameters. The fit for the ICAL resolution gives the parameter values A _ICAL = (64 ± 2)% and B _ICAL = (2 ± 1)%. The consistency of our simulated results with the beam results of BM testifies to the correctness of our approach.

In the test beam run of MINOS with aluminium proportional tube (APT) active detectors and 1.5 inch (4 cm) steel plates, a hadron energy resolution of $71\%/\sqrt {E}\pm 6\%$ was reported in the range 2.5–30 GeV [186]. ICAL simulation with 4 cm iron plates in the same energy range gave $61\%/\sqrt {E}\pm 14\%$. The results are compatible within errors, because the two detector geometries are very different.

Obviously the final choice of the plate thickness depends not only on the behaviour of hadrons but also on the energy range of interest to the physics goals of the experiment. There are also issues of cost, sensitivity to muons and even possibly electrons. The thickness dependence study summarized here provides one such input to the final design.