The proposed Iron Calorimeter (ICAL) detector at Indian Neutrino Observatory (INO) will be a large (50 kt) magnetized detector located 1270 m underground at Bodi West Hills in Tamilnadu. ICAL is capable of identifying the charge of the particles. In this paper its potential for the measurement of the muon charge ratio is explored by means of a detailed simulation-based study, first using the CORSIKA code and then comparing it with an analytical model (the pika model). The simulated muon charge ratio is in agreement with the existing experimental observations; its measure can be extended by INO-ICAL up to 10.50 TeV and up to 60 degrees.

1. Introduction

The India-based Neutrino Observatory (INO) is planned to be built at Theni district of Tamilnadu in southern India. This neutrino observatory will have a 52 kton magnetized Iron CALorimeter (ICAL) detector [1], which is designed to study the physics of atmospheric neutrino flavor oscillations. The main goal of the ICAL detector is to make precise measurements of neutrino oscillation parameters and to determine the neutrino mass hierarchy. Roughly 1.2 km underground, the INO-ICAL detector will be the biggest magnetized detector to measure the cosmic ray muon flux with a capability to distinguish between and .

Cosmic rays below energies of about eV are thought to be of Galactic origin, the most probable source being Supernova remnants, that constantly hit the outer surface of the Earth’s atmosphere. The magnetic field of the Earth tends to exclude low energy charged particles (E 1 GeV), but high energy particles easily manage to penetrate it. These cosmic rays are primarily composed of high energy protons, alpha particles and heavier nuclei, which interact with the top layer of Earth’s atmospheric nuclei and produce pions and kaons. These particles further decay to muons

The charge ratio of these energetic muons is an important measurable which carries information about (i) /K hadronic production ratio, (ii) the composition of cosmic ray primaries, (iii) the contribution of charmed hadrons at very high energy, and (iv) the neutrino flux at very high energies, etc. The measurement of the charge ratio of muons (), up to tens of TeV, has been made by several experiments (L3+C[2], MINOS[3], CMS[4], OPERA[5], etc.) and they have found that charge ratio increases with an increase in energy. A charge ratio larger than unity reflects abundance of protons over heavier nuclei in the primary cosmic radiation. As the energy increases, the contribution of muons from kaon decay also increases because the longer-lived pions are more likely to interact before decaying than the shorter-lived kaons, competing processes could be the production of heavy flavors or a heavier primary composition. The muon charge ratio increases with energy since the production of positive kaons is favored by the associated production of cosmic ray protons (p ). The charm production is expected to contribute at very high energies but not relevant in the energy range explored by the experiment discussed in this work.

In this work we present the response of ICAL detector for high energy muon measurements using the latest version of the INO-ICAL code, which is based on the Geant4 toolkit[12]. For simulating the charge ratio at the underground detector using the surface muon data, we have added the real topography of the hill to the INO-ICAL code. The hill topography is defined using the standard rock chemical composition in the Geant4 code. After adding the hill topography in the code energy threshold for INO is estimated by propagating muons vertically from the various depths in the rock to the detector. After setting the threshold, muon charge ratio is calculated using cosmic muon flux, which is first generated using CORSIKA (COsmic Rays SImulation for KAscade)[13] and then this analysis is repeated with the analytical model known as “pika model"[14]. In CORSIKA, for generating the muon flux only vertical primary protons spectrum is considered for our analysis. We have generated all the secondary particles fluxes at the top of the hill and then selected only those muons, which have energy above the threshold energy for our analysis. The generated muon flux is used by the Geant4 code as an input flux for muon charge ratio analysis by propagating it through the hill to the detector. Muons with energy greater than the threshold energy will enter the detector and produce tracks. Track reconstruction is performed using the Kalman Filter Technique [15]. In the central region of each module, the typical value of the magnetic field strength is about 1.5 T in the y-direction, which is obtained using the MAGNET6.26 software [16]. Finally, we estimate the energy and the number of and particles which pass through rock and reach the detector. For this analysis the ICAL detector momentum and charge resolution are considered. Our estimated results for muon charge ratio using CORSIKA and “pika" model are compared with each other as well as with the existing experimental data.

