# Atmospheric Muon Charge Ratio Analysis at the INO-ICAL Detector

Authors: Jaydip Singh; Jyotsna Singh
Journal: Advances in High Energy Physics
 (2019)

---

## Abstract

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.


---

## Body

## 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 about1017 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 (1)π-→μ-+ν¯μπ+→μ++νμK-→μ-+ν¯μK+→μ++νμ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 →Λ+K+). 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 Section2, 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.Figure 1
Schematic view of the INO-ICAL detector geometry.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 Figure2, at Bodi West Hills in Theni district at Tamilnadu.Figure 2
3D plot of mountain surrounding (range 10 km ) around the INO-ICAL cavern location.The latitude and longitude at the cavern location are as follows:(i)
Latitude:77∘15′0′′E-77∘30′0′′E(ii)
Longitude:9∘57′30′′N-10∘0′0′′NThe 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 Figure2. For the detailed analysis of underground muon charge ratio at INO-ICAL detector, standard rock overburden of density 2.65 gm/cm3 is used in our work. The rock type around the INO location is charnockite whose density is around 2.6 gm/cm3. The real rock density is roughly equal to the standard rock density.


### 2.3. Muon Energy Loss Analysis through Standard Rock

The muon charge ratioRμ+/μ- 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 Eμ,0 at the Earth’s surface lose energy [6] as they traverse a slant depth X (g/cm2) through the Standard rock to reach the detector. The muon energy loss is described by the well-know Bethe-Bloch formula:(2)-dEdX=aEμ+bEμEμ where a parameterizes the energy loss by collisions and b parameterizes the energy loss by three radiative processes, namely, Bremsstrahlung, Pair production, and Photo-nuclear interactions. The energy loss parameters a and b 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 Eμ,0 by [6](3)Eμ=Eμ,0+abebX-ab. The values of the energy loss parameters a and b 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 δ(rμ) in charge asymmetry can be written as(4)δrμ~3.7δEμEμwhereδEμ 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.Figure 3
Muon energy loss versus range in standard rock, data (blue line data) dirctly taken from the PDG(Particle Data Group) web page[6] and the red line values are obtained from our simulation using Geant4, is shown for comparison. The threshold for the vertical muon hitting at the top surface at the INO location is 1.6 TeV.


## 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.Figure 1
Schematic view of the INO-ICAL detector geometry.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 Figure2, at Bodi West Hills in Theni district at Tamilnadu.Figure 2
3D plot of mountain surrounding (range 10 km ) around the INO-ICAL cavern location.The latitude and longitude at the cavern location are as follows:(i)
Latitude:77∘15′0′′E-77∘30′0′′E(ii)
Longitude:9∘57′30′′N-10∘0′0′′NThe 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 Figure2. For the detailed analysis of underground muon charge ratio at INO-ICAL detector, standard rock overburden of density 2.65 gm/cm3 is used in our work. The rock type around the INO location is charnockite whose density is around 2.6 gm/cm3. The real rock density is roughly equal to the standard rock density.


## 2.3. Muon Energy Loss Analysis through Standard Rock

The muon charge ratioRμ+/μ- 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 Eμ,0 at the Earth’s surface lose energy [6] as they traverse a slant depth X (g/cm2) through the Standard rock to reach the detector. The muon energy loss is described by the well-know Bethe-Bloch formula:(2)-dEdX=aEμ+bEμEμ where a parameterizes the energy loss by collisions and b parameterizes the energy loss by three radiative processes, namely, Bremsstrahlung, Pair production, and Photo-nuclear interactions. The energy loss parameters a and b 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 Eμ,0 by [6](3)Eμ=Eμ,0+abebX-ab. The values of the energy loss parameters a and b 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 δ(rμ) in charge asymmetry can be written as(4)δrμ~3.7δEμEμwhereδEμ 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.Figure 3
Muon energy loss versus range in standard rock, data (blue line data) dirctly taken from the PDG(Particle Data Group) web page[6] and the red line values are obtained from our simulation using Geant4, is shown for comparison. The threshold for the vertical muon hitting at the top surface at the INO location is 1.6 TeV.


## 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 χ2/ndf< 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 (ϵrec) is defined as the ratio of the number of reconstructed events, nrec, to the total number of generated events, Ntotal. We have (5)ϵrec=nrecNtotal with uncertainty,(6)δϵrec=ϵrec1-ϵrecNtotal Figure 4 shows the muon momentum reconstruction efficiency as a function of input momentum for different cos⁡θ values. The left and right panels demonstrate detector response for low and high energy muon momentum respectively.Figure 4
Muon momentum reconstruction efficiency as a function of the input momentum for different zenith angles. The left panel refers to low energy, 1≤Pin≤20GeV muons, while the right panel refers to the high energy, 20≤Pin≤500 GeV muons.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 Figure4, 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. Figure5 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 (Pin-Prec) 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].Figure 5
Reconstructed muon momentum distribution for four different values(20, 40, 50 and 100 GeV) of the input muon momentum and at fix zenith angle,cos⁡θ=0.65.Figure 6
Shift in the mean of the reconstructed momentum as a function of the input momentum.


