MSc Physics and Astronomy
Gravitation and Astroparticle Physics
Master Thesis
Investigating post-LHC hadronic
interaction models and their
predictions of cosmic ray shower
observables
June 22, 2018
60 ECTS
Author:
Stephan Runderkamp BSc.
Student ID: 10447288
Examiners:
prof. dr. ing. Bob van Eijk
prof. dr. Patrick Decowski
Daily Supervisor:
Kasper van Dam MSc.
Master Thesis Physics and Astronomy, conducted between September 2016 and June 2018
Abstract
There are large systematic uncertainties in the interpretation of the data of cosmic ray experi-ments due to the unreliability of hadronic interaction models in MC simulations at high cosmic ray energies.
In this work the predicted values of the number and energy of secondaries produced after the first hadronic interaction are compared for the three post-LHC models EPOS LHC, SIBYLL 2.3c and QGSJET-II-04. This study focuses on proton-nitrogen, alpha-nitrogen and iron-nitrogen
collisions with primary energies up to 1020eV.
Large shower-to-shower fluctuations are found, but also significant differences in the predicted values between the three models. An important finding is the significant difference in the pre-dicted number of produced baryons and mesons and the dependence of the multiplicity on primary energy and primary mass number for the three models.
The position of shower maximum, the number of ground-level particles (for both photons, elec-trons, muons and hadrons) and the energy, arrival time and lateral distance to the shower core at ground level are compared for the three post-LHC models. Vertical proton, alpha and iron
initiated showers with primary energies up to 1016.5eV are considered. Significant differences
between the predictions of the models are already found at primary energies far below LHC energies.
The three models predict the same value of the position of the shower maximum Xmax. Large
deviations, up to 40%, are found (even below LHC energies) for the number of ground-level hadrons. The energy spectra of ground-level muons is found to be the hardest for QGSJET-II-04, followed by SIBYLL 2.3c, followed by EPOS LHC. Consequently, ground-level muons arrive latest for EPOS LHC, followed by SIBYLL 2.3c, followed by QGSJET-II-04. This is relevant for describing the lateral shower profile.
Contents
1 Introduction 4
2 Cosmic rays 6
2.1 Cosmic rays in the atmosphere . . . 8
2.1.1 Electromagnetic shower . . . 8
2.1.2 Hadronic shower . . . 9
2.1.3 Longitudinal profile . . . 10
2.1.4 Lateral profile . . . 14
3 HiSPARC 17 3.1 Charged particles in scintillators . . . 18
3.2 Signal . . . 22
3.3 Direction reconstruction . . . 23
3.4 Energy reconstruction . . . 27
4 Hadronic interaction models 29 4.1 Post-LHC hadronic interaction models . . . 30
4.1.1 EPOS LHC v3400 . . . 30
4.1.2 QGSJET-II-04 . . . 31
4.1.3 SIBYLL 2.3c . . . 33
4.1.4 Summary . . . 34
5 Analysis 34 5.1 First hadronic interaction . . . 34
5.1.1 Number of secondaries . . . 35
5.1.2 Energy of secondaries . . . 49
5.2 Cosmic ray shower observables . . . 52
5.2.1 Shower maximum . . . 53
5.2.2 Number of particles at ground level . . . 54
5.2.2.1 Photons . . . 54
5.2.2.2 Electrons . . . 56
5.2.2.3 Muons . . . 58
5.2.2.4 Hadrons . . . 60
5.2.3 Distance to shower core . . . 62
5.2.4 Arrival time . . . 63
5.2.5 Energy distribution . . . 64
7 Discussion 68
1
Introduction
One of the main challenges of astroparticle physics is to understand the source and propaga-tion of cosmic rays. Cosmic rays are charged particles from extra-terrestrial origin bombarding
Earth’s atmosphere. The energies of these cosmic rays range from less than 1 GeV up to 1020eV.
Cosmic rays can be measured directly by detectors in balloons and satellites. However, above
1015eV the flux of high-energy cosmic rays becomes too low to measure cosmic rays directly. At
energies E ≥ 1014eV cosmic rays can be measured indirectly.
Cosmic rays interact high in the atmosphere with air nuclei (mostly nitrogen and oxygen) and
cosmic rays with energies E ≥ 1014eV initiate cascades, showers of particles, of which a sufficient
amount of particles reaches the ground to be measured by ground-level cosmic ray experiments. These showers are called extensive air showers (EAS).
To be able to determine the properties of the incoming cosmic ray, such as the energy, direction and particle type, measurements of cosmic ray experiments have to be compared with predic-tions of MC simulapredic-tions of showers. Calibration of cosmic ray experiments with high-energy test beams is impossible and therefore we can only compare measured data with data from lots of simulated events. Currently, the largest source of uncertainty of cosmic rays experiments is the unreliability of these MC simulations [1]. This is because at the highest cosmic ray energies,
E ≥ 1017eV, hadronic interactions, and especially the hadronic multiparticle production, are
poorly understood, since these energies are higher than that reached by current collider experi-ments.
There are currently three hadronic interaction models, models that describe hadronic
interac-tions, tuned to LHC data with center-of-mass energies √s = 0.9TeV and 7TeV, therefore called
the post-LHC models: EPOS LHC [2], SIBYLL 2.3c [3, 4] and QGSJET-II-04 [5]. Proton-proton
collision energy √s = 7TeV corresponds to about 2.5 × 1016eV in the laboratory frame. This
means that simulations of hadronic interactions with higher energies are highly uncertain. The predicted values of properties of hadronic interactions at these extrapolated energies depend on the parametrization of the chosen model and the values differ widely as will be shown in this thesis.
My work in this thesis has two goals. One goal is to investigate and compare the model-dependent predicted properties of the secondaries produced in proton-nitrogen, alpha-nitrogen and
iron-nitrogen collisions with lab energies up to the highest cosmic ray energies (E = 1020eV). The
focus will be on the number of produced secondaries and how the energy is distributed between the different produced particles.
The other goal is to investigate and compare the predicted values of cosmic ray observables
for vertical (zenith angle = 0◦) proton, alpha and iron initiated showers with primary energies
up to E = 1016.5eV for the three post-LHC models. Furthermore, it is discussed if and how
the differences between these model-dependent predicted values of shower observables could be measured to investigate which one of the models describes the showers best. Analysis is done
for the following shower observables: position of the shower maximum, the number of particles at ground level (photons, electrons, muons and hadrons) and the arrival times, energy spectra and lateral distance to the shower core of ground-level electrons and muons.
In Section 2 Cosmic Rays will be described in more detail with emphasis on the properties of extensive air showers (EAS). An example of a cosmic ray experiment, the HiSPARC experiment, and its dependence on MC simulations will be described in Section 3. In Section 4 hadronic interactions and the theoretical difficulties to describe these interactions will be explained. Fur-thermore, in this section the different theoretical approaches to describe the hadronic interactions for the three post-LHC models EPOS LHC, QGSJET-II-04 and SIBYLL 2.3c will be described. In Section 5 the results of MC simulations for the three different hadronic interaction models are compared. First, the predicted number of produced secondaries in the first hadronic in-teraction and the predicted energy distribution between these secondaries are analyzed for the three post-LHC hadronic interaction models. This is done for proton-nitrogen, alpha-nitrogen
and iron-nitrogen collisions with primary energies up to E = 1020eV. Further in Section 5 the
model-dependent predicted values of the position of the shower maximum, the number of par-ticles at ground level (photons, electrons, muons and hadrons) and the arrival times, energy spectra and lateral distance to the shower core of ground-level electrons and muons will be compared for the three models. This analysis is performed for vertical proton, alpha and iron
initiated showers with energies up to E = 1016.5eV. The findings will be concluded in Section 6
2
Cosmic rays
Cosmic rays are particles from extra-terrestrial origin bombarding Earth’s atmosphere. The
energies of cosmic rays range from less than 1GeV up to 1020eV, which is orders of magnitude
higher than the energy that is reached in man-made collider experiments. These cosmic rays can improve our understanding of high energy processes in both galactic and extra-galactic astrophysical sources and the propagation of these particles through the Galaxy. It is one of the main challenges of astroparticle physics to understand how these particles are accelerated to such high energies and in what sources, where these sources are (e.g. galactic or extra-galactic) and how they propagate through the insterstellar medium (or even intergalactic medium). Currently, there are different theories for this and these predict different chemical composition changes of cosmic rays over primary energy. Therefore, investigating the chemical composition of cosmic rays for different cosmic ray primary energies is crucial. From direct measurements it is known that the low primary energy cosmic rays (10.6GeV/nucleon) are mainly protons (about 93%) and alpha particles (about 6%) [6]. About one percent of the cosmic rays are heavier nuclei, such as iron, and electrons and gamma rays (high energy photons), although gamma rays are not necessarily a part of cosmic rays, since they are not charged. As already described, this particle composition of cosmic rays is expected to change over primary energy.
The cosmic ray flux depends on primary energy and can be seen in Figure 1. Note that the flux
is scaled with E2.5 to make the differences in the slope visible. The cosmic ray spectrum spans
multiple orders of magnitude, from several MeV’s (rest masses of the particles) up to 1020eV,
and even further. The cosmic ray spectrum follows the power law F (E) ∼ E−γ with spectral
index γ ∼ 2.7. This power law is associated with the acceleration processes in astrophysical sources. However, as can be seen in Figure 1, the spectral index is not constant.
At an energy of ∼ 4 × 1015eV the slope steepens, this is called the knee. There are multiple
explanations for this feature in the spectrum. One of the explanations is that cosmic rays leak away from our Galaxy since they are less contained by the Galactic magnetic fields at these energies. The gyromagnetic radius R for a relativistic charged particle with proton number Z is
proportional to EZ. Therefore, particles with low Z of similar energies have larger gyromagnetic
radii and are more likely to escape our Galaxy and not reach Earth. Another explanation is that the maximum energy to which galactic sources can accelerate particles is reached. Cosmic rays
with energies up to about 1017eV are most likely accelerated by shock acceleration in supernova
remnants. Because these shock waves have a finite life time, also a maximum energy can be reached for the particles. This maximum energy is proportional to the proton number Z. For both explanations a transition in cosmic ray chemical composition from protons to heavier nuclei (such as iron) at these energies is expected. A detailed review of the different theories to describe the knee can be found in [7].
At an energy of about 4 × 1018eV the slope flattens again, which is known as the ankle. Most
the dip model [8] predicts that the slope change comes from energy losses by interactions of protons with CMB photons in extragalactic propagation, thereby producing electron-positron
pairs (p + γCM B → p + e−e+), which can only be true if cosmic rays with primary energies
E ≥ 1018eV are mostly protons from extragalactic sources.
For the highest energies E ≥ 1019eV the cosmic ray flux is suppressed. There are multiple
explanations for this flux suppression. One explanation is that it is an imprint of the interactions
of protons and nuclei with CMB photons, producing pions via the ∆ resonance: γCM B+ p →
∆+→ p + π0 or γCM B+ p → ∆+→ n + π−. This process is called the Greisen-Zatsepin-Kuzmin
(GZK) [9, 10] energy loss effect and this theory only holds if the sources are extragalactic. The
threshold energy for this process is 0.5 × 1020eV×A, with A the nuclear mass of the cosmic ray.
Therefore, this theory predicts extragalactic sources and a transition in cosmic ray chemical composition from protons to heavier nuclei (such as iron) at these extremely high energies. Another explanation is that the maximum energy for (extra-)galactic accelerators is reached. If the changes in cosmic ray chemical composition can not be explained by astrophysical processes, this can be a hint of new particle physics at ultrahigh energies. A detailed review of the different theories for ultrahigh energy cosmic rays can be found in [11].
Figure 1: Cosmic ray spectrum (all particles) with data from different experiments. The flux is scaled with E2.5 to make subtle differences of the slope visible. The c.m. energies at the top-axis for proton-proton collisions correspond to the energies at the bottom-axis for cosmic ray protons. Figure from [12].
2.1 Cosmic rays in the atmosphere
Before galactic cosmic rays reach the Earth’s atmosphere they have traveled through the inter-stellar medium of our Galaxy. In [13] a calculation is done that cosmic rays originating from
the Galactic centre traverse a column density ≤ 10 g/cm2. Extragalactic cosmic rays travelled
through the interstellar medium of the galaxy it is produced in, the intergalactic medium and the interstellar medium of our Galaxy, thereby traversing a higher column density than galactic
cosmic rays. The interaction length of protons is about 90g/cm2, therefore some of the
cos-mic rays reaching Earth have already interacted. In our atmosphere the density of particles is much higher than in the interstellar medium. Our atmosphere has a column density of about
1000g/cm2, which is more than 11 interaction lengths long. Therefore, all of the cosmic rays
interact with air nuclei (mostly nitrogen and oxygen) in the atmosphere before reaching sea
level. We usually don’t refer to column density, but to atmospheric depth X (in g/cm2) at a
given altitude:
X =
Z ∞
z
ρ(h)dh
On average the first hadronic interaction takes place at an altitude of about 20 kilometers, but this also depends on the inclination and mass number of the primary particle. When cosmic rays do not enter the atmosphere vertically (from zenith), but with a certain angle θ, the path length is longer and the particles traverse more column density and will therefore interact higher in the atmosphere. Heavier nuclei, with a higher mass number, have shorter interaction lengths.
The interaction length of iron is ∼ 5g/cm2, for instance. These particles will interact earlier and
therefore higher up in the atmosphere.
2.1.1 Electromagnetic shower
In the first hadronic interaction of the incoming cosmic ray with an air nucleus (mostly nitrogen
and oxygen) a lot of secondary particles are produced, mostly pions. The neutral pion, π0, is
unstable and has a short mean lifetime of 8.52 ± 0.18 · 10−17s and will therefore immediately
decay before interacting. The π0decay mode with the highest branching ratio is the decay
into two photons, π0 → γ + γ , with BR = (98.823 ± 0.034) %. In the vincinity of an atomic
nucleus these photons (and photons of other processes) will produce an electron-positron pair, this process is called pair production:
γ → e−+ e+
Charged light particles like electrons and positrons radiate photons when interacting with elec-tromagnetic fields of nuclei, termed Bremsstrahlung:
Of course there is also the positron equivalent of this process. The produced photons by Bremsstrahlung will create new electron-positron pairs by pair production, and these electrons and positrons which will later radiate photons by Bremsstrahlung, and so on. More and more particles will be produced and an electromagnetic shower with lots of electrons, positrons and photons is formed. The average energy of the particles decreases with every interaction and after some interactions the energy of the particles is so low that no new particles can be produced anymore. The low-energetic electrons and positrons will be scattered out of the shower or lose their remaining energy by ionization and eventually get absorbed by the atomic nuclei. The low energy photons will get absorbed by the photoelectric effect and Compton scattering. This process is known as an electromagnetic shower.
2.1.2 Hadronic shower
In the first hadronic interaction of the incoming cosmic ray with an air nucleus also charged pions
and other particles (but in lower numbers) such as kaons are produced. The charged pions, π±,
have a mean lifetime of 26ns at rest. The interaction length of pions is ∼ 120g/cm2, which means
that only charged pions with energies higher than about 30GeV will interact before decaying, due to their high Lorentz-factors. But this also depends on the altitude where these pions are produced. High in the atmosphere the density is low and particles are more likely to decay before interacting than at lower altitudes. High-energy charged pions interact with air nuclei,
mostly nitrogen and oxygen (π±+ air → . . .). This hadronic interaction produces, depending
on the c.m. energy of the collision, multiple secondaries, mostly pions. With every interaction with air nuclei more and more (charged) pions are created and these form a so-called hadronic shower. The average energy of the charged pions decreases with every hadronic interaction with air and after a while the energies are too low for the charged pions to interact before decaying. The pions decay to muons and neutrinos according:
π−→ µ−+ νµ
π+→ µ++ νµ
The same mechanism is true for kaons, which have similar mean lifetimes and interaction lengths, but these can decay differently. But since kaons are produced in lower amount with respect to pions these processes are less important for understanding the overall shower. Muons are about 200 times heavier than electrons and have a mean lifetime of 2.2µs. Muons with energies of only a few GeV can, due to time dilation, travel several kilometers through the atmosphere, and even reach sea level (the ground). Since the muons are much heavier than electrons the lost energy due to Bremsstrahlung is negligible compared to that of electrons (total radiated power goes
with m−4 or m−6). Muons on average lose about 2GeV to ionization and have a mean energy
and neutrinos, thereby feeding the electromagnetic shower:
µ+ → e++ νe+ νµ
µ− → e−+ νe+ νµ
2.1.3 Longitudinal profile
An important feature of cosmic ray showers that cosmic ray experiments use to determine the properties of cosmic rays, such as primary energy, particle type and inclination, is the longitudinal development of the shower: the longitudinal profile. The longitudinal profile, the
evolution of the number of particles in the atmosphere (over g/cm2), of showers firstly depends
on the altitude of the first interaction, which depends on mass number of the primary and the inclination of the shower (as described in Section 2.1). This is the point where the growth of the number of particles starts. The longitudinal profile depends also on the properties of the first interaction itself: the number of secondaries produced in the first interaction. Therefore, the initial growth of the particles in the longitudinal profile, depends on the energy of the collision and the mass number of the incoming particle (and of the target nucleus). The more secondaries produced, the less average energy per particle, which influences the longitudinal
development of the shower. In the first few interactions, the first few g/cm2 after the first
interaction, the number of particles increases almost exponentially. After every interaction the average energy of the particles decreases and after some interactions more and more low-energy electrons and positrons are absorbed by atomic nuclei or are scattered out of the shower. Also, low energy photons get absorbed by the photoelectric effect and Compton scattering. Therefore,
at a certain point, at the shower maximum Xmax, the number of particles starts decreasing
again. The shower maximum Xmax increases logarithmically with primary energy, therefore
it is a good characteristic of a shower to determine the primary energy. An example of the
longitudinal profile of electrons (e− and e+) of a vertical 1015eV proton shower is shown in
Figure 2. The number of electrons/positrons of this particular example starts increasing from
the first interaction at ∼ 100g/cm2 (∼ 22km) and starts decreasing from the shower maximum
Xmax at about 575g/cm2 (∼ 5km). At the ground (sea level), corresponding to ∼ 1030g/cm2,
Figure 2: Longitudinal profile of electrons (e−and e+) of a vertical 1015eV proton initiated shower. Figure from [14].
However, showers with the same initial conditions are not the same. There are, often large, shower-to-shower fluctuations. For instance, the column depth between two interactions is not always equal to the (average) interaction length, it is subject to probability. Also, the number of produced secondaries per interaction and the energy/momentum distributions of these sec-ondaries are all a matter of probability (following cross-sections). Therefore, the longitudinal development of showers with the same initial conditions may differ widely. For cosmic ray exper-iments it is important that these shower-to-shower fluctuations are known very precisely. This is the reason why for a cosmic ray experiment lots of simulations of showers have to be done, for instance to determine the uncertainty of the experiment. An example of such shower-to-shower
fluctuations is depicted in Figure 3, where the longitudinal profile of electrons (e− and e+) of
50 vertical 1015eV proton showers are plotted. Most of the 50 showers start developing before
200g/cm2 and have Xmax at about 500g/cm2, but there are three showers that develop slower
and have Xmax at ∼ 800g/cm2. The number of electrons/positrons at the ground (sea level)
ranges from ∼ 5 × 104 to ∼ 5 × 105.
Figure 3: Longitudinal profile of electrons (e−and e+) of 50 vertical 1015eV proton initiated showers. Notice the large fluctuations in position of first interaction, position of shower maximum and number of electrons at different atmospheric depths. Figure from [14].
For showers with primary energies lower than 1014eV almost all the particles (mostly photons
the ground. For showers with primary energies higher than 1014eV the number of electrons on the ground exceeds the number of muons. For these showers enough particles reach the ground to be able to observe them with ground-level detectors. These showers are called Extensive Air Showers (EAS). An example of such an Extensive Air Shower is presented in Figure 4, where the
particle tracks of a vertical 1015eV proton shower are plotted. Notice that the first interaction
takes place at an altitude of about 20km. A high number of particles reaches sea level (z = 0km), therefore it is an EAS.
Figure 4: Particle tracks of a vertical 1015eV proton initiated shower. Left: electron (e− and e+) tracks (red). Right: hadron tracks (blue) drawn over the muon (µ−and µ+) tracks (green) drawn over the electron (e−and e+) tracks (red). The number of electrons, muons and hadrons at ground level (z = 0) are about 105, 104 and 2 × 103, respectively. Figure from [14].
The longitudinal profile of different particles in an EAS, in this case a vertical 1015eV proton
shower, is shown in Figure 5. It can be seen that the most abundant particles in an EAS are photons, followed by electrons, followed by muons. The number of muons at this primary energy is quite stable for large atmospheric depths, whereas the number of photons and electrons decreases because of absorption in the atmosphere.
Figure 5: Longitudinal profile of photons, electrons (e− and e+) and muons (µ− and µ+) of a vertical 1015eV proton initiated shower. Figure from [15].
As described earlier, for non-EAS showers, with primary energies lower than 100TeV (1014eV),
most of the electrons and photons will be absorbed in the atmosphere. Therefore, the number of muons at sea level is higher than that of electrons and photons. Since the flux of these low-energy showers is much higher than that of high-energy (extensive air) showers, as can be seen again in Figure 1, the muon flux at sea level is about four times higher than the electron flux. The particle fluxes for muons, electrons and protons as a function of atmospheric depth
are presented in Figure 6. The muon rate at the ground is ∼ 1cm−2min−1.
Figure 6: Particle flux for muons (both µ− and µ+), electrons (both e− and e+) and protons as a function of atmospheric depth. Figure from [15].
2.1.4 Lateral profile
Another important feature of cosmic ray showers that experiments use to determine properties of the cosmic ray is the lateral development of the shower: the lateral profile. With every hadronic and electromagnetic interaction the particles in the shower obtain some transverse momenta. In electromagnetic interactions transverse momentum is obtained by bremsstrahlung and pair production due to the angles between the particles. Particles can also be deflected by Coulomb forces and Compton scattering. The (produced) particles in hadronic collisions also obtain transverse momentum, as do the particles produced by particle decay. Due to these processes the shower widens.
The footprint of an EAS on the ground can be several kilometers wide. The properties of this footprint depend on the primary energy of the cosmic ray, the type of the incoming cosmic ray particle (e.g. proton, alpha, iron) and the inclination of the shower. On the ground the particle density is the highest in the shower core, which follows the axis of the incoming cosmic ray, and this particle density decreases outwards. The number of particles that can travel to the ground increases with the primary energy of the shower, therefore the particle density also increases with primary energy. Cosmic rays with higher mass numbers interact higher in the atmosphere because their interaction length is shorter (i.e. probability to interact is higher).
As a result, Xmax is also lower (higher in the atmosphere) and particles are more attenuated,
therefore decreasing the particle density on the ground. The particle densities on the ground of a proton shower is thus higher than that of an iron or alpha shower with the same primary energies. The particle densities and the form of the footprint on the ground also depend on
the inclination of the shower. An inclined shower with zenith angle θ traversed cos θ1 times more
column depth than a vertical shower. Therefore particles are more attenuated, decreasing the particle densities on the ground. The transverse momenta cause the shower to develop, more or less isotropically, perpendicular to the shower axis (lateral direction). This means that for a vertical shower the footprint is circular around the shower core. But for inclined showers the footprint becomes elliptical. This effect is illustrated in Figure 7, where an inclined shower with zenith angle θ reaches the ground. Particles (right side) at distance r from shower core arrive ∆t later than the core of the shower (front). Particles on the left arrived earlier than the shower core.
Figure 7: Side view of an inclined shower with zenith angle θ at ground level. Particles (right side) at distance r from shower core arrive ∆t later than the core of the shower (front). Particles on the left side of the shower core arrived earlier than the shower core. Figure from [15].
The lateral density of electrons (e− and e+), muons (µ− and µ+) and the sum of both on the
ground (horizontal) of a vertical 1016eV proton shower and a 1016eV proton shower with zenith
angle θ = 45◦ are shown in Figure 8. The density of the vertical shower on the left is dominated
by the electrons. At distances larger than several hundred meters from the shower core the muon density start to become a significant fraction of the sum density and at a distance larger than 1km the muon density starts to overtake the electron density. The sum density of the
θ = 45◦ shower on the right is also dominated by electrons, but here the muon density starts to
overtake the electron density much closer to the shower core, at about 300m. The differences between the densities in the shower core of the vertical and the inclined shower is because of the
1
cos θ times higher column depth. The large decrease of the electron density with respect to the
decrease of the muon density is because electrons are more attenuated than muons. The muon density starts to overtake the electron densities at some (large) distance from the shower core because at these distances the traversed column depth is higher than that to the shower core and electrons are more attenuated than muons. Another reason is the differences in expected lateral distances between electrons and muons. Electrons undergo many interactions in the shower,
which is a random walk process. Therefore their expected lateral distance follows hxi ∝ √h,
with h the vertical distance. Muons travel nearly unattenuated to the ground with few (or none) interactions and their expected lateral distance follows hxi ∝ h.
Figure 8: Horizontal density at sea level of electrons (e−and e+) in red, muons (µ−and µ+) in green and the sum of both as a dashed line of a vertical (θ = 0◦) 1016eV proton initiated shower (left) and a 1016eV proton initiated shower with zenith angle θ = 45◦(right). Figure from [14].
Note that in Figure 8 the photon density is not taken into account. Photons are, however, the
dominant particle at sea level. In Figure 9 the lateral profile of photons, electrons (e− and e+)
and muons (µ−and µ+) at sea level is plotted of a vertical 1015eV proton initiated shower. The
primary energy of this cosmic ray is an order of magnitude lower than that in the left panel of 8.
Note that the densities of the electrons and muons (at all distances) of the vertical 1015eV proton
initiated shower are much lower than that of the vertical 1016eV proton initiated shower. It can
also be seen that for the 1015eV proton initiated shower the muon density starts to overtake the
electron density much closer to the shower core. This is because of the lower Xmax(higher up in
the atmosphere) of the 1015eV shower and the fact that the on average lower-energy electrons
in the 1015eV shower are more attenuated in the atmosphere.
Figure 9: Lateral profile of photons, electrons (e− and e+) and muons (µ− and µ+) at sea level of a vertical (θ = 0◦) 1015eV proton initiated shower. Figure from [15].
The total number of particles on the ground, the shape of the footprint and the lateral densities can be used by cosmic ray experiments to investigate cosmic rays. An example of such a cosmic ray experiment, the HiSPARC experiment, will be described in the next chapter. Also will be explained how this experiment uses the properties of the particles at ground level to investigate cosmic rays and how the experiment relies on simulations of showers.
3
HiSPARC
The HiSPARC Cosmic Ray Experiment aims to investigate cosmic ray air showers and re-construct the primary energy and direction of showers. The scientific purpose of HiSPARC is to study large-scale correlated effects like measuring the Gerasimova-Zatsepin effect. The Gerasimova-Zatsepin effect is when a solar photon interacts with a cosmic ray nucleus, emitting a particle, mostly a proton. A detailed description of this effect can be found in [16]. This effect can be measured by detecting both the shower from the residual cosmic ray nucleus and the shower from the emitted proton. HiSPARC is suitable for measuring this effect because it has a relatively large surface area and (some of) the detectors are close together, much closer than that of the Pierre Auger observatory in Argentina [17], for instance. Because the detectors are closer together lower energy showers can be detected, from which there are many, and the possibility of measuring both the showers is higher. Apart from the scientific questions, another purpose of HiSPARC is to reach out to High School students (and teachers) and introduce them to research in modern physics. Currently, HiSPARC has over 100 detection stations spread out over the roofs of high schools and (scientific) institutions in the Netherlands, Denmark and England, with the vast majority in the Netherlands (about 90%). A HiSPARC detector is composed of a 2cm thick 100cm × 50cm plastic scintillator, a plastic light guide and a photo multiplier tube (PMT) and is illustrated in Figure 10. The detector stations can have two or four detectors and are put in separate ski-boxes on the roofs of the participating high schools and institutions. There are two configurations for four-detector stations: an equilateral triangle with a fourth detector in
the barycenter or a 60◦ diamond configuration. These two configurations are depicted in Figure
11. In Figure 12 the locations of the HiSPARC detector stations in the Netherlands can be seen.
Figure 10: Illustration of HiSPARC detector with scintillator in white, light guide in gray and PMT in black. Figure from [14].
Figure 11: Two configurations of four-detector stations. Triangle station on the left and 60◦diamond station on the right. Figure from [14].
Figure 12: HiSPARC detector station locations in the Netherlands. Figure from [14].
3.1 Charged particles in scintillators
High-energy charged particles in matter can lose energy by collisions with atoms in the medium (ionization and excitation) and by radiation (Bremsstrahlung). In the HiSPARC detector this medium is a plastic scintillator. It consists of polyvinyltoluene, a solid plastic solvent, with a small amount of the fluor anthracene dissolved into it. The electrons in the atoms of the polyvinyltoluene absorb energy from the charged particles, get excited and emit photons in their re-excitation which will excite the fluor (anthracene). π-molecular orbital free valence electrons in the fluor can get excited to singlet or triplet excited states, see Figure 13. After
internal degradation, which is decay without radiation, the molecule de-excites from the S1state
to the ground state or one of the vibrational sublevels of the ground state S0. This process is
called fluorescence and in this process photons are emitted. In the indirect de-excitation of the
lowest excited triplet state (T0) to the ground state also photons are emitted. However, this
process, which is called phosphorescence, is slower than the singlet de-excitation (fluorescence). The emitted scintillation photons will not be absorbed in the scintillator and can travel through the scintillator and the plastic light guide to the PMT (transmission). The PMT signal from
the photoelectrons is converted to a digital value (ADC counts) and this is the detector signal.
Figure 13: π-orbital energy levels, both singlet S and triplet T states of a fluor molecule. In the de-excitation from S1 to the ground state (or its vibrational sublevels S01, S02 etc.) fluorescence takes place and photons are emitted. In the indirect de-excitation from T0 to the ground state (or its vibrational sublevels) phosphorescence takes place and photons are emitted, although slower. There are no photons emitted by internal degradation. Figure from [15].
The detector signal depends on the amount of scintillation photons produced and therefore the amount of energy lost by the particle in the medium by ionization/excitation. The mean energy loss of a particle in a medium, the stopping power, can be calculated with the Bethe-Bloch equation. − dE dx = K 2 Z A 1 β2 ln 2mec2β2γ2Tmax I2 + F − δ − C ! , (1) where −dE dx
is the mean energy loss per thickness of the medium, the stopping power, in
MeV g−1cm2. K is a constant depending on the electron mass and classical electron radius, Z
is the proton number and A is the mass number of the medium. β is the relative velocity, the
velocity of the particle relative to the speed of light, vc, and is close to 1 for most of the incoming
particles. γ = √1
1−β2 is the Lorentz factor and can become very high for particles with velocities
close to the speed of light. Tmax is the maximum kinetic energy transfer in one collision and
depends on the mass of the incoming particle. Therefore, this parameter is different for muons and electrons/positrons, because of the mass difference. I is the mean excitation energy of the medium. F is another parameter that is different for different types of particles and depends on the velocities of the incoming particles. δ is the density correction and is different for different regimes of the value βγ. C is a so-called shell correction and depends on βγ of the particles, the mean excitation energy I of the medium and proton number Z of the medium. In [14]
the stopping power is calculated with the values of Z, A, δ, I and C for polyvinyltoluene, the scintillator material. The stopping power of the scintillator for muons, electrons and positrons is shown in Figure 14 as a function of βγ. However, in this figure only the stopping power due to collisions (ionization) is shown. The contribution of radiation (Bremsstrahlung) to the total stopping power, which becomes important for values βγ > 10, is not shown because the produced photons in this process do not interact with the medium, therefore leave the scintillator, and will not be detected.
Figure 14: Stopping power from ionization for electrons, positrons and muons in the HiSPARC scintillator detector (polyvinyltoluene). Figure from [14].
The mean energy loss −dE
dx can be calculated by the Bethe-Bloch equation, but in reality the
energy loss follows a distribution because of a process called energy straggling. This distribution is called the Landau energy loss distribution and follows:
f (∆) = 1
ξφ(λ), (2)
with ∆ the energy loss in MeV and ξ equal to:
ξ = KZ
2Aβ2x0, (3)
where x0 is the thickness of the medium in g cm−2, which for our 2cm thick scintillator with
density 1.03 g cm−3 equals to 2.06g cm−2. However, this is only true for vertically incoming
particles. Particles incoming under an angle θ have traveled a thickness x0
cos θ in the medium.
φ(λ) is the Landau probability density function with λ given by:
λ = ∆ − ¯∆ ξ − ln ξ Tmax − β2− 1 + CE, (4)
where CEis the Euler-Mascheroni constant and ¯∆ is the mean energy loss (from stopping power).
and muons are calculated for HiSPARC’s 2cm thick polyvinyltoluene-based scintillator. These distributions can be seen in Figure 15. The most probable energy losses, corresponding to the location of peaks in the distributions, are lower than the mean energy losses calculated from the Bethe-Bloch equation. According to [14], the Landau distribution is different for different particles. However, there is disagreement, in [18] only very small differences are found between the Landau distribution of different particles.
Figure 15: Landau energy loss distribution for vertically incoming high-energy electrons, positrons (βγ > 68) and muons (βγ > 4.3) for the HiSPARC scintillator detector. Figure from [14].
As described earlier, particles coming in with zenith angle θ have traveled a larger distance in
the medium, following x0
cos θ. Consequently, the Landau energy loss distribution of the particles
changes with θ. In Figure 16 the energy loss distributions for high-energy muons (βγ ≥ 100) for three different zenith angles is presented. It can be observed that the most probable value of energy loss (location of the peak) increases with θ and also the distributions broaden.
Figure 16: Landau energy loss distribution for high-energy muons (βγ ≥ 100) for zenith angles 0◦, 22.5◦and 45◦. Notice that the most probable value of energy loss (location of peak) increases with θ and also the distributions broaden. Figure from [14].
3.2 Signal
In the following subsections a summary will be given of the HiSPARC signal and how the direc-tion and energy of cosmic rays are reconstructed. A more detailed and more actual descripdirec-tion can be found in [18]. For an event to be stored at least two or more detectors of the same (four-detector) station have to have had a signal within a 1.5µs timeframe. The idea is that if the event is an EAS more than one detector will have detected a particle. Events where only one detector detects a signal are associated with background muons from low-energy cosmic rays
(e.g. E < 1014eV) and have to be rejected. The threshold for a signal is if two detectors of
the same station have a signal larger than 70mV or three detectors detect a signal larger than 30mV.
An example of an event where all four detectors in a four-detector station have a signal can be seen in Figure 17. It is clear that there can be a time difference of several nanoseconds between signals of the same event and that the ADC counts peak (pulse height) is at different heights for the same event in the four different detectors. The peaks of the signal (in ADC counts) is called the pulse height and the integral of the signals (in ADC counts ns) is the pulse integral. The pulse height and pulse integral of a full week of detecting events of a four-detector station are shown in Figure 18 on the left and right, respectively. The peak of the pulse integral is associated with the peak of the Landau distribution. It is the most probable energy loss, and the minimum of the stopping power distribution. Therefore, the value at the peak of the pulse integral is taken as the value for 1 MIP (minimum ionizing particle) instead of in ACD ns. This is, however, not entirely accurate because it depends on the angle of incidence where this MIP peak should exactly be.
Figure 17: Signals (in ADC counts) of an event with a signal in all four detectors of a four-plate station. Figure from [14].
Figure 18: Pulse height and pulse integral of a full week of events for four detectors in a station. Figure from [14].
Due to the uncertainty in quantum efficiency of the PMT and the uncertainty of how many of the produced photons travel to the PMT (transmission efficiency) the real signal (the pulse inte-gral) is a convolution of the Landau energy loss distribution and the detector resolution, which roughly follows a Gaussian distribution. Therefore, the shape of the pulse integral is somewhat broader than the Landau energy loss distribution.
3.3 Direction reconstruction
If at least three or more detectors of a four-detector station detect signals from one shower, the direction of the incoming cosmic ray can be reconstructed from the arrival times of the particles in the detectors. This is also the case if three or more stations from the same cluster, for instance the cluster at the Science Park in Amsterdam, have detected a signal. A detailed description of the direction reconstruction of the HiSPARC experiment can be found in [15].
The direction reconstruction is more accurate for larger arrival time differences between the de-tectors. Then, the timing uncertainties due to the transport in the detector and the shower front shape are relatively smaller. The arrival time differences are larger if the distances between the detectors are larger. Therefore, for the four-detector station in the form an equilateral triangle with a fourth detector in the barycenter, on the left in Figure 11, the direction reconstruction is most accurate if the three corner detectors detect a signal.
In Figure 19 the fraction of simulated events for which all three corner detectors detect at least one charged particle is plotted as a function of core distance for 1PeV proton initiated showers
for zenith angles θ = 0◦, 22.5◦ and 35◦. The detection efficiency is larger for smaller zenith
angles. This is because the number of ground-level particles is larger for showers with smaller zenith due to the shorter distance traversed in the atmosphere. Therefore, the probability to measure particles is larger for these showers. The detection efficiency decreases over distance from the shower core. This is because near the shower core the particle density is higher and it is therefore more probable to measure particles. If the arrival times between two detectors is
larger than the light time between two detectors (∼ 33ns) the detected event is rejected and no direction reconstruction is done because the outcome would be unphysical. Especially events at large distances from the shower core, where the number of detected particles is small and the shape of the shower front is curved more of the detected events will be rejected. This can be seen in Figure 20, where the fraction of detected events is shown for which a direction reconstruction can be done as a function of distance to the shower core for 1PeV initiated proton showers for three zenith angles. The reconstruction efficiency is larger for smaller zenith angles, because of the higher number of ground-level particles for such showers.
Figure 19: Fraction of simulated events for which three HiSPARC detectors at the corners of a 10m-sided equilat-eral triangle detect at least one charged particle as a function of core distance for 1PeV initiated proton showers for three zenith angles. The detection efficiency is larger for smaller zenith angles and decreases over distance from the shower core. Figure from [15].
Figure 20: Fraction of detected events for which a direction reconstruction can be done as a function of the distance to the core for 1PeV initiated proton showers for three zenith angles. The reconstruction efficiency is larger for smaller zenith angles and decreases over distance from the shower core. Figure from [15].
In Figure 21 the precision of the azimuth reconstruction is shown as a function of simulated
azimuth angle. This is done for 1PeV initiated proton showers with zenith angle θ = 22.5◦. It
is found that the azimuth reconstruction is more precise if at least two charged particles have to be detected in all corner detectors instead of one. Note that the uncertainty range in terms of standard deviations is larger than these 50% regions.
The precision of the zenith reconstruction as a function of simulated zenith angle is shown for 1 PeV initiated proton showers in Figure 22. The zenith reconstruction is also more precise for the stronger condition that at least two charged particles have to be detected in the three corner detectors. The precision of the zenith reconstruction increases with zenith angle. This is because for showers with large zenith angles the arrival time difference between the detector is larger. Therefore the relative timing uncertainties due to the transport in the detector and the shower front shape are relatively smaller, improving the precision of the zenith reconstruction. Note that the uncertainty range in terms of standard deviations is larger than these 50% regions.
Figure 21: The precision of the azimuth reconstruction as a function of the simulated azimuth angle for 1PeV initiated proton showers with zenith angle θ = 22.5◦. The white dots are the median values and the grey range contains 25% of the events both above and below the median (50%). For the figure on the right at least two charged particles (instead of one) have to be detected in all corner detectors, which leads to a more precise azimuth reconstruction. Figure from [15].
Figure 22: The precision of the zenith reconstruction as a function of the simulated zenith angle for 1PeV initiated proton showers. The white dots are the median values and the grey range contains 25% of the events both above and below the median (50%). For the figure on the right at least two charged particles (instead of one) have to be detected in all corner detectors, which leads to a more precise zenith reconstruction. The precision of the zenith reconstruction increases with zenith angle. Figure from [15].
3.4 Energy reconstruction
For an Extensive Air Shower (EAS) the lateral profile on the ground, such as the ones in Figure 8 is often approximated by the Nishimura-Kamata-Greisen (NKG) function [19, 20], which depends on the distance from the shower core r and is given by:
ρ(r) = kN r RM s−2 1 + r RM s−4.5 , (5)
with ρ(r) the particle density at a distance r from the shower core, k a constant, N the number
of particles, RM the Moli`ere radius and s a paramater for the age of the shower.
In reality the position of the shower core and the total number of particles N in the shower are not known. These have to be determined by minimizing the following function:
χ2= n X i=1 wi− ρi σi 2 , (6)
with wi the real signal in detector i, ρi the expected signal in detector i from the NKG-function
and σi =√ρi.
If we now write the NKG-function as ρ(r) = N v(r) the χ2 function changes to:
χ2 = n X i=1 (wi− N vi)2 N vi . (7)
The value N for which χ2 is minimized, in other words for which δχδN2 = 0, then equals:
N = v u u tP( wi2 vi) P vi (8)
With this condition for N the χ2 function becomes:
χ2 = 2 v u u tP( w2i vi) P vi X vi− 2Xwi (9)
The real difficulty is then finding the right shower core position for which χ2 is minimized. After
trying lots of different shower core positions the optimal shower core position is found for which
χ2 is minimized. From this the corresponding number of particles, or shower size, N can be
determined.
With the expected shower size N and the expected zenith angle θ from the direction reconstruc-tion the cosmic ray energy can be determined. This is done with the help of a large number of simulated cosmic ray air showers of different energies and with different inclinations. The
dependence of the shower size Ne+µ on cosmic ray energy for three zenith angles 0◦, 45◦ and
60◦ is depicted in Figure 23.
Figure 23: Shower size Ne+µ (with errors) as a function of cosmic ray energy for three zenith angles 0 ◦
, 45◦and 60◦from cosmic ray simulations for proton primaries. Figure from [14].
The error bars describe the shower-to-shower fluctuations between the showers with the same initial conditions. It is important to emphasize that this energy reconstruction can only be done by comparing with lots of simulated events. However, these simulations, and therefore also the link between shower size and the estimated cosmic ray energy for different inclinations, depend on the chosen hadronic interaction model in these simulations. In the next section will be explained what an hadronic interaction model is and why it is so difficult to model high-energy hadronic interactions. Also, three of the main hadronic interaction models will be described in depth, what their differences are and what the advantages and disadvantages of these models are.
4
Hadronic interaction models
Hadrons are particles consisting of quarks and gluons. These quarks and gluons have color charge. When these color charges are exchanged between two separate hadrons a hadronic in-teraction takes place. Inside a hadron there are both valence quarks and an indefinite amount of virtual (sea) quarks. Valence quarks determine the quantum number of the hadron and sea quarks do not. There are two types of hadrons: baryons and mesons. Baryons contain three valence quarks, mesons contain two. Protons and neutrons are examples of baryons. Protons consist of two up quarks and one down quark (uud) and neutrons consist of two down quarks and
one up quark (udd). Pions are examples of mesons. π+, for instance, consist of one up quark
and one down antiquark (u¯d). Hadronic interactions, like proton-proton collisions in hadron
colliders, can be described by quantum chromodynamics (QCD). This theory describes strong interactions, the interactions between quarks and gluons.
However, perturbative QCD can only be used for hard interactions: processes with large mo-mentum transfers. It is known that valence quarks carry most of the momo-mentum of the hadron, often over 50% [21]. In hadronic interactions we can distinguish two kinds of produced particles: leading particles, the particles which contain one or more valence quarks, and the large majority of particles which do not contain any of these valence quarks. Most of the produced particles without valence quarks have a low momentum fraction of the initial hadron. Therefore a lot of these can be produced. On the contrary, leading particles have a flat distribution in momentum. These particles often have large momentum fraction of the initial hadron and only a few of them will be produced. QCD can well describe the production of hard (often leading) particles, but inaccurately describes the production of soft particles. Therefore, to understand multiparticle production, we need to combine QCD with phenomenological models.
At high energies other processes that cannot be described by QCD, for instance, parton satu-ration, take place. At high energies the parton densities, the density of the gluons and quarks in the hadrons, can become so high that the wave functions of the partons start to overlap. This process is called parton saturation [22]. This process decreases the growth of individual partons and therefore decreases the multiplicity, the amount of secondary particles produced, over energy. This process becomes more important with increasing center-of-mass (c.m.) en-ergy. To predict hadronic interactions at extrapolated energies, at center-of-mass energies much higher than that of LHC, a flexible model that describes the processes above for proton-proton, proton-nucleus and nucleus-nucleus interactions is necessary. In the next section three hadronic interaction models are described. These models are called post-LHC models, since they are tuned to the newest LHC data.
4.1 Post-LHC hadronic interaction models
4.1.1 EPOS LHC v3400
The Monte Carlo event generator EPOS [23, 2, 24] was originally used to describe heavy-ion collisions at RHIC, a heavy-ion collider experiment in the United States. EPOS LHC is a much faster version of the original EPOS model, but has more parameters and is therefore less precise than the original one for heavy-ion collisions. For instance, it does not take into account the 3D hydrodynamic flow calculations [24] of the original EPOS model for heavy-ion collisions. But, since these precise calculations are less relevant (less collective hydrodynamic flow) for hadron-hadron and hadron-hadron-nucleus collisions, EPOS LHC is a very good model to describe minimum bias proton-proton (hadron-hadron) and proton-nucleus (hadron-nucleus) collisions. The parameters of the model are tuned to large sets of hadronic interaction data of both hadron, hadron-nucleus and hadron-nucleus-hadron-nucleus collisions. EPOS is based on the Gribov-Regge effective field theory [25, 26, 27, 28]. This describes a hadronic interaction with multigluon diagrams, Pomerons, for both soft and (semi)hard and hard interactions. In a hadronic interaction many interactions are taking place in parallel. An example of such an elementary scattering is depicted in Figure 24 and is called the parton ladder or the cut Pomeron. If the first parton of the projectile/target
(or nucleon in Figure 24) is a sea quark or a gluon and the momentum transfer q2 is lower than
Q20 (few GeV2), the cutoff above which pQCD is applicable, a soft Pomeron is exchanged. These
soft Pomeron emissions lead to a soft preevolution; a nonperturbative parton cascade takes place. This gives rise to the large number of low transverse momentum hadrons. The semi-hard
scatterings with q2 > Q20 can be calculated perturbatively and the parton evolution follows
the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations. For hard interactions the ladder partons eventually form a flux tube, see Figure 24, where fragmentation takes place. Quark-antiquark pairs are formed, which later form the hadrons.
Figure 24: Elementary interaction in the EPOS model. Figure from [24].
get high enough for gluon fusion to occur. This causes elastic and inelastic rescattering of the ladder partons, see resp. Figure 25 and Figure 26, where the ladder parton rescatters of the projectile or target nucleon/parton. The screening due to this elastic scattering leads to a reduction of the total and inelastic cross sections and therefore decreases the particle production. The parameters of elastic scattering depend on nuclear mass number and the c.m. energy of the hadron collision and saturate at very high energies. The inelastic scattering parameters are very important for reproducing collider data (e.g. (charged) particle production and transverse momentum spectra).
Figure 25: Elastic ladder parton rescattering. Figure from [24].
Figure 26: Inelastic ladder parton rescattering. Figure from [24].
Unlike other models, EPOS LHC is on the microcanonical level (parton-level) rather than the canonical level (hadron-level). EPOS LHC calculates all the Pomeron-Pomeron couplings of the multiple scatterings with conservation of energy, flavour, etc. Furthermore, particular for EPOS LHC is the implementation of collective hadronization. In EPOS LHC, a differentiation is made between a core and a corona. In the core the density of the string segments exceeds some density
ρ0. The low density areas are termed the corona. In the core particles will lose momentum due
to multiple scattering. Above some cutoff pcutt particles will lose some of its momentum but live
on as independent particles. Below this cutoff pcutt the particles are absorbed in the core and
these particles will hadronize collectively. For cores heavier than 3 MeV/c2 a collective radial
and longitudinal flow is then calculated. These flows depend on the mass of the core, rather than the c.m. energy. A mass rescaling is done to conserve energy, reducing the particle production in the core.
4.1.2 QGSJET-II-04
QGSJET-II-04 [29, 5] is a Monte Carlo generator for hadronic interactions. QGSJET-II-04 is, like EPOS, based on the Gribov-Regge effective field theory [25, 26]. Not surprising, because the treatment of semihard processes in EPOS is based on QGSJET. However, in EPOS the Gribov-Regge theory is applied to partons (microcanonical), where in QGSJET it is applied to
hadrons (canonical). As a consequence, in QGSJET energy is not conserved at the amplitude level [12]. With the Gribov-Regge theory the soft and semihard scattering contributions can be combined into a general elementary scattering. The soft scatterings with momentum transfer
q2 lower than Q20 are associated with an exchange of a soft Pomeron, which is a representation
of a nonperturbative parton cascade (preevolution). The semihard processes with q2 > Q20 are
associated with semihard Pomeron exchange and the parton cascade evolves perturbatively, ac-cording to the DGLAP equations.
As opposed to EPOS LHC, QGSJET-II-04 does not parametrize the treatment of non-linear effects. QGSJET-II-04 uses the Pomeron-Pomeron coupling from the Gribov-Regge theoretical framework to account for the non-linear effects. The Pomeron-Pomeron coupling represents the interaction between two parton cascades when they overlap in phase space. To calculate the total and elastic cross sections QGSJET-II-04 does a resummation of all the enhanced diagrams, of all the possible Pomeron interactions, in Figure 27, where the cut Pomerons represent the parton cascade and the uncut Pomerons the elastic rescatterings off the cascade, projectile and
target partons. The non-linear effects depend on the triple-Pomeron coupling G3P, which
de-pends on impact parameter and parton density (or c.m. energy).
The low number of parameters in QGSJET do not depend on nuclear mass number (as in EPOS LHC) and only depend on the Pomeron diagrams. As opposed to EPOS LHC, no additional parameters are necessary to go from hadron-hadron to hadron-nucleus or nucleus-nucleus in-teractions. The extrapolation follows from the resummation of the enhanced cut diagrams and extending Gribov-Regge with Glauber theory [27, 28]
One of the disadvantages of the model is that the multi-Pomeron coupling is dominated by soft processes, so there is only saturation for soft processes in this approach. The multi-Pomeron couplings for hard processes are not taken into account by QGSJET-II-04 and therefore there is no screening/saturation for hard processes in the model. As a countermeasure QGSJET-II-04
uses a relatively high Q20 cutoff, so there is more phase space for soft interactions. Another
dis-advantage is that in QGSJET-II-04 energy is not conserved between the multiple scatterings as in EPOS LHC. The advantage of this model over EPOS LHC is the low number of parameters.
4.1.3 SIBYLL 2.3c
SIBYLL 2.3c [30, 3, 31] is the easiest model of the three and is based on the Dual Parton Model [32, 33, 34], Lund MC string fragmentation [35, 36] and the minijet model [37, 38, 39, 40]. In the newer versions of SIBYLL [3, 31] elements of Gribov-Regge theory are implemented. Diffractive interactions are modeled by a two-channel eikonal model similar to the Good-Walker [41] used by QGSJET-II-04. Extrapolation from hadron-hadron interactions to hadron-nucleus interac-tions is done by Glauber scattering [28] (as well as the other models), but for nucleus-nucleus interactions a combination between Glauber theory and superposition model is used: the semisu-perposition model [42].
The basic Dual Parton Model (DPM) is based on quarks and diquarks. Two strings are formed between the quark of the target and the diquark of the projectile, and vice versa. These strings, if the mass of the strings is above some threshold mass, fragment following the Lund string
fragmentation (as do QGSJET-II-04 and EPOS LHC) and q ¯q and qq − ¯q ¯q pairs are created.
At c.m. energies higher than hundreds of GeV this model can not describe very well the increase of multiplicity, mean transverse momentum and the density of central rapidity particles with energy. In SIBYLL 2.3c these effects are associated with low transverse momentum jets from hard interactions and are described by the perturbative minijet model (instead of the DGLAP equations in QGSJET-II-04 and EPOS LHC). To prevent singularities in the cross section from collinear factorization and saturation an energy-dependent cutoff for the transverse momentum,
similar to Q20 in the other models, is introduced in this new version of SIBYLL. The minijet
model is not completely different from the Gribov theory, it follows the same cutting rules as QGSJET-II-04 and EPOS LHC. The number of interactions and minijets is defined as a func-tion of energy s and impact parameter b. A minijet pair is defined as two strings between two
gluons. At the two sides of the strings a q ¯q-pair is created and these form leading particles.
String fragmentation is then done following the Lund model.
Because of the energy-dependent cutoff, SIBYLL 2.3c has a relative large phase space for soft interactions. For these soft interactions some elements of the Gribov-Regge theory are im-plemented into the model to account for multiple scattering. The soft cross section has an energy-dependent part for Pomeron and one for Reggeon exchange, which are parametrized. As the soft interactions cannot be calculated by perturbative QCD the parameters in the soft cross sections are tuned to collider data. Just like the minijets, the soft interactions are treated as fragmenting gluon pairs. However, the transverse momentum distribution of these gluons and the size of this transverse momentum is different from the minijet approach. Note that for multiple interaction at least one valence quark is required.
A disadvantage of the model is that the semisuperposition model for nucleus-nucleus interac-tions is less precise than the full Glauber model as used in EPOS LHC and QGSJET-II-04. A full Glauber model will also make a consistent treatment of diffraction dissociation possible for hadron-nucleus and nucleus-nucleus interactions. Another disadvantage is that the saturation
parameter, in this case the transverse momentum cutoff, only depends on energy and not on the impact parameter like in QGSJET-II-04. Therefore, the parameter is independent of the gluon density, thus has no physical meaning. The advantage of this model is that it is the easiest.
4.1.4 Summary
Model Diffraction by Saturation Effects hh → hA or AA Extra
EPOS LHC Pomerons Parametrized (c.m.
energy and nuclear mass number)
Parametrized Collective flow and
Energy
conserva-tion (microcanoni-cal)
QGSJET-II-04 (cut) Enhanced
Pomeron graphs Pomeron-Pomeron coupling (impact parameter and c.m. energy) Glauber theory SIBYLL 2.3c two-channel eikonal model Energy-dependent cutoff Q0(s) Parametrized (Semisuperposi-tion model)
Table 1: Properties of the three post-LHC hadronic interaction models.
5
Analysis
5.1 First hadronic interaction
The first hadronic interaction is the interaction of the incoming cosmic ray particle with a nucleus in the atmosphere. This is the most energetic interaction in the shower and the properties of this interaction (e.g., the amount of secondaries produced, energy distribution amongst these secondaries) have a big influence on the development of the shower.
To investigate the differences between the model-dependent predictions of properties of produced secondaries in the first hadronic interaction, simulations have been done with the Extensive Air Shower simulation program CORSIKA version 7.6400 [43].
Simulations have been made of showers with proton, alpha and iron primaries with primary
energies in the range of 1012eV to 1020eV in steps of 0.5 log(E) for both EPOS LHC [2], SIBYLL
2.3c [3, 4] and QGSJET-II-04 [5]. The low-energy hadronic interactions (E ≤ 80GeV) are treated by GHEISHA [44] and the electromagnetic component is treated by EGS4 [45]. The first interaction was set at a fixed altitude with a nitrogen nucleus, because this is the most abundant
particle in the atmosphere and to keep similar conditions for all simulated showers. The energy cuts were 300MeV for hadrons and muons and 3MeV for electrons, positrons and photons. The default decay settings were used for all three hadronic interaction models. Also, the default settings for fragmentation of remainder nuclei were used (NFRAGM = 2). For QGSJET-II-04 the nuclei fragment (with evaporation) with transverse momentum of fragments following from experimental data. For SIBYLL 2.3c the fragments have no transverse momenta. For EPOS LHC there is no fragmentation, only evaporation of remainder. For all different primaries, energies and models 1000 events were simulated and the properties of the secondary particles of the first hadronic interactions were stored with the OUTFILE option of the CORSIKA program.
5.1.1 Number of secondaries
An interesting property to investigate of the first hadronic interaction is the number of produced secondaries. In Section 4 is desribed that parton saturation takes place at extrapolated energies
(E ≥ 1017eV), reducing the rise of number of produced secondaries in high-energy hadronic
interactions. It is interesting to compare the model-dependent predicted number of secondaries produced in such high-energy hadronic interactions to investigate what the consequences are of the different parametrizations for these high-energy effects between the models. Also, if one model were to predict a significantly higher number of particles, this is expected to have a large influence on the development of the shower. For instance, if one model predicts a higher number produced charged pions in the first interaction(s), this may result in a higher number of muons at ground level.
In Figure 28 the mean number of secondaries for proton-nitrogen collisions can be seen as a function of primary energy. It can be seen that for larger primary energies QGSJET-II-04 has, on average, the most secondaries, followed by EPOS LHC, followed by SIBYLL 2.3c. It can also be observed that at lower primary energies, at energies well-understood by collider data, there are already differences in the mean number of secondaries between the models. This is because the particle decay of some short-living particles, for instance ρ mesons, is treated differently between the models. This will be further explained later.
For alpha-nitrogen collisions, see Figure 29, QGSJET-II-04 also has the highest mean number of secondaries. But here the difference between the mean number of secondaries of EPOS LHC and SIBYLL 2.3c is smaller, although the multiplicity is still a bit higher for EPOS LHC. For iron-nitrogen collisions SIBYLL 2.3c has, on average, more secondaries than EPOS LHC. This can be seen in Figure 30. Here QGSJET-II-04 has the highest mean number of secondaries as well.
It can be observed that the mean number of secondaries increases with mass number and that the mean number of secondaries of EPOS LHC relative to QGSJET-II-04 and SIBYLL 2.3c decreases with mass number.
12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0
E
prim[log(E[eV])]
0
200
400
600
800
1000
<Secondaries>
EPOS LHC
SIBYLL 2.3c
QGSJET-II-04
Figure 28: Mean number of secondaries as a function of primary energy for proton-nitrogen collisions.
12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0
