# Cosmogenic particles at ground level and their correlations with primary
particles

Hariom Sogarwal  and Prashant Shukla

###### Abstract

The particle showers produced in the atmosphere due to the interactions of
primary cosmic particles require a thorough understanding in the backdrop of
searches for rare interactions. While the showers encompass the physics of
strong, weak and electromagnetic interactions, the very first interactions are
strong interactions producing hadronic showers which could introduce
uncertainties in the estimates of particle yields. In this work, we made a
comprehensive study of air shower simulations using various combinations of
hadronic models and particle transport code of the CORSIKA package. The
hadronic particles, mostly pions and kaons decay to muons which are the most
abundant charged particles on Earth. The primary proton and helium
distributions are taken as power law which are scaled to match the measured
flux in balloon experiments at the top of atmosphere. The shower simulation
includes production, transport, and decays of secondaries up to the ground
level. In this study, we focus on the bulk of the spectra and particles which
is computationally intensive and hence parallel processing of events is done
on computer cluster. We provide a way to normalize the simulation results to
be compared with the ground-based measurements namely, single and multiple
muon yields and their charge ratios as a function of zenith angle and
momentum. This provides a basis for comparisons among the six model
combinations used in this study and the differences are outlined. Most of the
hadronic models in CORSIKA produce the bulk ground based measurements fairly
well. We use one of the best model combinations to quantitatively predict the
absolute and relative yields of various particles at ground level as well as
their correlations with primaries and with each other. The leptonic ratios are
obtained as a function of energy and zenith angle which are important inputs
for the neutrino oscillation physics. These distributions are important both
for the measurements and validation of models.

##  1 Introduction

The cosmic rays consisting of high energy charged particles continuously
bombard the Earth from all directions and have remained one of the most
interesting and active areas of research [1]. The particle showers produced in
the atmosphere due to the interactions of these primary cosmic particles have
provided a natural laboratory of physics of standard and beyond standard
models which led to many discoveries. Accurate knowledge of both the primary
as well as secondary particles is an important prerequisite for search of rare
interactions. The primary cosmic rays consist of protons (90%), helium nuclei
(9%) and other heavy nuclei (1%) which have been produced and accelerated in
cosmological processes to high energies as large as $10^{8}$ TeV. Upon
entering the atmosphere, the primary cosmic particles interact with air
molecules (mainly oxygen and nitrogen nuclei) roughly 12 km above Sea level.
These interactions lead to the production of shower of hadrons, mainly pions
and kaons with a small percentage of heavier particles which move downwards
while losing energy and decaying. Although neutrinos are the most abundant
shower particles at the Earth surface, muons are the most abundant charged
particles and are responsible for half of the terrestrial radiation dose taken
by a body on Earth. Measurements of muon distributions at ground is an active
field and the improved theoretical understanding and detection technologies
are providing more insights in the particle production in showers [2, 3]. In
the year 2023, LHC is planning a dedicated experiment of proton-oxygen
collisions to understand the air shower physics. To specify, it will measure
hadron multiplicity and the fraction of energy that goes into neutral pions in
each interaction, which is required to resolve the ambiguity of mass
composition for ultra-high energy air showers [4].

The shower simulation codes employ hadronic interaction models which have
energy dependent hadronic cross-section calculations as the main body of the
simulation. CORSIKA is a Monte Carlo program for detailed simulation of air
showers initiated by high-energy primary cosmic particles such as protons [5]
using various hadronic interaction models which can be selected in the CORSIKA
options (version 7.7100) [6]. For the hadronic interactions at high energies,
the most popular models are QGSJET and SIBYLL. QGSJET [7, 8, 9] (Quark Gluon
String model with JETs) describes high-energy hadronic interactions using the
quasi-eikonal pomeron parameterization for the elastic hadron-nucleon
scattering amplitude. The hadronization process is treated in the quark gluon
string model. The more improved version is QGSJET-II-04 [10, 11] which
includes pomeron loop and the cross-sections are tuned to LHC data. SIBYLL
[12, 13, 14, 16, 15] is a program developed to simulate hadronic interactions
at extremely high energies based on the QCD mini-jet model. SIBYLL also
activates the inelastic hadronic interaction cross-sections at higher energies
which are based on QCD calculations. We use SIBYLL 2.3 version [16].

For low-energy (below 100 GeV) interactions, the popular packages are GHEISHA,
UrQMD and FLUKA. GHEISHA (Gamma Hadron Electron Interaction SHower code) is a
model which describes hadronic collisions up to some 100 GeV in many
experiments. A detailed description of the physics processes covered by
GHEISHA may be found in Ref. [17]. UrQMD (Ultra-relativistic Quantum Molecular
Dynamics) is another package designed to treat low-energy hadron-nucleus
interactions. A detailed description of this model may be found in Ref. [18,
19]. UrQMD 1.3 is used in CORSIKA to perform the elastic and inelastic
interactions of hadrons below 80 GeV in air. FLUKA (FLUctuating KAscade) [20,
21] is a package of routines which transport energetic particles through
matter by the monte carlo method. In combination with CORSIKA, only the low-
energy hadronic interaction part is used which calculates the inelastic hadron
cross-sections with the components of air and produce secondary particles
including many details of the de-excitation of the target nucleus. CORSIKA is
also used for generating showers for small local arrays [22]. There are also
plans concerning the design and development for the new generation of CORSIKA
for various scientific applications making an efficient use of computational
resources [23].

CORSIKA is widely used for generating muons for detector simulations [24, 25].
There have been several studies to compare and validate different CORSIKA
models with experimental data. In one of the initial studies [26], muon and
neutrino flux ratios were obtained using various models of CORSIKA and
comparison was made with experimental data and other theory models. The work
in Ref. [27] used CORSIKA to obtain the muon momentum distributions at ground
and at an altitude using VENUS and two low energy models GHEISHA and UrQMD
which are compared with the measurements. Recently, work in Ref. [28] used
four high energy models of hadronic interactions in CORSIKA to compare the
distributions of different particles. The study outlines the differences
between different models. The work in Refs. [29, 30] used CORSIKA to study the
variations of shower angle distribution as a function of atmospheric depth
from tens of TeV to PeV energies of protons.