This paper is organized as follows: In Section 2, we discuss the ICAL detector geometry, the hill topography, and the threshold energy of muons. In Section 3, we explain the response of high energy muons at INO-ICAL detector and have incorporated the detector efficiency to reconstruct the muon momentum and charge. Section 4 briefly explains the cosmic muon flux generation at the surface and the estimation of the muon charge ratio at the underground detector using CORSIKA. This section also explains the estimation of the charge ratio using the “pika" model. In Section 5 the systematic uncertainties addressed in this work are explained. In Sections 6 and 7, we discuss our results and conclusions, respectively.

2. THE ICAL Detector Geometry and the Bodi Hill Profile

2.1. Detector Geometry

The INO-ICAL detector will have a modular structure. The full detector consists of three modules, each of size 16 m 16 m 14.45 m. These modules are placed along the x-direction as shown in Figure 1, with a gap of 20 cm between them and the origin is taken to be the center of the central module, while the remaining horizontal and transverse direction are labelled as y-axis and z-axis, respectively. The z-axis points vertically upwards; hence the polar angle is equal to the zenith angle . Each module comprises 151 horizontal layers of iron plates of area 16.0 m 16.0 m and thickness 5.6 cm with a vertical gap of 4.0 cm, interleaved with Resistive Plate Chamber (RPCs) as an active detector elements.

The total mass of the detector is about 52 kton, excluding the weight of the coils. The basic RPC units of size 1.84 m 1.84 m 2.5 cm are placed in grid-format within the air gap [1]. More details of the detector elements and its parameters are discussed in [1]. This geometry has been simulated by the INO collaboration using GEANT4 [12] and its response for muon with energy of few tens of GeV is studied in [17].

2.2. Digital Data of the INO Peak at Pottipuram Site

We have generated the elevation profile of the mountain around INO cavern as shown in Figure 2, at Bodi West Hills in Theni district at Tamilnadu.

The latitude and longitude at the cavern location are as follows:(i)Latitude: (ii)Longitude:

The average height of the peak around the tunnel region is 1587.32 m from the sea level and the height of the detector location from the sea level is around 317 m; hence the rock cover above the detector is around 1270 m.

We are using the mountain profile provided by Geological Survey of India, as shown in Figure 2. For the detailed analysis of underground muon charge ratio at INO-ICAL detector, standard rock overburden of density 2.65 gm/ is used in our work. The rock type around the INO location is charnockite whose density is around 2.6 gm/. The real rock density is roughly equal to the standard rock density.

2.3. Muon Energy Loss Analysis through Standard Rock

The muon charge ratio measured in underground experiments is the ratio of to at the detector. Under the hypothesis of identical energy losses for positive and negative muons, the muon charge ratio on the surface would not be modified by propagation of muons in the rock down to the detector. At higher energies, the fluctuations in energy losses by and are important and an accurate calculation of muon charge ratio requires a simulation that accounts for stochastic loss processes [18]. Muons of energy at the Earth’s surface lose energy [6] as they traverse a slant depth () through the Standard rock to reach the detector. The muon energy loss is described by the well-know Bethe-Bloch formula: where parameterizes the energy loss by collisions and parameterizes the energy loss by three radiative processes, namely, Bremsstrahlung, Pair production, and Photo-nuclear interactions. The energy loss parameters and for the standard rock, as a function of energy, are taken from reference [6]. The muon energy reconstructed at the detector can be related to the surface muon energy by [6] The values of the energy loss parameters and depend on the average composition of the rock. Energy loss processes of and are almost the same but there are corrections of the order of the fine structure constant , which were shown to depend almost only on the ionisation energy loss [14]. The change in charge asymmetry can be written as

where is the difference in energy loss for and . As shown in [14], a change from standard rock composition to a different rock composition produces a negligible effect on the muon charge ratio. As such, our use of standard rock in place of charnockite is justified. A good knowledge of the rock overburden is required to estimate the muon momentum at the underground detector from the surface muon momentum. The variation in the composition of rock requires the muon energy loss estimation in all differential segments of rock, where the density is constant, and the addition of these losses gives the complete muon energy loss. After feeding the mountain geometry within the ICAL code, 10000 muons were propagated vertically to the detector through different depths of the rock shielding. For energy loss analysis in the rock overburden the different depths (from the origin of the central module) considered in our analysis are 1, 2, 3, 3.37, 4, 5, and 6 km.w.e.; the value 3.37 km.w.e. refers to the highest peak value or maximum rock overburden from the origin of the detector. The muon energy thresholds for increasing slant depths are shown in Figure 3; a threshold of 1.6 TeV corresponds to the experimental depth of 3.37km.w.e.

