Simple hadronic cascade simulations

POSTERS

Simple hadronic cascade simulations

Fernando Sepúlveda; Claudio Dib

Dept. of Physics, Universidad Técnica Federico Santa María, Valparaíso, Chile

ABSTRACT

We obtain results for the average number of muons at sea level in a proton-initiated vertical atmospheric cascade using a simple model of hadronic interactions based on the Hillas splitting algorithm. We study the muon yield at sea level as a function of the proton primary energy, varying the parameters of the interaction model in order to see the behavior of our results. We find that our results are in agreement with experimental data and with those of more sophisticated simulation models for some particular values of the model parameters.

Keywords: Hadronic cascade; Hillas splitting; Muon yield

I. INTRODUCTION

The hadronic component of an air shower is the one that continuously feeds the electromagnetic and muonic components through the decay of neutral and charged mesons respectively. The muon yield at detector level is a key signal to determine the energy and composition of the primary cosmic rays [1] and that is the motivation of this study.

In the simulation of hadronic cascades, the most important ingredient is, of course, the hadronic interaction model, several of which have been developed and are currently used in extensive air shower simulations [2, 3]. We simulate the longitudinal development of proton-initiated vertical cascades within a simple hadronic interaction model based on the Hillas splitting algorithm (HSA) which is computationally very fast [4]. According to this model only pions are produced as secondaries; their energies are generated as we explain in section II. Muons, on the other hand, arise from the decay of charged pions. This work is focused on studying the number of muons at sea level as a function of the proton primary energy. Varying the parameters of the interaction model allows us to see the behavior of the muon yield and we can select the values of the parameters for which our results best agree with known ones, either from experiments or other simulations.

In the next section the interaction model is explained and the average multiplicity in single interactions is studied. In sections III and IV we briefly discuss pion decay and the atmosphere model used. Section V is devoted to the cascade simulation and in section VI results regarding the muon content of the shower are discussed.

II. HADRONIC INTERACTION MODEL

As mentioned above, we study the development of the hadronic component of the shower using a HSA type model for the generation of secondaries in an inelastic hadronic (p-Air or p^±-Air) interaction. Here the secondaries are all pions and each kind is produced with a probability of one third. Hence, the interactions we consider in the simulations are of the following kind:

where X represents any number of p^± and p⁰.

The energies of the secondaries at each interaction, whether the incoming particle is a proton or a charged pion, are generated by the HSA as follows:

1. Split the available energy of the incoming particle in two parts A and B at random.

2. Assign A as the kinetic energy of the leading particle. This particle is of the same kind as the incoming one.

3. Split B at random in J = 2^Nbranches (N a fixed positive integer).

4. Split one of the J branches at random in two parts A' and B'.

5. Assign A' as the kinetic energy of a pion.

6. From B' subtract m_p, split the rest at random in two parts and assign one piece as the kinetic energy of a pion. Call the other piece B' again.