There have been few recent studies of muon observables obtained with CORSIKA.
The CORSIKA is used to study the dependence of the power index of angular
distribution of muons on the energy and atmospheric depth using the QGSJET and
GHEISHA models [31]. Bahmanabadi [32] measured the muon charge ratio in the
low energy range and made a comparison with two low energy models (GHEISHA and
UrQMD) in CORSIKA. Bahmanabadi and Fazlalizadeh [33] studied the east-west
asymmetry as well as the muon charge ratio sensitivity to the low-energy
GHEISHA and UrQMD in combination with the high energy QGSJET-II model. In Ref.
[34], the momentum and angle distributions of muons obtained using CORSIKA are
compared with the analytical formulae.

Many of the studies in the literature focussed on the hadronic models for very
high-energy showers. Also many studies concentrated on the ratios of the
particles rather than absolute yields. We planned our work to make a
comprehensive study of air shower simulations using different combinations of
low and high energy models in CORSIKA. In this study, we focus on the bulk of
the particles at ground level and obtain the absolute yields using normalized
simulation results. We have used data for muon distributions namely, zenith
angle distribution, momentum distributions (at $0^{\circ}$ and $75^{\circ}$)
and charge ratios to compare with six model combinations. The predictions are
made for multiplicity and East-west asymmetry for muons for all model
combinations. The leptonic distributions and ratios form an important input
for neutrino oscillation studies.

Contribution in the different muon observables coming from different energy as
well as zenith angle intervals of primary protons and heliums is obtained to
correlate the ground based measurements with the initial particle direction
and energy. We predict the Zenith angle and momentum distributions of all
other particles and their ratios which could be measurable at ground in
addition to muons. The distributions of other hadronic and leptonic particles
are obtained. Distributions of muons and neutrinos as a function of $\theta$
and $\phi$ is presented to quantify East-west asymmetry.

We perform CORSIKA simulation using various hadronic interaction models and
obtain various distributions of particles at the ground. The low-energy
hadronic interaction models used in this work are GHEISHA, UrQMD and FLUKA,
whereas the high-energy interaction models are QGSJET, QGSJET-II and SIBYLL.
The initial cosmic proton and helium distributions have been taken as power
law which are matched with experimental data from balloon measurements. The
muon distributions obtained are then compared with the experimental data
measured at the ground level and with some of the recent parametrizations. In
view of a dedicated run of proton oxygen collisions planned at LHC, such study
becomes more important and timely.

##  2 Parameterizations of muon distributions

In this section, we describe the parametrizations which are used to compare
with the CORSIKA simulations. The energy distribution of primary cosmic rays
follows power law $E^{-n}$. The pion and the muon distributions also follow
the same power law which is modified in the low-energy region [35] as

$\displaystyle
I(E)=I_{\circ}N(E_{\circ}+E)^{-n}\bigg{(}1+\frac{E}{\epsilon}\bigg{)}^{-1}.$
(2.1)

Here, $N=(n-1)\,(E_{\circ}+E_{c})^{(n-1)}$ is the normalization in which
$E_{c}$ (= 0.5 GeV) is the cut-off value of the data while and $E_{\circ}$ and
$\epsilon$ are parameters. The zenith angle distribution of energy integrated
flux in terms of $I_{\circ}=\Phi(\theta=0)$ is obtained as

$\displaystyle\Phi(\theta)=I_{\circ}\ D(\theta)^{n-1},$ (2.2)

where D($\theta$) is the ratio of pathlength of a muon from inclined direction
to that of a muon from the vertical direction and is given by [35]

$\displaystyle
D(\theta)=\sqrt{\dfrac{R^{2}}{d^{2}}\cos^{2}\theta+2\dfrac{R}{d}+1}-\dfrac{R}{d}\cos\theta.$
(2.3)

Here, $d$, $R$ and $\theta$ are the vertical depth, radius of the Earth, and
the zenith angle respectively. For the case when the Earth is assumed to be
flat, $D(\theta)$ = $1/\cos\theta$, then the Eq. 2.2 leads to

$\displaystyle\Phi(\theta)=I_{\circ}\cos^{n-1}\theta.$ (2.4)

With $n\simeq 3$, this gives the usual $\cos^{2}\theta$ distribution which is
empirically known to describe Zenith angle distribution. The values of the
parameters for muon energy distributions corresponding to Eq. 2.1 are listed
in the Table 1 [35]. Table 2 lists the values of parameters of muon zenith
angle distribution corresponding to Eqs. 2.2 and 2.4. Table 3 lists the values
of parameters of proton and helium energy distributions corresponding to Eqs.
2.1.

Table 1: Parameters of Eq. 2.1 for muon energy distributions [35]. | $I_{\circ}$ | $n$ | $E_{\circ}$ | $1/\epsilon$
---|---|---|---|---
| (m-2s-1sr-1) |  | (GeV) | (GeV)-1
$\mu$ at $0^{\circ}$ | 70.7$\pm$0.2 | 3.01$\pm$0.01 | 4.29$\pm$0.04 | 1/854
Sea level (E$>$0.5 GeV) |  |  |  |
$\mu$ at $75^{\circ}$ | 65.2$\pm$1.5 | 3.00$\pm$0.02 | 23.78$\pm$0.30 | 1/2000
Sea level (E$>$1.0 GeV) |  |  |  |
Table 2: Parameters for Zenith angle distribution functions of muons [35]. Function | $I_{\circ}$ | $n$ | $R/d$
---|---|---|---
$\Phi(\theta)=I_{\circ}D(\theta)^{n-1}$ | 88.0 $\pm$ 2.4 | 3.09 $\pm$ 0.03 | 174 $\pm$ 12
$\Phi(\theta)=I_{\circ}\cos^{n-1}\theta$ | 85.6$\pm$ 2.4 | 3.01 $\pm$ 0.03 | -
Table 3: Parameters of Eq. 2.1 for proton and helium energy distributions. | $I_{\circ}N$ | $n$ | $E_{\circ}$ | $1/\epsilon$
---|---|---|---|---
|  |  | (GeV) | (GeV)-1
Proton (E$>$0.5 GeV) | 16200$\pm$513 | 2.760$\pm$0.032 | 0.52$\pm$0.11 | 0.0
Helium (E$>$0.5 GeV) | 759$\pm$292 | 2.703$\pm$0.002 | 0.31$\pm$0.55 | 0.0