3. High Energy Muon Response at ICAL

In this work we have done simulation for low energy (1-20GeV) muons as well as for the high energy (20-500GeV) muons using the latest version of INO-ICAL code. Details of the detector simulation for low energy muons with older version of INO-ICAL code are already published in the reference [17]. In this work we have followed the same approach as in reference [17], for muon response analysis up-to energies 500 GeV inside the detector. Detector response for the peripheral region of the detector was performed in [19], while in this work detector response for the central region is performed. For simulating the response of high energy muon in the ICAL detector, 10000 muons were propagated uniformly from a vertex, randomly located inside the central volume of the detector of dimensions 8 m 8 m 10 m. In our simulation we have considered a uniform magnetic field of strength 1.5T.

3.1. Tracks and Event Selection

The muon track reconstruction is based on a Kalman filter algorithm that takes into account the local magnetic field. This algorithm is briefly described in [17]. It is used to fit the tracks based on the bending of the tracks in the magnetic field, by collecting the hitting points in multiple layers of iron plates. Only tracks for which the quality of fit is better than 10 are used in this analysis. Each muon track is further analysed in order to identify its direction, charge, and momentum. Both fully contained and partially contained muons events are included in this work.

3.2. Momentum Reconstruction Efficiency of ICAL

The momentum reconstruction efficiency () is defined as the ratio of the number of reconstructed events, , to the total number of generated events, . We have with uncertainty, Figure 4 shows the muon momentum reconstruction efficiency as a function of input momentum for different values. The left and right panels demonstrate detector response for low and high energy muon momentum respectively.

One can see that the momentum reconstruction efficiency depends on the incident particle momentum and the angle of propagation. The rise in efficiency of muons, as shown in the left panel of Figure 4, is due to the fact that higher energy muons can cross more number layers. A muon with fixed energy hitting the detector at larger zenith angle will give lesser no. of hits as compared to the same muon hitting the detector vertically (angle dependence). At higher energies the reconstruction efficiency of muons, as shown in the right panel of Figure 4, initially drops and then becomes almost constant (around 80%). The reason for the drop in momentum reconstruction efficiency of muons at high energies (above about 50 GeV) is that the muons travel nearly straight without being deflected in the magnetic field of the detector.

Charge ratio of muons typically depends on the reconstructed muon momentum at the detector; hence a change in momentum value will lead to an uncertainty in the charge ratio. Figure 5 indicates a shift in the mean reconstructed muon momentum distribution at fixed zenith angle and this shift increases with increasing input momentum. Figure 6 shows the shift in the muon momentum (-) as function of input muon momentum for various zenith angles. This shift can arise due to the non-linearity, etc., at the larger momentum, the effect is simply more prominent, and a detailed discussion can be found in the reference [20].

3.3. Relative Charge Identification Efficiency of ICAL

The charge of the particle is determined from the direction of curvature of the particle track in the magnetic field. Thus it is crucial for the determination of the neutrino mass hierarchy, which is the main aim for the INO project, and for the estimation of the atmospheric muon charge ratio. Relative charge identification efficiency is defined as the ratio of the number of events with correct charge identification, to the total number of reconstructed events, of same charge, i.e., with uncertainty,