7. Repeat step 6 until the energy remaining (B') is less than some preassigned threshold E_th. This value must be at least as large as m_p.

8. Once E_th has been reached, go back to step 4 with the next branch. Continue until all branches have been used.

This is an energy priority algorithm, which means that energy is conserved at each step, while the multiplicity is a result of the splitting process.

From the description of the algorithm given above, it is clear that the average multiplicity will increase with incoming energy. The algorithm is sensitive to the number of branches J. Tests were run for single interactions with E_th = m_p and J = 1,2,4 and 8 in order to see how the average multiplicity varies as a function of the incoming energy. These results are shown in Fig. 1.

Observing Fig. 1 we can see that, for a fixed incoming energy, the average multiplicity grows with J. This occurs because if the threshold is reached in the splitting process (step 7 in the description above) when J = 1, pion production ceases altogether. On the other hand, if the threshold is reached for one of the branches when J > 1, pion production stops in that particular branch, but we have more branches where pions are created from. Hence the average multiplicity is higher. Since muons come from the decay of charged pions, we expect the muon number at sea level to increase notoriously with growing J. The value J = 4 is of special interest because it is suggested as a good interpolation between two more complex hadronic interaction models [3], GEISHA and QGSJET, in the energy region where neither of these models apply: GEISHA is valid for energies below a few tens of GeV and QGSJET is valid above several hundreds of GeV.

III. PION DECAY

The lifetime of the charged pions is of order 10^-8s, which makes the decay and interaction processes compete in the development of the shower. Hence the simulation must include both processes . The decay channels considered are:

p^± ® m ^± + n _m(_µ)

Since the secondary particles are isotropic in the CM, the energies of the products have flat distributions in their kinematic ranges [5]:

That is, the muon takes an energy with probability distributed uniformly within the ranges given in equation (2) and the neutrino carries the rest. These kinematic ranges are valid in the approximation of a highly relativistic pion and massless neutrinos. In our simulation only muons are tracked down, not the neutrinos.

The lifetime of the neutral pions is of order 10^-16s, so even at the highest energies they decay almost instantaneously. Their decay mode is:

p⁰ ® 2g.

Since the photons are isotropic in the CM, their energies have a flat distribution in the range (0, E_p) for highly relativistic pions.

IV. ISOTHERMAL ATMOSPHERE

Since the natural units for interaction and decay paths are g / cm² and km respectively, we must have a model for the atmosphere in order to transform among these units. In the simulations we consider an isothermal atmosphere. Within this simple model, the ratio between vertical penetration depth (in g / cm²) and density is constant, this constant represents a height scale and is defined as h₀:

In this case, we have an exponential relation between vertical penetration depth and altitude [5]:

with h in km. This last relation, and its inverse, are used in the simulations to convert the units in pion decay distances from km to g / cm².

V. SIMULATION OF THE CASCADES

In this section we discuss in detail the method used for the simulation of the cascades. As was mentioned in section I, we simulate the longitudinal development of vertical cascades initiated by a high energy proton. The processes considered are the following: for protons we only consider inelastic interactions with atmospheric nuclei; for charged pions we consider interactions and decays; finally, for neutral pions only decays are taken into account. All distances are measured in g / cm². The basic idea of the algorithm can be summarized as follows:

1. For each particle with energy above the threshold of interest and altitude above sea level, determine the distance traveled until a decay or an interaction occurs, depending on the case at hand.

2. At that point, generate the corresponding secondaries. Record the depth at which each particle is created along with its energy and kind.

3. Go back to step one for each particle.

A. Protons in the Cascade

In the case of protons we generate the distance traveled between interactions from an exponential probability distribution using the inversion method. We use an energy-dependent interaction length for protons (in g / cm²) as suggested in ref. [6]:

where E is the energy of the incoming proton. Once the interaction point is known, the secondaries are generated using the HSA as described in section II. When the energy of the proton falls below 10 GeV, no more interactions are considered and the proton is only propagated until it reaches sea level.

B. Charged Pions in the Cascade

In the case of charged pions, we must consider both interaction and decay processes. This is done by sampling two distances, one for each process, and selecting the shortest one. If the selected process is interaction, the secondaries are generated at that point using the HSA as described in section II. Instead, if the selected process is decay, a muon and a neutrino are created using a flat distribution as discussed in section III. For interactions, the distance traveled is generated in the same way as for protons, but with an interaction length (in g / cm²) [6]:

where E is the energy of the incoming pion. In an isothermal atmosphere, the decay length of charged pions, in g / cm², is given by [5]:

where e _p = 115 GeV, E is the energy and X is the vertical penetration depth of the pion. Since d_p depends on X, the distance traveled until decay is not so easily sampled. The natural units for the decay distance are kilometers. Hence, we must generate the decay distance in these units. This is achieved using the inversion method with an exponential probability distribution of mean path l_d = g c t, where g and t are the lorentz factor and the mean life of the pion, respectively, together with unit conversion between g / cm² and km and viceversa using the exponential relation of equation (5).

The algorithm for the generation of the pion decay distance can be described as follows:

1. Take the initial position of the pion X₁, which is in units of g/cm², and transform it to a height h₁ in km.

2. Calculate l_d = g c t for the pion in question.

3. Generate the decay distance l_d from an exponential distribution of mean path l_d.

4. With the decay distance l_d, calculate the height h₂ at which the pion decays: h₂ = h₁-l_d.

5. Transform h₂ to g/cm² to obtain the final position of the pion X₂.

6. The decay distance in g/cm² is thus X₂-X₁.

The probability for a pion to decay as it traverses a small DX of atmosphere is DX/d_p. As can be seen from equation (8), this probability is negligible if E e_p. On the other hand, decay is the dominating process if E e_p. In the simulations we separate three energy regions for pions:

1. When E > 1 TeV only interactions are considered. The distance traveled is generated from an exponential probability distribution with mean free path l_p and secondaries are generated using the HSA.

2. When E < 10 GeV only decays are considered. The distance traveled until decay is generated as discussed above. The secondaries are created as discussed in section III.

3. When 10 GeV < E < 1 TeV, two distances are sampled, one for each process just like above, and the process corresponding to the shortest distance is the one that occurs.

C. Neutral Pions in the Cascade

In the case of neutral pions, due to their short lifetime, only decays are considered. The decay is simulated right where the neutral pion is created. The photons produced share the energy of the parent p⁰ at random, but here we do not treat this part of the process, as it feeds into the electromagnetic component of the shower, which does not subsequently affect the hadronic development.

D. Muons in the Cascade

Muons above threshold are propagated from their point of creation down to sea level without considering any interaction or decay processes. Since we use thresholds above 10 GeV, disregarding decays does not affect the results.

VI. RESULTS AND CONCLUSIONS

We concentrate on the muon content of the shower at sea level. As mentioned in section II the key parameter in the simulation is the number of branches J in the HSA. We study the muon number at sea level as a function of primary energy and its behavior as J is varied. The simulations were carried out for primary protons with energies between 10¹⁴ and 10¹⁶ eV using J = 1, 2, 4 and 8 for each primary energy. The results are shown in Fig. 2.

Irrespective of the value of J used, we find that the average number of muons with energies above 10 and 30 GeV at sea level as a function of the primary energy is very well fitted by a power law:

where all correlation coefficients satisfy R² > 0.999 and 0 < b < 1. The values of b and C depend on the muon energy threshold and on the particular J used in the HSA. These values are shown in table I. The general result expressed in equation (9) is expected from what is known from more complex models and from experimental results [2, 5].

In table I we can see that for each value of J the exponent b is larger for the 10 GeV threshold and grows smoothly as J goes from 1 to 8. As expected, Fig. 2 shows a marked increase in the number of muons at sea level as J goes from 1 to 8 for every primary energy. However, each time J is doubled this increase is progressively slower. For example, take the case of the 10 GeV threshold and E₀ = 100 TeV. When J goes from 1 to 2, N_µ grows in a factor of 1.4; when J goes from 2 to 4, N_µ grows in a factor of 1.17; and when J goes from 4 to 8, N_µ grows only in a factor of 1.12. This behavior occurs because as the average multiplicity in the HSA grows with J, the average energy per pion decreases, which means the muons produced in the decays are less energetic.

The values of C and b for J = 4 and 8 are in very good agreement with those obtained by Stanev [2] with the SYBILL model of hadronic interactions [7] and by the simple, deterministic model of Matthews [8]. Furthermore, these results are in very good agreement with the experimental results reported by the Akeno group [9]. On the other hand; the results obtained with J = 1 and 2 have values for b that are too low, especially for J = 1. This is due to the extremely low multiplicity given by the HSA for these values of J, as can be seen in Fig. 1.

This simple model for the simulation can be improved in several ways. We can consider a more realistic model for the atmosphere. Also, we can improve on the hadronic interaction model modifying the HSA algorithm by setting the number of branches J to be energy-dependent and fit hadronic data more precisely.

It is natural to conclude that the parameters in equation (9) are sensitive to the hadronic interaction model. The advantage of this approach lies in the computational speed of the HSA and in its flexibility.

Acknowledgments

This study has been partially funded by Fondecyt, Chile grant no. 1030254.

[1] A. Haungs, H. Rebel, and M. Roth, Rep. Prog. Phys. 66, 1145 (2003).

[2] Todor Stanev, High Energy Cosmic Rays, Springer-Praxis 2004.

[3] J. Knapp, D. Heck, S.J. Sciutto, M.T. Dova, and M. Risse, Astropart. Phys. 19, 77 (2003).

[4] A. Hillas, Proc. 17th Int. Cosm. Ray Conf. 8, 193 (1981).

[5] T. K. Gaisser, Cosmic Rays and Particle Physics, Cambridge University press 1990.

[6] l is taken as constant for incoming E < 0.1 TeV, as suggested in ref. [5], and decreasing logarithmically for larger E. This logarithmic term is fitted to the values of Table 5.1 in ref. [5]. These values agree with the p-p and p-p collisions listed in S. Eidelman et al., Phys. Lett. B592, 1 (2004) and converted to p-Air and p-Air collisions, respectively, according to B. Kopeliovich et al., Phys. Rev. D 39, 769 (1989).

[7] R.S. Fletcher, T.K. Gaisser, P. Lipari, and T. Stanev, Phys. Rev. D 50, 5710 (1994).

[8] J. Matthews, Proc. 27th Int. Cosm. Ray Conf. 1, 261 (2001).

[9] M. Nagano et al., J. Phys. G 10, 1295 (1984).

Received on 20 December, 2006

Publication Dates

• Publication in this collection
  13 Aug 2007
• Date of issue
  July 2007

History

• Received
  20 Dec 2006