Gaisser had given analytical formula [36] for muon energy distribution
assuming flat Earth and which is valid for high-energy
$E_{\mu}>100/\cos\theta$ GeV, is given by

$\displaystyle{\frac{dN_{\mu}}{dE.d\Omega}}$ $\displaystyle\approx
1400E^{-2.7}_{\mu}/({\rm m^{2}sGeVsr})$ (2.5)
$\displaystyle\Bigg{(}\frac{1}{1+\frac{1.1E\cos\theta}{\epsilon_{\pi}}}+\frac{0.054}{1+\frac{1.1E\cos\theta}{\epsilon_{\kappa}}}\Bigg{)},$

where the two terms in the bracket give the contributions of pions and kaons
in terms of the two parameters $\epsilon_{\pi}\approx$ 115 GeV and
$\epsilon_{\kappa}\approx$ 850 GeV.

##  3 CORSIKA simulation inputs

In our simulations, the atmospheric model was chosen as the Central European
atmosphere for October month out of the available atmospheric models in the
CORSIKA simulation. The conditions are taken as central Europe which are also
a close approximation to conditions at Tsukuba. The horizontal component
(towards the North) and vertical component (downwards) of the Earth’s magnetic
field are taken as 30.07 $\mu$T and 35.32 $\mu$T respectively. The observation
level is a sphere following the Earth’s curvature at an altitude of 30 m. The
muon multiple scattering angles is selected for large steps by Moliere’s
theory and for small steps by adding many single Coulomb scattering events.

The primary proton/helium energy distribution is chosen as the power law of
the form

$\displaystyle I=A\,E^{-2.7}.$ (3.1)

The energy is generated in the range between 10 GeV and $5\times 10^{5}$ GeV
in our simulations. The simulation process is computationally intensive hence
parallel processing of events at computer cluster has been performed. The
primary particle is generated in nine intervals of energy given by [10, 50]
GeV, [50, $10^{2}$] GeV, [$10^{2}$, $5\times 10^{2}$] GeV, [$5\times 10^{2}$,
$10^{3}$] GeV, [$10^{3}$, $5\times 10^{3}$] GeV, [$5\times 10^{3}$, $10^{4}$]
GeV, [$10^{4}$, $5\times 10^{4}$] GeV, [$5\times 10^{4}$, $10^{5}$] GeV and
[$10^{5}$, $5\times 10^{5}$] GeV. In each energy interval, 100000 events are
generated except in the last interval in which 1000 events are generated.

The spectrum in each energy range is weighted by the factor $f$ given as

$\displaystyle f=\frac{E_{1}^{-1.7}-E_{2}^{-1.7}}{E_{3}^{-1.7}-E_{4}^{-1.7}},$
(3.2)

where $E_{1}$ = 10 GeV, $E_{2}$= 50 GeV and $E_{3}$, $E_{4}$ are the lower and
upper bound of a given energy interval. With this method, then we obtain the
equivalent number of proton showers generated over the whole distribution in
the energy range 10 - $5\times 10^{5}$ GeV as "106,933".

The default primary intensity distribution of protons and heliums with the
zenith angle $\theta$ in CORSIKA goes as

$\displaystyle I(\theta)\sim\,\sin\theta.$ (3.3)

The sine term respects the solid angle element of the sky. The $\phi$
distribution is taken as uniform between $-\pi$ to $\pi$.

Figure 1: (a) The proton flux as a function of momentum falling at the top of
the atmosphere as measured by balloon experiment [37, 38, 39, 41, 42, 43]
along with the initial proton spectrum used in the simulation and the
parameterization given by Eq. 2.1. (b) The zenith angle distribution of
protons at the top of the atmosphere as generated by CORSIKA.

Figure 1(a) shows the proton flux as a function of momentum falling at the top
of the atmosphere as measured by balloon experiment [37, 38, 39, 41, 42, 43]
along with the proton spectrum used in the simulation and the parameterization
given by Eq. 2.1. The momentum distribution of initial protons used in the
simulation is represented here by a histogram having logarithmic bins. The
primary proton spectrum is matched to the measured spectrum and a scale factor
is obtained which is used to obtain normalization for all the spectra obtained
by protons in CORSIKA. Figure 1(b) shows the zenith angle distribution of
protons at the top of the atmosphere as generated by CORSIKA.

Figure 2: (a) The helium flux as a function of momentum per nucleon falling at
the top of the atmosphere as measured by balloon experiment [37, 38, 40, 41,
42, 43] along with the initial helium spectrum used in the simulation and the
parameterization given by Eq. 2.1. (b) The zenith angle distribution of helium
at the top of the as atmosphere generated by CORSIKA.

Figure 2(a) shows the helium flux as a function of momentum per nucleon
falling at the top of the atmosphere as measured by balloon experiment [37,
38, 40, 41, 42, 43] along with the initial helium spectrum used in the
simulation and the parameterization given by Eq. 2.1. The momentum
distribution of initial heliums used in the simulation is represented here by
a histogram having logarithmic bins. The primary helium spectrum is matched to
the measured spectrum and a scale factor is obtained which is used to
normalize all the spectra obtained by heliums in CORSIKA. The normalized
spectra obtained from protons and heliums are then added. Figure 2(b) shows
the zenith angle distribution of heliums at the top of the atmosphere as
generated by CORSIKA.

##  4 Results & discussions

In CORSIKA, among low energy models, we choose GHEISHA, UrQMD and FLUKA and
among high energy models we choose QGSJET, QGSJET-II and SIBYLL. We make six
model combinations as follows:

  1. 1.

[A] GHEISHA & QGSJET

  2. 2.

[B] UrQMD & QGSJET

  3. 3.

[C] UrQMD & QGSJET-II

  4. 4.