Figure 7 shows the relative charge identification efficiency as a function of input momentum for different values. Here the left and right panels demonstrate the detector response for low and high energy muon momentum, respectively. Muons propagated with small momentum will cross lesser number of layers and will give lesser number of hits. This may lead to an incorrect reconstruction of the direction of bending and wrong charge identification. Hence the charge identification efficiency is relatively poor at lower energies, as can be seen from the left plot of Figure 7. As the energy increases, the length of the track also increases due to which the charge identification efficiency improves. Beyond a few GeV/c, the charge identification efficiency becomes roughly constant at 98-99%. But as the energy increases further above 50 GeV the charge identification efficiency starts decreasing because in the constant magnetic field the deflection of highly energetic muon will be less, and at a certain energy the muon will leave the detector without getting deflected by magnetic field. Due to this limitation the charge identification efficiency falls to 70% at muon momentum 400 GeV/c. In our work we have simulated the muon charge ratio using CORSIKA up to muon momenta of 400 GeV/c, where the detector charge identification efficiency is 70% for the central region. The 70% efficiency of the detector is only for the muons generated by vertical protons using CORSIKA, whereas the same studies have been extended in Section 4.2 for 90% efficiency.

4. The Muon Charge Ratio Evaluation

The muons charge ratio is estimated using two approaches. In the first approach the estimation is performed using the muon flux generated using CORSIKA package, while in the second approach the muon flux is estimated using analytical expression. We discuss both the cases in the later two sections.

4.1. Muon Charge Ratio Analysis Using CORSIKA Data

CORSIKA package is used here for the air shower simulation. For simulating the air shower, we exploited QGSJET II-04 [21] for high-energy interactions, i.e., for the energy range 100 GeV up to 100 EeV, and GHEISHA[22] for low energy interactions, i.e., energy up to few 100 GeV. Preliminary tests with a different low energy interaction code showed negligible difference in the final determination of the muon charge ratio. We have generated the muon flux at an observational level of 2.0 Km. above the sea level, which is the top surface of the hill at the INO location as shown in Figure 2. This muon flux is used as an input for GEANT4 based INO-ICAL code. Muon energy loss analysis in the rock suggests that only high energy muons will be able to reach the detector, so for the cosmic muon analysis only high energy muons are generated for this work. For the generation of primary proton flux above the atmosphere, a power law parameterization is considered (). Proton flux is generated at zenith angle in the energy range 1-400 TeV. The entire energy range is divided into 100 equal bins with the bin size of 200 GeV and each bin has 1000000 protons. The generated muon spectrum from the selected primary spectrum is shown in Figure 8 as a function of muon momentum using a log scale. The expected time exposure for collecting this data is around 3 years [23]. To validate the surface muon flux data generated by us using CORSIKA, we have compared our flux with the muon flux evaluated by T. Gaisser [11] and the available experimental data. The flux presented in Figure 9 with the solid blue line consists of several primaries species, with the corresponding spectra inferred from measurements still comparable with our data because of the abundance of primary protons in the primary spectra. This is illustrated in Figure 9 as a function of muon momentum (only vertical muon flux is used). Figure 9 is reproduced by taking the data form the existing experimental fluxes ( [2],[7], +[8], [9], [10]) and the flux (solid blue line) from [11]. The generated muon flux at the top of the hill surface will spread in the large area of X-Y plane as shown in Figure 10 where the observational plane is at 2.0 km. above the sea level. Most of the muons will be distributed in the range -100 m X 100 m and -100 m Y 100 m of X-Y plane and the origin of the plane will be the shower axis.

Density of the muons near the shower axis is very high as shown in Figure 10 and the uncertainty in the distribution of the flux is considered in the statistical error bars of the final results.

In this section we determine the muon charge ratio, , as a function of the reconstructed muon momentum at the INO-ICAL detector.

From the generated muon data sample only energetic muons with energies in the range 1600GeV to 2000GeV are selected. This muon flux sample is then divided into 8 equal data samples and each has a bin size of 50 GeV. Each data sample is propagated through the rock to the detector. The reconstructed muon momentum distribution for one of the sets (1600-1650GeV) is shown in Figure 11. This figure shows a clear separation of and flux as a function of reconstructed muon momentum at the detector, corresponding energy of the muons at the earth surface is 1600-1650GeV.

The vertical muon charge ratio has been evaluated up to 2 TeV and plotted in Figure 12 and the corresponding values are listed in the second columns of Table 1, whereas the third column shows the vertical muons charge ratio at the surface. The error bar shown in the Figure 12 represents statistical errors only.

4.2. Charge Ratio Analysis Using “Pika" Model