### 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,ncid to the total number of reconstructed events, nrec of same charge, i.e.,(7)ϵcid=ncidnrec with uncertainty,(8)δϵcid=ϵrec1-ϵcidnrec.Figure7 shows the relative charge identification efficiency as a function of input momentum for different cos⁡θ 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.Figure 7
The relative charge identification efficiency as a function of the input muon momentum for differentcos⁡θ values. The left panel refers to low energy, 1≤Pin≤20 GeV muons while the right panel refers to the high energy, 20≤Pin≤500 GeV muons.


## 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 χ2/ndf< 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 (ϵrec) is defined as the ratio of the number of reconstructed events, nrec, to the total number of generated events, Ntotal. We have (5)ϵrec=nrecNtotal with uncertainty,(6)δϵrec=ϵrec1-ϵrecNtotal Figure 4 shows the muon momentum reconstruction efficiency as a function of input momentum for different cos⁡θ values. The left and right panels demonstrate detector response for low and high energy muon momentum respectively.Figure 4
Muon momentum reconstruction efficiency as a function of the input momentum for different zenith angles. The left panel refers to low energy, 1≤Pin≤20GeV muons, while the right panel refers to the high energy, 20≤Pin≤500 GeV muons.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 Figure4, 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. Figure5 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 (Pin-Prec) 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].Figure 5
Reconstructed muon momentum distribution for four different values(20, 40, 50 and 100 GeV) of the input muon momentum and at fix zenith angle,cos⁡θ=0.65.Figure 6
Shift in the mean of the reconstructed momentum as a function of the input momentum.


## 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,ncid to the total number of reconstructed events, nrec of same charge, i.e.,(7)ϵcid=ncidnrec with uncertainty,(8)δϵcid=ϵrec1-ϵcidnrec.Figure7 shows the relative charge identification efficiency as a function of input momentum for different cos⁡θ 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.Figure 7
The relative charge identification efficiency as a function of the input muon momentum for differentcos⁡θ values. The left panel refers to low energy, 1≤Pin≤20 GeV muons while the right panel refers to the high energy, 20≤Pin≤500 GeV muons.


## 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 (E-2.7). Proton flux is generated at 0∘ 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.Figure 8
Differential Muon flux generated using CORSIKA at the top of the hill as a function of muon momentum at log scale.Figure 9
Muon flux generated using CORSIKA at the top of the hill (i.e. 2 k.m. above the sea level) with primary protons in energy range 1-20 TeV(brown dots) and its comparison with the existing experimental fluxes (★ [2],△[7], +[8], □[9], •[10]) and the flux (solid blue line) from [11].Figure 10
The distribution of secondary muons produced 2.0 Km. above the sea level from the vertical primary proton shower, in X-Y plane.Density of the muons near the shower axis is very high as shown in Figure10 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,Rμ+/μ-≡Nμ+/Nμ-, 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 Figure11. 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.Figure 11
Differential muon flux(red lineμ+ and blue line μ-) at the ICAL detector as a function of reconstructed muon momentum after passing through the rock.The vertical muon charge ratio has been evaluated up to 2 TeV and plotted in Figure12 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.Table 1
Data Table for Charge Ratio Analysis from Vertical Muons.
Energy (TeV) R μ + / μ - at ICAL R μ + / μ - S u r f a c e C O R S I K A R μ + / μ - S u r f a c e P I K A 1.60-1.65 1.32±0.013 1.31 1.41 1.65-1.70 1.31±0.022 1.34 1.41 1.70-1.75 1.32±0.014 1.34 1.41 1.75-1.80 1.32±0.034 1.37 1.41 1.80-1.85 1.41±0.008 1.40 1.41 1.85-1.90 1.36±0.008 1.37 1.41 1.90-1.95 1.41±0.026 1.37 1.41 1.95-2.00 1.32±0.030 1.36 1.41  Rμ+/μ- 1.35±0.019 1.36 1.41Figure 12
Muon charge ratio with statistical errors at the INO-ICAL detector as a function of surface muon energy, obtained using the CORSIKA data.


### 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 ofEμcos⁡θ. 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(9)dNμdEμ=0.14Eμ-2.7cm2ssrGeV×11+1.1Eμcos⁡θ/ϵπ+η1+1.1Eμcos⁡θ/ϵKThis formula is valid when muon decay is negligible(Eμ>100/cos⁡θ GeV) and the curvature of the Earth is neglected (θ<700). 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 K+/K- 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 fπ and fK,(10)dNμ+dEμ=0.14Eμ-2.7cm2ssrGeV×fπ1+1.1Eμcos⁡θ/115GeV+η×fK1+1.1Eμcos⁡θ/850GeV(11)dNμ-dEμ=0.14Eμ-2.7cm2ssrGeV×1-fπ1+1.1Eμcos⁡θ/115GeV+η×1-fK1+1.1Eμcos⁡θ/850GeVIn the above equationϵπ and ϵK 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 π/K 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:(12)Rμ=fπ/1+1.1Eμcos⁡θ/115GeV+η×fK/1+1.1Eμcos⁡θ/850GeV1-fπ/1+1.1Eμcos⁡θ/115GeV+η×1-fK/1+1.1Eμcos⁡θ/850GeVThe charge ratio of muon from pion decay isrπ = fπ/(1-fπ) and from kaon decay is rK=fK/(1-fK). Equation (12) is a function of Eμcos⁡θ only; all other parameters are extrapolated from the available experimental data; the parameters value used in our work are fπ= 0.5510±0.0006 and fK = 0.7006±0.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].Figure 13
Muon charge ratio with systematic errors at the INO-ICAL detector as a function of surface muon energy, obtained using the “pika" model.