[D] GHEISHA & SIBYLL

  5. 5.

[E] UrQMD & SIBYLL

  6. 6.

[F] FLUKA & SIBYLL

We describe our studies in the following three subsections.

###  4.1 Muon observables

Here, we describe the results for the studies done on the muon observables
which are either measured or can be easily measured at ground level detectors.
The zenith angle and momentum spectra of muons obtained by CORSIKA are scaled
by $1/\sin\theta$ to remove the effect of solid angle dependence on zenith
angle $\theta$. The plots are further scaled with the factor
$F=(\pi/2)/\Delta\theta$, where $\Delta\theta$ is the $\theta$ bin size. Here,
$\pi/2$ is the total zenith angle coverage.

We calculate the $\chi^{2}$ (per degree of freedom) for each of the spectra
given by

$\displaystyle\chi^{2}=\frac{1}{n}\sum_{i=1}^{n}w_{i}\Bigg{(}\frac{\phi_{i}^{\rm
data}-\phi_{i}^{\rm corsika}}{\Delta\phi_{i}^{\rm corsika}}\Bigg{)}^{2}.$
(4.1)

Here, $\phi_{i}^{\rm corsika}$ is the value of the flux in CORSIKA and
$\phi_{i}^{\rm data}$ is the same in data for bin $i$. The extrapolated value
between data points is used to match the bins in both cases to obtain the
$\chi^{2}$.

Figure 3 shows the zenith angle distributions calculated with CORSIKA using
various model combinations for muon momentum $p>0.5$ GeV/$c$ along with the
data measured by various experiments. The parameterization given by Eq. 2.2 is
also shown [35]. The data are taken from the collection of Ref. [44] with the
original references [45, 46, 47, 48]. The different datasets have different
muon energy thresholds and we take the normalized data from the review [44].
The plots are scaled by $F=(\pi/2)/\Delta\theta$, where
$\Delta\theta=2^{\circ}$.

The $\chi^{2}$ calculated (with Eq. 4.1) for all the models is given in column
2 of Table 4. All the model combinations produce the shape of zenith angle
distribution given by $\cos^{2}\theta$. Since the Zenith angle distribution is
dominated by low energy muons we can conclude about low energy models. The
combination [D] and [A] of models, where GHEISHA is used as a low energy model
and the combination [F], where FLUKA is used as a low energy model give very
good descriptions of the data. The mismatch at higher angles is due to flat
Earth approximation used in CORSIKA and this mismatch exists for all model
combinations. The parametrization given by Eq. 2.2 produces the shape of the
distribution at all angles.

Figure 3: The zenith angle distributions calculated from CORSIKA using various
model combinations for muon momentum $p>0.5$ GeV/$c$ along with data measured
by various experiments. The parameterization is given by Eq. 2.2 is also shown
[35]. The data are taken from the collection of Ref. [44].

Figure 4: The momentum distribution of vertical atmospheric muons at the
ground level calculated by CORSIKA using different hadronic interaction model
combinations for muon momentum $p>0.5$ GeV/$c$ along with the experimental
data [49, 50]. The parametrizations given by Eq.2.1 [35] and Eq.2.5 [36] are
also shown.

Figure 5: The momentum distribution of atmospheric muons at $75^{\circ}$
zenith angle at the ground level calculated by CORSIKA using different
hadronic interaction model combinations for $p>0.5$ GeV/$c$ along with the
experimental data [51]. The parametrizations given by Eq. 2.1 and Eq. 2.5 are
also shown.

Figure 4 shows the momentum distribution of vertical atmospheric muons at the
ground level calculated by CORSIKA using different hadronic interaction model
combinations for $p>0.5$ GeV/$c$ along with the experimental data [49, 50].
The parametrizations given by Eq.2.1 and Eq.2.5 are also shown. The plots in
Figure 4 are obtained by integrating all muons falling between 0 to 10 degrees
and are scaled by $F=(\pi/2)/\Delta\theta$, where $\Delta\theta=10^{\circ}$.
All model combinations give a good description of experimental data although
all underestimate the muon flux in the very low energy region. The
parametrization by Eq.2.1 gives a good description of the data.

Figure 5 shows the momentum distribution of atmospheric muons at $75^{\circ}$
zenith angle at the ground level calculated by CORSIKA using the different
hadronic interaction model combinations for $p>0.5$ GeV/$c$ along with the
experimental data [51]. The parametrizations given by Eq. 2.1 and Eq. 2.5 are
also shown. The plots in Figure 5 are obtained by integrating all muons
falling between 71 to 81.8 degrees and are scaled by $F=(\pi/2)/\Delta\theta$,
where $\Delta\theta=10.8^{\circ}$. All model combinations give a good
description of experimental data although all underestimate the muon flux in
the very low energy region. The parametrization by Eq.2.1 gives a good
description of the data.

The $\chi^{2}$ calculated (with Eq. 4.1) for all the models corresponding to
Figures 3, 4 and 5 are given in Table 4. From this table it is noted that all
models other than model combination [C] give good values of $\chi^{2}$. The
model QGSJET-II is an improved version of the QGSJET and might be suitable for
ultra-high energies. Overall, the model combinations [E] and [F] give the best
results as compared to the data.

Table 4: Calculated $\chi^{2}$ for zenith angle distribution and momentum
distributions at $0^{\circ}$ and $75^{\circ}$ from CORSIKA models. The points
above muon momentum $p>1$ GeV/$c$ are included in the calculations.

Model | $\chi^{2}$
---|---
Combination | for $\theta$ distribution | for momentum at $0^{\circ}$ | for momentum at $75^{\circ}$
[A] GHEISHA & QGSJET | 1.19 | 1.38 | 3.25
[B] UrQMD & QGSJET | 2.84 | 0.86 | 2.07
[C] UrQMD & QGSJET-II | 4.88 | 2.29 | 6.18
[D] GHEISHA & SIBYLL | 1.10 | 1.45 | 4.52
[E] UrQMD & SIBYLL | 2.34 | 1.02 | 2.30
[F] FLUKA & SIBYLL | 2.04 | 2.14 | 2.31