Since the muon charge ratio also depends on the zenith angle; hence we study it as a function of . In this section, the muon charge ratio at INO-ICAL detector is estimated using atmospheric muon energy spectrum [14]. A formula for atmospheric muon energy spectrum given by Gaisser [11] is

This formula is valid when muon decay is negligible GeV) and the curvature of the Earth is neglected . The two terms inside the curly braces give the contributions to muon flux from charged pion and kaons. The rise in the muon charge ratio can be understood from the properties of and K mesons, and the observation of the rise in muon charge ratio can be used to determine the and ratio. A method to study separately the positive and negative muons intensities has been developed in reference [14] and it is known as “pika" model. This model provides the muon energy spectrum for +ve and -ve muons by using the positive fraction parameters and ,

In the above equation and have been replaced by their numerical values. The parameter specifies the relative role of kaon decays, compared to the ones of pions. It depends upon the ratio, branching ratios, and kinematic factors, which arise due to difference between the and K masses. The numerical value of this parameter ( = 0.054) is discussed in [24]. Equations (10) and (11) give the expression for the surface muon charge ratio:

The charge ratio of muon from pion decay is = and from kaon decay is . Equation (12) is a function of only; all other parameters are extrapolated from the available experimental data; the parameters value used in our work are = 0.55100.0006 and = 0.70060.0061; a detailed discussion can be found in [14]. The pika model describes the data well, from the experiment L3+C [2], MINOS Near and Far Detector[3], in comparison to the calculations made by Honda [25], CORT [26], and Lipari [27]. The charge ratio estimated at the INO detector is shown in Figure 13, for charge ratio calculation with the pika model muons up to momentum of 500 GeV/c at the detector. The fourth column of Table 1 shows the charge ratio for vertical muons only. Muon flux generation using the Gaisser formula at INO and muon charge ratio estimation is performed in [28].

5. Systematic Uncertainty

The simulation technique used for the estimation of the muon charge ratio tends to cancel out many systematic uncertainties; however few uncertainties which do not cancel out remain in our analysis. The systematic uncertainties considered in our work are (I) those arising due to the range estimation in the standard rock for and as discussed in Section 2.3: (II) those arising due to the shift in the mean of the reconstructed muon momentum as discussed in Section 3.3: (III) those arising due to the track reconstruction algorithms in magnetic field, which comes mainly due to the non-uniformity in the magnetic field and available dead spaces as discussed in detail in reference [17]: (IV) those arising due to variation in the surface height, (as stated, we did our analysis for a fixed height (1300m) and then varied this height by 50 m to get the error in charge ratio results due to this variation). The total average systematic uncertainties as discussed above are 0.055, which is included in our error bars as shown in Figure 13.

To evaluate the systematic uncertainties listed above, we exploited a Monte Carlo code for the estimation of the slant depth, in unit of kilometer water equivalent (where 1 km.w.e = g/). For the estimation of the accurate slant depth, we have zoomed 3.0 km of the surrounding area about the origin of the detector location. We observe that the height of the rock cover in this region is around 1300 ± 50 m above the detector. Equation (3) shows a relationship between slant depth () and threshold energy of muon to reach the detector after travelling a distance . The values of parameters and in equation (3) depend on muon energy and the values of these parameters at different muon energies are taken from [6]. Using these energy dependent values of and we have estimated the threshold (for muon) for the corresponding slant depth covering in the zenith angle range . In our simulation we consider only downgoing muons and assume constant height (1300 ± 50 m). The azimuthal angle is varied to cover all the 5 planes, which include the top and four sides of the detector, excluding the bottom surface. This simulation is limited to downgoing muons; therefore the bottom surface of the detector is excluded. Detector efficiency is more than 80% in the energy range of 1-200 GeV for the zenith angle . In our work charge ratio estimation shown in Figures 13 and 14 corresponds to the zenith angle range . The surface near the region of cavern is almost constant in height (1300 m) but there is variation in the height along the higher latitude (9.97-9.98 deg) and this variation is around 50 m. We estimated the uncertainties in surface muon energy due to the variation in height and this error is incorporated in muon charge ratio results, as shown in Figure 13.

6. Discussion