## 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 (E-2.7). Proton flux is generated at 0∘ 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.Figure 8
Differential Muon flux generated using CORSIKA at the top of the hill as a function of muon momentum at log scale.Figure 9
Muon flux generated using CORSIKA at the top of the hill (i.e. 2 k.m. above the sea level) with primary protons in energy range 1-20 TeV(brown dots) and its comparison with the existing experimental fluxes (★ [2],△[7], +[8], □[9], •[10]) and the flux (solid blue line) from [11].Figure 10
The distribution of secondary muons produced 2.0 Km. above the sea level from the vertical primary proton shower, in X-Y plane.Density of the muons near the shower axis is very high as shown in Figure10 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,Rμ+/μ-≡Nμ+/Nμ-, 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 Figure11. 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.Figure 11
Differential muon flux(red lineμ+ and blue line μ-) at the ICAL detector as a function of reconstructed muon momentum after passing through the rock.The vertical muon charge ratio has been evaluated up to 2 TeV and plotted in Figure12 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.Table 1
Data Table for Charge Ratio Analysis from Vertical Muons.
Energy (TeV) R μ + / μ - at ICAL R μ + / μ - S u r f a c e C O R S I K A R μ + / μ - S u r f a c e P I K A 1.60-1.65 1.32±0.013 1.31 1.41 1.65-1.70 1.31±0.022 1.34 1.41 1.70-1.75 1.32±0.014 1.34 1.41 1.75-1.80 1.32±0.034 1.37 1.41 1.80-1.85 1.41±0.008 1.40 1.41 1.85-1.90 1.36±0.008 1.37 1.41 1.90-1.95 1.41±0.026 1.37 1.41 1.95-2.00 1.32±0.030 1.36 1.41  Rμ+/μ- 1.35±0.019 1.36 1.41Figure 12
Muon charge ratio with statistical errors at the INO-ICAL detector as a function of surface muon energy, obtained using the CORSIKA data.


## 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 ofEμcos⁡θ. 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(9)dNμdEμ=0.14Eμ-2.7cm2ssrGeV×11+1.1Eμcos⁡θ/ϵπ+η1+1.1Eμcos⁡θ/ϵKThis formula is valid when muon decay is negligible(Eμ>100/cos⁡θ GeV) and the curvature of the Earth is neglected (θ<700). 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 K+/K- 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 fπ and fK,(10)dNμ+dEμ=0.14Eμ-2.7cm2ssrGeV×fπ1+1.1Eμcos⁡θ/115GeV+η×fK1+1.1Eμcos⁡θ/850GeV(11)dNμ-dEμ=0.14Eμ-2.7cm2ssrGeV×1-fπ1+1.1Eμcos⁡θ/115GeV+η×1-fK1+1.1Eμcos⁡θ/850GeVIn the above equationϵπ and ϵK 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 π/K 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:(12)Rμ=fπ/1+1.1Eμcos⁡θ/115GeV+η×fK/1+1.1Eμcos⁡θ/850GeV1-fπ/1+1.1Eμcos⁡θ/115GeV+η×1-fK/1+1.1Eμcos⁡θ/850GeVThe charge ratio of muon from pion decay isrπ = fπ/(1-fπ) and from kaon decay is rK=fK/(1-fK). Equation (12) is a function of Eμcos⁡θ only; all other parameters are extrapolated from the available experimental data; the parameters value used in our work are fπ= 0.5510±0.0006 and fK = 0.7006±0.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].Figure 13
Muon charge ratio with systematic errors at the INO-ICAL detector as a function of surface muon energy, obtained using the “pika" model.


## 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 =105 g/cm2). 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 (XI) and threshold energy of muon to reach the detector after travelling a distance XI. The values of parameters a and b 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 a and b we have estimated the threshold (for muon) for the corresponding slant depth covering in the zenith angle range 00-600. 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 (00-700). In our work charge ratio estimation shown in Figures 13 and 14 corresponds to the zenith angle range 00-600. 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.Figure 14
Muon charge ratio for a wide range of zenith angles(θ=00-600) at INO as a function of surface muon energy and their comparison with standard existing experimental results [2–5, 7]. The charge ratio results estimated from CORSIKA are also shown in this Figure (red blocks) till 2 TeV lying in the ratio range around 1.3 to 1.4.


## 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 Sections3.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 0-600 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 Figure13. The average value of the charge ratio observed at underground detector is 1.347±0.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 ofEμcos⁡θ 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 m2) 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.



---

*Source: 9585234-2019-07-01.xml*