The integrated muon fluxes for all the model combinations calculated from
zenith angle distribution and from momentum distributions at $0^{\circ}$ and
$75^{\circ}$ are given in the Table 5. The integrals obtained from theta
distribution are multiplied with the $2\pi.(\pi/2)$ to cover the whole solid
angle. We find that the integrated muon flux given by all interaction models
reasonably match with each other except for combination [C] UrQMD & QGSJET-II
which is higher in all 3 cases.

Table 5: Integrated muon flux for different models at ground level for muon
momentum $p>0.5$ GeV/$c$.

| Integrated muon flux
---|---
Model | for $\theta$ distribution | for momentum at $0^{\circ}$ | for momentum at $75^{\circ}$
Combination | (m-2sec-1) | (m-2sec-1sr-1) | (m-2sec-1sr-1)
[A] GHEISHA & QGSJET | 184.80 | 75.55 | 1.32
[B] UrQMD & QGSJET | 191.49 | 73.20 | 1.28
[C] UrQMD & QGSJET-II | 197.59 | 80.38 | 1.35
[D] GHEISHA & SIBYLL | 183.29 | 74.97 | 1.33
[E] UrQMD & SIBYLL | 190.00 | 73.19 | 1.25
[F] FLUKA & SIBYLL | 189.85 | 80.74 | 1.26

Figure 6: The charge ratio $\mu^{+}/\mu^{-}$ as a function of momentum of
atmospheric muons at the ground level calculated by CORSIKA using different
hadronic interaction model combinations along with the experimental data
collection [3].

Figure 6 shows the charge ratio $\mu^{+}/\mu^{-}$ as a function of the
momentum of atmospheric muons at the ground level calculated by CORSIKA using
the different hadronic interaction model combinations along with the
experimental data collection [3]. The original references of experimental data
are BESS [49], CMS [52], MINOS [53], OPERA [54] and L3+C [55]. The combination
of models [F] FLUKA $\&$ SIBYLL gives the best results for the muon charge
ratio in the given momentum range as compared to other model combinations.

Figure 7: The charge ratio $\mu^{+}/\mu^{-}$ as a function of zenith angle of
atmospheric muons for momentum $p>0.5$ GeV/$c$ at the ground level calculated
by CORSIKA using different hadronic interaction model combinations, [A]
GHEISHA & QGSJET, [B] UrQMD & QGSJET, [C] UrQMD & QGSJET-II, [D] GHEISHA &
SIBYLL, [E] UrQMD & SIBYLL, [F] FLUKA & SIBYLL. Integer numbers are added to
models other than [A] as shown to separate the results.

Figure 7 shows the charge ratio $\mu^{+}/\mu^{-}$ as a function of zenith
angle of atmospheric muons for momentum $p>0.5$ GeV/$c$ at the ground level
calculated by CORSIKA using different hadronic interaction model combinations,
[A] GHEISHA & QGSJET, [B] UrQMD & QGSJET, [C] UrQMD & QGSJET-II, [D] GHEISHA &
SIBYLL, [E] UrQMD & SIBYLL, [F] FLUKA & SIBYLL. Integer numbers are added to
models other than [A] as shown to separate the results. Most of the models
give charge ratio around 1.2 and which remains almost flat as a function of
zenith angle but slightly increases towards higher angles (60-80 degrees) for
the models [A], [C], [D] and [E]. The results can be compared with measured
data.

Figure 8: The multiplicity of muon in a shower calculated for muon momentum
$p>0.5$ GeV/$c$ using the (a) [F] FLUKA $\&$ SIBYLL models of CORSIKA. (b)
Multiplicity ratio of model combinations ([A] GHEISHA & QGSJET, [B] UrQMD &
QGSJET, [C] UrQMD & QGSJET-II, [D] GHEISHA & SIBYLL and [E] UrQMD & SIBYLL)
with respect to [F] FLUKA $\&$ SIBYLL model. An integer number is added to
each model result as shown in figure.

Figure 8 shows the multiplicity of muons in a shower calculated for muon
momentum $p>0.5$ GeV/$c$ using (a) [F] FLUKA $\&$ SIBYLL model of CORSIKA. (b)
Shows the multiplicity ratios of model combinations ([A] GHEISHA & QGSJET, [B]
UrQMD & QGSJET, [C] UrQMD & QGSJET-II, [D] GHEISHA & SIBYLL and [E] UrQMD &
SIBYLL) with respect to [F] FLUKA $\&$ SIBYLL model. An integer number is
added to each model result as shown in figure. It can be seen from the ratios
that all models gives similar shape of multiplicity distributions. Models [A]
and [C] give higher values of multiplicities as compared to other models.

###  4.2 Correlations of muons and primary particles

Contributions in the different muon observables coming from different energy
as well as zenith angle intervals of primary protons and heliums are obtained.
Figures 9 to 13 are calculated by CORSIKA using [F] SIBYLL $\&$ FLUKA model.

Figure 9: The zenith angle distribution of muons with momentum $p>0.5$ GeV/$c$
at the ground level. Contributions in total flux are shown from various energy
ranges for the primary particles, (a) Protons and (b) Heliums. Note that the
energy intervals are not the same but increase with energy. The relative
contribution of each bin is normalized as described in section 3.

Figure 9 shows the zenith angle distribution of muons with momentum $p>0.5$
GeV/$c$ at the ground level. Contributions in total flux are shown from
various energy ranges for the primary particles, (a) Protons and (b) Heliums.
Note that the energy intervals are not the same but increase with energy. The
energy bins are normalized as described in section 3. The flux especially at
lower angles is dominated by lower energy showers. Higher zenith angle gets
more contribution from higher energies of primary particles.

Figure 10: Vertical muon flux as a function of momentum for muon momentum
$p>0.5$ GeV/$c$ at the ground level. Contributions in the total flux from
various energy ranges are shown for the primary particles, (a) Proton and (b)
Helium.