This work describes the method used to simulate the muon flux and charge ratio at the underground INO-ICAL detector. Using the real topography of the rock overburden, we have estimated the threshold energy of the muon to reach the detector from the top surface of INO rock overburden. Our work in this paper is divided into three parts: in the first part muon charge ratio is estimated using CORSIKA flux, in the second part muon charge ratio is estimated using the “pika" model, and finally muon charge ratio results estimated by CORSIKA are compared with the results obtained by “pika" model and available experimental data. The analysis is performed only for the downgoing atmospheric muons. The muon charge ratio analysis at INO-ICAL using CORSIKA flux is limited to the flux entering from the top surface of the detector only. The same analysis is repeated using the “pika" model which incorporates five out of six surfaces of the detector. The muon charge ratio estimation is based on the strength of the detector’s magnetic field.

Due to its large dimensions, ICAL detector can detect muons coming from large zenith angles, as discussed in Sections 3.2 and 3.3. Momentum reconstruction and charge identification efficiency is more than 90% in 1-100 GeV energy range but falls with increasing energy. The ICAL charge identification efficiency drops from 80% to 70% for the input momentum range 200 GeV/c to 400 GeV/c. Our analysis done with CORSIKA flux is limited to momentum reconstruction and charge identification efficiency of the detector more than 70%, as shown in Figures 4 and 7. The analysis of detector efficiency is performed at various fixed energies and zenith angles while the azimuthal angle is uniformly averaged over the entire range -. A detailed simulation is peformed with real topography of the surface surrounding the INO location for the zenith angle range and azimuthal angle -.

Estimation of muon charge ratio at the underground INO-ICAL detector as well as at the top of the hill surface is performed using CORSIKA flux having surface muon energy in the range 1.60-2.00 TeV. The estimated muon charge ratio varies between 1.3 and 1.4 as shown in Figure 13. The average value of the charge ratio observed at underground detector is 1.3470.019 and corresponding average muon charge ratio at the surface is 1.358. Charge ratios observed in the same energy range by L3+C[2], MINOS[3], CMS[4], and OPERA[5] experiments are shown in Figure 14.

The muon charge ratio analysis at higher energies and higher zenith angles is included in our work by using “pika" model. For the estimation of the muon charge ratio, the muon momentum reconstruction efficiency and charge identification efficiency of INO-ICAL detector are taken into consideration. And charge ratio for surface muon energy in the range 1.6-10.50 TeV surface muon energy is estimated. The estimated muon charge ratio as a function of is shown in Figure 13. The estimated average muon charge ratio is found to be 1.437 ± 0.055, which is slightly higher than the muon charge ratio results estimated by available experimental data.

7. Conclusion

In this work we are presenting the cosmic rays physics potential of INO-ICAL and demonstrate that INO-ICAL will have a wide surface (48 16 ) area providing large acceptance with respect to the other experiments mainly OPERA, having similar depth 3.8 km water equivalent, and thus will reduce the statistical Uncertainty. The estimated results using CORSIKA can be used for relating the secondary spectrum corresponding to the primary cosmic rays spectrum, which will be helpful for studying the primary cosmic rays spectrum. Charge ratio estimation is further extended for larger zenith angle by using the analytical model where the muon’s will travel through more rock cover and this will enable us to go to probe charge ratio at much higher energy range. Muon charge ratio estimation using both the approach for vertical as well as larger zenith angle can be used for estimating the neutrino flux in the lower and higher energy range at the detector.

Data Availability

the data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work is partially supported by Department of Physics, Lucknow University, Department of Atomic Energy, Harish-Chandra Research Institute, Allahabad and INO collaboration. Financially it is supported by Government of India, DST Project no-SR/MF/PS02/2013, Department of Physics, Lucknow University. The authors thank Prof. Raj Gandhi for useful discussion, providing the cluster facility for data generation and accommodation to work in HRI. Real hill topography is provided by IMSC group so we thank them all for giving this useful information and the discussion specially with Meghna, Lakshmi, Sumanta, Prof. D. Indumathi and Prof. M.V.N. Murthy. The authors thank a lot Prof. G. Majumder for code-related discussion and implementation and many clarifications.