Figure 10 shows the vertical muon flux as a function of momentum at the ground
level for muon momentum $p>0.5$ GeV/$c$. Contributions in the total flux from
various energy ranges are shown for the primary particles, (a) Protons and (b)
Heliums. Note that the energy intervals are not the same but increase with
energy but they are normalized as per the power law distribution. Contribution
to lower energy of muons distribution comes from lower energy part of primary
particles. The higher energy part of primary particles give contribution to
the higher energy part of the muon distribution and thus there is a crossover
between distributions of two successive energy histograms of primaries.

Figure 11: Muon flux distribution as a function of momentum at the ground
level for $p>0.5$ GeV/$c$ at 75∘ zenith angle. Contributions in the total flux
from various energy ranges are shown for the primary particles, (a) Proton and
(b) Helium.

Figure 11 shows the muon flux distribution as a function of momentum for
$p>0.5$ GeV/$c$ at 75∘ zenith angle. Contributions in the total flux from
various energy ranges are shown for the primary particles, (a) Protons and (b)
Heliums. One can notice that the contribution from higher energy bins
[100-500] GeV is more dominant as compared to the lower energy bins. It shows
that the muons at higher zenith angles get contributions from higher energy
primary particles. One can compare it with Figure 10, where the lower energy
bins were more dominating in the lower energy part of muon distribution.

Figure 12: Multiplicity variation of cosmic muons above momentum $p>0.5$
GeV/$c$ from various energy ranges of the primary particles, (a) Protons and
(b) Heliums. Each energy interval corresponds to 0.1 million events of
primaries.

Figure 12 shows the multiplicity variation of muons above momentum $p>0.5$
GeV/$c$ from various energy ranges of the primary particles, (a) Protons and
(b) Heliums. Each energy interval corresponds to 0.1 million events of
primaries. The hadronic interaction of primary particles at high energy
produces larger number of secondary hadrons which enhance the muon
multiplicity at ground level. This study is useful to correlate the muon
multiplicity at the ground with the primary particle energy.

Figure 13: Muon flux distributions at ground level for muon momentum $p>0.5$
GeV/$c$. The contributions of primary particles from different angle intervals
is shown in (a) Zenith angle distribution of muons and (b) Momentum
distribution of muons.

Figure 13 shows the muon flux distributions at ground level for muon momentum
$p>0.5$ GeV/$c$. The contribution of primary particles from different angle
intervals is shown in (a) Zenith angle distribution of muons and (b) Momentum
distribution of muons. It is clear from Fig. 13(a) that the secondary
particles (muons in this case) keep the memory of the direction of the primary
particles and Fig. 13(b) shows that the contribution to the lower energy part
of the muon distributions comes from primary particles which are at lower
zenith angles.

The explanations to Figure 9 to Fig. 13 are easy to understand and give us a
quantitative understanding of correlation between the primary particle
direction and energy with muon direction and energy.

###  4.3 Flux and ratios of various particles

Figure 14: The flux of dominant cosmic ray shower particles reaching at ground
level for particle momentum $p>0.5$ GeV/$c$ as a function of (a) zenith angle
and (b) momentum of $\nu_{\mu}$, $\bar{\nu_{\mu}}$, $\nu_{e}$,
$\bar{\nu_{e}}$, $\mu^{+}$, $\mu^{-}$, $n$ and $p$.

We also looked at the distributions of other particles and their contributions
at ground except the electromagnetic part. Figures 14 to 18 are calculated by
CORSIKA using [F] SIBYLL $\&$ FLUKA model. Figure 14 shows the flux of the
dominant cosmic ray shower particles reaching at ground level for particle
momentum $p>0.5$ GeV/$c$ as a function of (a) zenith angle and (b) momentum of
$\nu_{\mu}$, $\bar{\nu_{\mu}}$, $\nu_{e}$, $\bar{\nu_{e}}$, $\mu^{+}$,
$\mu^{-}$, $n$ and $p$. The figure in (a) shows that although the muon flux
falls with increasing angle, the neutrino flux remains flat and then increases
towards higher angles. Neutrons and protons remain the third dominant
component at the ground after neutrinos and muons. The distributions of
neutrons and protons are similar but the number of neutrons is slightly more
than the number of protons.

Figure 15: Flux of the various particles as a function of the zenith angle
from the cosmic ray shower reaching at ground level for particle momentum
$p>0.5$ GeV/$c$: (a) $\mu^{+}$, $\mu^{-}$, $\pi^{+}$, $\pi^{-}$, $K^{+}$,
$K^{-}$ and $K^{\circ}_{L}$ and (b) $\mu^{+}$+$\mu^{-}$, $n$, $p$, $\bar{p}$
and $\bar{n}$.

Figure 15 shows the flux of the various particles as a function of the zenith
angle from the cosmic ray shower reaching at ground level for particle
momentum $p>0.5$ GeV/$c$: (a) $\mu^{+}$, $\mu^{-}$, $\pi^{+}$, $\pi^{-}$,
$K^{+}$, $K^{-}$ and $K^{\circ}_{L}$ and (b) $\mu^{+}$+$\mu^{-}$, $n$, $p$,
$\bar{p}$ and $\bar{n}$. These distributions are important and if compared
with the measurements will help in validation of models.

Table 6 shows the integrated particles flux from zenith angle distribution at
ground level for particle momentum $p>0.5$ GeV/$c$.

Table 6: Integrated particles flux obtained from zenith angle distribution at
ground level for particle momentum $p>0.5$ GeV/$c$.

Particles (at G.L) | Area of $\theta$ distribution | Particles (at G.L) | Area of $\theta$ distribution
---|---|---|---
| (m-2sec-1) |  | (m-2sec-1)
$\nu_{\mu}$ | 1657.97 | $\pi^{+}$ | 0.0237
$\bar{\nu_{\mu}}$ | 1590.28 | $\pi^{-}$ | 0.0178
$\nu_{e}$ | 802.83 | $\bar{p}$ | 0.0040
$\bar{\nu_{e}}$ | 679.52 | $\bar{n}$ | 0.0034
$\mu^{+}$ | 106.12 | $K^{+}$ | 0.0001
$\mu^{-}$ | 85.07 | $K^{-}$ | 0.0001
$n$ | 6.07 | $K^{\circ}_{L}$ | 0.0012
$p$ | 2.26 | $K^{\circ}_{S}$ | 8.2244$\times 10^{-7}$

Figure 16: Ratio of various leptons for momentum $p>0.5$ GeV/$c$ at ground
level (a) zenith angle distribution and (b) momentum distribution. Figure 17:
Ratio of numbers of various leptons for momentum $p>0.5$ GeV/$c$ at ground
level (a) zenith angle distribution and (b) momentum distribution.

Figure 16 and Figure 17 show the ratio of numbers of various leptons at ground
level for momentum $p>0.5$ GeV/$c$ as (a) zenith angle distribution and (b)
momentum distribution. The $\nu_{\mu}/\overline{\nu_{\mu}}$ ratio remains flat
over angle. The ratio $\nu_{e}/\overline{\nu_{e}}$ is similar to the ratio
$\mu^{+}/\mu^{-}$. Thus, measurement of ratio $\mu^{+}/\mu^{-}$ is useful to
get the neutrino ratios in both flavours ($e$ and $\mu$). The ratio
$(\nu_{\mu}+\overline{\nu_{\mu}})/(\nu_{e}+\overline{\nu_{e}})$ decreases with
increasing zenith angle since muons at higher zenith angles have more
probability to decay giving more number of electron neutrinos. Figure 17 shows
the ratio $(\nu_{\mu}+\overline{\nu_{\mu}})/(\nu_{e}+\overline{\nu_{e}})$
increases with increasing energy since muons at higher energy are less
probable to decay giving less number of electron neutrinos.

Figure 18: East-west asymmetry at ground level for particle momentum $p>0.5$
GeV/$c$ for (a) Neutrinos ($\nu_{\mu}$, $\bar{\nu_{\mu}}$, $\nu_{e}$ and
$\bar{\nu_{e}}$) and (b) Muons ($\mu^{+}$ and $\mu^{-}$).

Figure 18 shows the East-west asymmetry at ground level for particle momentum
$p>0.5$ GeV/$c$ for (a) Neutrinos ($\nu_{\mu}$, $\bar{\nu_{\mu}}$, $\nu_{e}$
and $\bar{\nu_{e}}$) and (b) Muons ($\mu^{+}$ and $\mu^{-}$). Asymmetry is
visible in fig (a) of neutrinos at higher zenith angles while it is not
clearly seen in the case of fig (b) of muons. At the ground, more particles
will reach from the west as compared to east direction because the shower of
the secondary particles follows the direction of the primary cosmic particles
(proton and alpha etc). These primaries are positively charged and hence are
deflected by Earth’s magnetic field in the west direction. It was measured by
Johnson and Street in 1933.

##  5 Conclusions

We presented a comprehensive study of air shower simulations using various
combinations of hadronic models in the CORSIKA package. Six model combinations
are used one from the low energy models GHEISHA, UrQMD and FLUKA and one from
the high-energy models QGSJET, QGSJET-II and SIBYLL. In this study, we focus
on the bulk of the particles at ground level and obtain the absolute yields
using normalized simulation results. We have used data for muon distributions
namely, zenith angle distribution, momentum distributions (at $0^{\circ}$ and
$75^{\circ}$) and muon charge ratios and recent parametrizations to compare
with six model combinations. These distributions from various models are
quantified using the $\chi^{2}$ method and integrated flux is obtained from
all the distributions can be compared among the models. All the model
combinations produce the experimental shape of zenith angle distribution and
underestimate the flux at very high zenith angles. Overall, GHEISHA and FLUKA
give the best results for zenith angle distribution. All model combinations
give a good description of momentum distribution although in the low energy
region, they underestimate the flux. Overall the model combinations, UrQMD &
SIBYLL and FLUKA & SIBYLL give best results. The charge ratio of muons as a
function of the momentum is obtained with the different hadronic models and
compared with the experimental data. The combination of models FLUKA & SIBYLL
gives the best results for the muon charge ratio in the given momentum range
as compared to other model combinations. The zenith angle dependence of the
charge ratio of muons is predicted and compared among different models. The
multiplicity distributions from different models have been obtained. The
differences in the results among various models have been outlined.

Contribution in the different muon observables coming from different energies
as well as zenith angle intervals of primary protons and heliums is obtained.
Higher zenith angle gets more contribution from higher energies of primary
particles. The muons at ground keep the memory of the direction of the primary
particles.

We predict the Zenith angle and momentum distribution of all other particles
$\nu_{\mu}$, $\bar{\nu_{\mu}}$, $\nu_{e}$, $\bar{\nu_{e}}$, $\mu^{+}$,
$\mu^{-}$, $n$, $p$, $\pi$ and $K$ and their ratios which could be measurable
at ground in addition to muons. Neutrons and protons remain the third dominant
component at the ground after neutrinos and muons. These distributions are
important and if compared with the measurements will lead to the validation of
models. Various leptonic ratios are obtained as a function of zenith angle and
momentum. These are important for a quantitative understanding of fluxes of
muons and neutrino flavours and their correlations which is an important
prerequisite for neutrino oscillation physics. Distribution of muons and
neutrinos as a function of $\theta$ and $\phi$ are presented to quantify East-
west asymmetry.

## Acknowledgements

The authors acknowledge fruitful discussions and help in simulations from Dr.
Vineet Kumar, Mr. Raman Sehgal and Sundaresh Sankrith.

## References

  * [1] P. K. F. Grieder, Cosmic rays at earth: Researcher’s reference manual and data book, Elsevier (2001).
  * [2] H. Sogarwal and P. Shukla, JCAP 07 11, (2022).
  * [3] R. L. Workman et al. [Particle Data Group], PTEP 083C01 (2022).
  * [4] H. P. Dembinski, R. Ulrich and T. Pierog, "Future Proton-Oxygen Beam Collisions at the LHC for Air Shower Physics", PoS ICRC2019, 235 (2020).
  * [5] D. Heck, G. Schatz, T. Thouw, J. Knapp, J. N. Capdevielle, CORSIKA: A Monte Carlo code to simulate extensive air showers, FZKA-6019 (1998).
  * [6] D. Heck, T. Pierog, Extensive air shower simulation with CORSIKA: A users guide (Version 7.7100) (2016).
  * [7] N. N. Kalmykov, and S. S. Ostapchenko Phys. Atom. Nucl., 56, 346. Yad. Fiz. 56N3 (1993).
  * [8] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, (1994), Bull. Russ. Acad. Sci. Phys., 58, (1966).
  * [9] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, Nucl. Phys. Proc. Suppl., 52, 17 (1997).
  * [10] S. Ostapchenko, Phys. Rev. D 83, 014018 (2011).
  * [11] S. Ostapchenko, Phys. Rev. D 89, 074009 (2014).
  * [12] R. S. Fletcher, T. K. Gaisser, P. Lipari, and T. Stanev, Phys. Rev. D 50, 5710 (1994).
  * [13] J. Engel, T. K. Gaisser, T. Stanev and P. Lipari, 1992, Phys. Rev. D 46, 5013 (1992).
  * [14] Ahn, Eun-Joo, R. Engel, T. K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 80, 094003 (2009).
  * [15] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, T. Stanev,, PoS, ICRC2015, 558 (2016).
  * [16] F. Riehn, Hadronic multiparticle production with Sibyll, KIT, Karlsruhe, Dept. Phys. https://publikationen.bibliothek.kit.edu/1000052699 (2015).
  * [17] H. Fesefeldt, The simulation of hadronic showers: Physics and Applications, PITHA-85-02 (1985).
  * [18] S. A. Bass, et. al, Prog. Part. Nucl. Phys., 41, 255 (1998).
  * [19] M. Bleicher, et. al, J. Phys. G 25, 1859 (1999).
  * [20] T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov and V. Vlachoudis, Nucl. Data Sheets 120, 211 (2014).
  * [21] A. Ferrari, P. R. Sala, A. Fasso and J. Ranft, CERN-2005-10, (2005).
  * [22] T. Wibig, Chin. Phys. C 45, 085001 (2021).
  * [23] R. Engel, D. Heck, T. Huege, T. Pierog, M. Reininghaus, F. Riehn, R. Ulrich, M. Unger and D. Veberič, Comput. Softw. Big Sci. 3 no.1,2 (2019).
  * [24] A. S. Barber, D. B. Kieda, Springer, R. W., Proceedings, 35th International Cosmic Ray Conference ICRC2017, (2017).
  * [25] S. Mueller, and M. Roth, PoS, ICRC2015, 419 (2016).
  * [26] J. Wentz, I. M. Brancus, A. Bercuci, D. Heck, J. Oehlschlager, H. Rebel and B. Vulpescu, Phys. Rev. D 67, 073020, (2003).
  * [27] T. Djemil, R. Attallah and J. N. Capdevielle, Int. J. Mod. Phys. A 20 6950 (2005).
  * [28] R. V. Nikolaenko, A. G. Bogdanov, R. P. Kokoulin and A. A. Petrukhin, Phys. At. Nucl. 84, 1011 (2021).
  * [29] S. M. H. Halataei, M. Bahmanabadi, M. K. Ghomi and J. Samimi, Phys. Rev. D 77, 083001 (2008).
  * [30] P. Patgiri, D. Kalita and K. Boruah, Indian J. Phys. 91, 351 (2017).
  * [31] M. Bahmanabadi and L. Rafezi, Phys. Rev. D 98, 103003 (2018).
  * [32] M. Bahmanabadi, Nucl. Instrum. Meth. A 945, 162635 (2019).
  * [33] M. Bahmanabadi and M. Fazlalizadeh, Phys. Rev. D 100, 083004 (2019).
  * [34] A. Cohu, M. Tramontini, A. Chevalier, J. Marteau and J. C. Ianigro, JAIS 2022, 250 (2022).
  * [35] P. Shukla and S. Sankrith, Int. J. Mod. Phys. A 33, 1850175 (2018).
  * [36] Gaisser, T. K., Astropart. Phys., 16, 285 (2002).
  * [37] K. Abe, et. al, Astrophys. J.822, 65 (2016).
  * [38] O. Adriani et al. (PAMELA), Science332, 69 (2011).
  * [39] M. Aguilar et al. (AMS), Phys. Rev. Lett.114, 171103 (2015).
  * [40] M. Aguilar et al. (AMS), Phys. Rev. Lett.115, 211101 (2015).
  * [41] K. Asakimori et al., Astrophys. J. 502, 278 (1998).
  * [42] A.D. Panov et al. (ATIC Collab.), Bull. Russian Acad. of Science, Physics, 73, 564 (2009).
  * [43] V.A. Derbina, et. al (RUNJOB), Astrophys. J. 628, L41 (2005).
  * [44] S. Cecchini, and M. Spurio, 2012, arXiv:1208.1171 (astro-ph).
  * [45] J. N. Crookes, and B. C. Rastin, Nucl. Phys. B 39, 493 (1972).
  * [46] A. N. Dmitrieva, et. al, Phys. Atom. Nucl. 69, 865 (2006).
  * [47] R. W. Flint, and R. B. Hicks, and S. Standil Can. J. Phys. 50, 843 (1972).
  * [48] M. Gettert, and J. Unger and R. Trezeciak, and J. Engler and J. Knapp, Proc. 23rd International Cosmic Ray Conference (ICRC), Calgary, 4, 394 (1993).
  * [49] Haino, S. et. al, Phys. Lett. B 594, 35 (2004).
  * [50] B. C. Rastin, J. Phys. G 10, 1609 (1984).
  * [51] H. Jokisch, K. Carstensen, W. D. Dau, H. J. Meyer and O. C. Allkofer, Phys. Rev. D 19, 1368 (1979).
  * [52] V. Khachatryan et al. [CMS], Phys. Lett. B 692, 83 (2010).
  * [53] P. Adamson et al. [MINOS], Phys. Rev. D 76, 052003 (2007).
  * [54] N. Agafonova et al. [OPERA], Eur. Phys. J. C 67, 25 (2010).
  * [55] P. Achard et al. [L3], Phys. Lett. B 598, 15 (2004).