Small shower CORSIKA simulations

The first and undoubtedly most important parameter of a shower is its size, which is understood here as the number of particles at the observation level, as it is typical for surface shower arrays. In experiments analyzing the light (Cherenkov or fluorescent light) produced when a shower passes through the atmosphere, the shower size is defined as the number of charged particles in the shower maximum. CORSIKA allows us to precisely count the number of tracked electrons (and positrons), as well as positive and negative muons, that eventually reach the observation level. Certainly, all thinning options in the program have to be switched off. Repeating simulations continuously for the same particle with the same energy and the same angle of arrival, we will obtain a different result each time. Fluctuations in the size of the shower are an intrinsic property of the simulations, and they result from the probabilistic nature of multi-particle production processes. Figure 1 presents examples of electron and muon shower size spreads for cases of showers initiated by protons with energies of 10 $^{13}$ eV and 10 $^{15}$ eV. The lines correspond to the Gaussian distribution (actually Log-normal) fitted to the iluustrated histograms. As can be observed, in all cases, the Log-normal distribution optimally describes the scattering of the sizes.

Zoom In Zoom Out Reset image size

For showers smaller than those illustrated in Fig. 1, when the sizes are approximately a few particles, the Log-normal distribution is obviously no longer the best. At small number of particles, we should expect an important correlation between them triggered by their common origin from subcascades, which developed accidentally just above the level of observation. In addition, in such cases, the Poisson-type statistical fluctuations are superimposed on fluctuations in their number caused by the probabilistic nature of the phenomenon itself.

The discrete nature of the size variable means that we expect a number of shower cases with zero charged particles. The analysis of the simulation outputs indicates that there are more of such cases than would originate from the Poisson nature of the process. The fraction of "empty" showers increases rapidly with a decrease in the energy of the primary particles. This relationship is important if you want to integrate the observed fluxes of particles appropriately. Obviously, this truncates the flux of primary cosmic ray particles on the low energy side. Fig. 2 shows this cut for vertical showers. It is represented as black dots for the results of simulations with the CORSIKA program for primary protons. As can be observed, the truncation for all charged particle size starts to work below the energy of 10 $^{11}$ eV, and at the energy of 10 $^{10}$ eV only approximately 1 shower in 100 contains at least one charged particle (Fig. 2c). The same figure also presents results obtained in the case when the primary particle is an iron nucleus. As expected, the truncation in this case starts an order of magnitude earlier (below 10 $^{12}$ eV) and is significantly more abrupt. We will discuss effects connected to the mass of the primary particle below.

Zoom In Zoom Out Reset image size

A comparison of the truncation for the electron (Fig. 2a)) and muon (Fig. 2b)) sizes indicates that for very low energies, muons are the particles that manage to reach the ground, and they determine the counts in the individual detectors. In the following, we will analyze this effect quantitatively.

As already mentioned, the size of the shower initiated by a primary particle of a given energy is the most important (for integration of surface particle flux) parameter, which the CORSIKA program provides. We can define the size of simulated showers separately as electron and muon sizes, as the number of electrons or muons at the observation level. Both of these quantities are important for our purposes, and both are further analyzed in parallel.

Figure 3 presents the average values obtained from simulations for primary protons (filled symbols), where the proper energy scale is represented on the bottom axis and the values on the left axis, and for the iron nuclei empty symbols, on the top and right scales. First, it is important to note that, practically, in the entire energy range presented from 10 $^{11}$ eV to 10 $^{15}$ eV, the dependence definitely satisfies the power-law with different indices for muon and electron sizes.

Zoom In Zoom Out Reset image size

First, the simplest assumption concerning the relationship between quantities describing showers initiated by protons and by the complex atomic nuclei is the assumption of simple superposition. According to this concept, a nucleus is considered a set of single nucleons, which behave like protons in their interaction with the nuclei of atoms of the Earth's atmosphere. Consequently, a shower initiated by, for example, an iron nucleus is the same as 56 proton showers. This assumption is relatively natural and correct to a large extent, as can be observed in Fig. 3, where the average sizes of proton showers are compared with the sizes of iron showers divided by 56 and shifted on the energy scale to the same energy per nucleon.

The superposition model was also adopted to compare the number of showers with no particle at the observation level in the proton and iron events illustrated in Fig. 2. The possibility that none of the $56$ showers, which are parts of a shower initiated by an iron nucleus, will have any particle at the observation level, i.e., electrons Fig. 2a), muons Fig. 2b), or none of them Fig. 2c), corresponds to a $56$-fold observation of an "empty" proton shower with the same energy per nucleon. The result of such an assumption is presented in Fig. 2 by empty squares connected with a dotted line.

Both of these observations (of the CORSIKA results) indicate that the superposition assumption is correct; however, this is not entirely and exactly true. Another important characteristic of EAS, which is indispensable for carrying out correct calculations of particle fluxes in small showers, is the size of fluctuations of shower sizes at small energies of particles that are initiating them. Examples of such fluctuations are presented in Fig. 1. As we have already demonstrated, in the first moments, the average values of these distributions agree with the superposition assumption; however, in the case of the second moments, this phenomenon is not as significant.

Figure 4 presents the dispersion of the logarithm of the electron (a) and muon (b) size distributions as a function of the energies of the protons and iron nuclei initiating the showers. The energy scales, abscissa, at the top and bottom are for the iron and proton showers, respectively, and the ordinates are depicted on the left for proton and on the right for iron showers (see Fig. 3). Energy scales for iron showers correspond to the same energy per nucleon as respective scales for proton induced showers. The ordinate for iron showers (right) is scaled down by $\sqrt {56}$. In the simple superposition figure, the dispersion of $56$ independent proton showers would correspond to the expected dispersion for iron induced showers.

Zoom In Zoom Out Reset image size

As can be observed, the points from simulation calculations with the CORSIKA program for proton initiated showers do not overlap with the correspondingly shifted values for iron showers. We ignore the moment of discrepancies for very small showers with sizes below approximately 10 particles (for irons below $\sim 5 \times 10^{13}$, for protons below $\sim 10^{12}$). Outside this area, the fluctuations for iron induced showers are definitely wider, by approximately 50% for electron size, and slightly smaller by about 30% for muon size from expectation corresponding to proton showers at corresponding energies. We will explain this discrepancy below. The lines in Fig. 4 represent our proposed solution, the results of the "fast small shower generator."

It has been known for many years that the transverse distributions of particles in extensive air showers are well described by a simple formula proposed by Greisen [19]. Its validity was confirmed by the theoretical considerations and numerical calculations with respect to electromagnetic cascades by Kamata and Nishimura [20]; hence, its commonly accepted name: Nishimura-Kamata -Greisen (NKG) function

$$\rho_{e/\mu}(r) = {N_{e/\mu} \over 2 \pi r_0^2}\: {\Gamma(4.5-s) \over \Gamma(s)\:\Gamma(4.5-2s)}\left({r\over r_0}\right)^{s-2} \left({1+{r\over r_0}}\right)^{s-4.5}\; \; \; , \tag{ 1 }$$

where $N_{e/\mu}$ is the electron/muon shower size, $r_0$ is a radial scale parameter, which is called the Molière unit in the case of the electromagnetic cascade theory, and is equal to approximately 2 radiation length units above the observation level. In addition, $s$ represents the age parameter.

Figure 5 presents transverse distributions of electrons and muons in small vertical showers initiated by protons and iron nuclei with relatively high energies (10 $^{13}$ and 10 $^{15}$ eV). As can be observed, the NKG function satisfactorily describes these distributions. It is evident that muons propagate to larger distances from the shower axis. The 'age parameter' $s$ characterizing the slope of the distributions, as well as the characteristic scale parameter $r_0$, is larger for muons than for electrons. The description of the lateral spread with the NKG function becomes problematic for very small showers. If the number of particles is small, one cannot refer to their distribution in a single event. Only the mean distribution can make some sense; however, even when applied to a single shower, it makes little sense because of the significant distortion via event-by-event fluctuations. Figure 6 illustrates this situation. It shows the averaged transverse distributions of electrons and muons in showers initiated by protons with energies $10^{10}$ eV and $10^{11}$ eV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For showers containing an average of approximately 1 particle, which corresponds to a (mean) central density of 10 $^{-6}$ per m $^2$ (100 GeV proton – Fig. 6b)), the NKG function does not yet perform poorly in describing the mean distributions of electrons and muons; however, at densities of 10 $^{-8}$ per m $^2$, which is the case as observed in Fig. 3 for 10 GeV protons and iron nuclei with a total energy of 1 TeV, the distributions can be considered to be almost uniform, and shower particles can appear almost everywhere, up to a distance of several hundred meters from the shower axis.

For the purposes of this study, the lack of a statistically verified mean distribution of particles in very small showers is not particularly important, because the integrals of these distributions we eventually aim for will be normalized by the total number of particles in such showers, which will be approximately 0.01, as shown in Fig. 3.

First, we verified a simple superposition model that allows us to discard the mass of the primary particle as an independent variable from our descriptions. As we have already demonstrated, this model works well for $Z$ and $\langle N \rangle$ variables. We have observed a significant discrepancy by comparing the dispersion (logarithmic) of the total number of particles for protons and iron nuclei. The observed value is greater than expected. The expectations obtained via scaling by the square root of the mass number $\sqrt{56}$ spread of the iron shower's size indicates the existence of a correlation between the "sub-shower components," provided some kind of superposition is assumed.

The simplest introduction of correlations among $A$ independent nucleon (proton) sub-showers with the same energy per nucleon as the complex nucleus they are supposed to mimic is to assume that some of them have identical actual sizes. This decreases the number of independent compounds, thus increasing dispersion, the spread of additive shower parameters.

It is well known that the main source of shower fluctuations is the height of the first interaction, which is related to the cross section of the interaction; hence, it is practical and theoretically justified to assume that the number of identical sub-shower components is proportional to the number of wounded nucleons in the interaction of the cosmic ray nucleus with the atomic nucleus of the atmosphere. Certainly, we are only referring to the wounded nucleon of the beam nucleus $A$.

In the case of iron and the atmospheric nucleus of nitrogen, the average number of wounded nucleons at the energies of interest in iron is approximately 12. Because we intend to use our superposition model also for other nuclei, we have approximately parameterized the number of wounded nucleons as a function of mass number, as presented in Fig. 7.

Zoom In Zoom Out Reset image size

The calculations performed indicate that the function presented in Fig. 7 does not provide a sufficient increase in the dispersion of the distribution (logarithmic) of the number of electrons and muons in iron showers. Because the assumption that the number of identical sub-showers depends on the number of wounded nucleons seems to be well-founded, in the next step, we assumed that it is not exactly equal but proportional to the number of wounded nucleons leaving the proportionality factor as a free parameter to be adjusted to the CORSIKA results. The actual integer number of identical sub-showers in each nucleus induced shower fluctuates according to a binomial distribution.

The results of the modified superposition model are shown in Fig. 4 for electrons and muons. As can be observed, in the region where the shower sizes are sufficiently large, we have obtained a sufficient agreement with CORSIKA simulations. As already mentioned, for very small showers, the integer number of particles in single showers and the combination of the physical spread with the effect from the Poisson distribution play a dominant role.

The modified superposition method provides an opportunity to determine the average size (electrons and muons) of the shower initiated by any nucleus from H (protons) to Fe (iron), which dominate the cosmic ray flux of energies in the range of our interest from 10 $^{10}$ eV to 10 $^{16}$ eV. The dispersion determined allows us to consider the random spread of the shower size at a fixed energy of the primary particle. Some examples are illustrated in Fig. 8.

Zoom In Zoom Out Reset image size

There is a weak positive correlation between moun and electron sizes of CORSIKA showers, which is represented in Fig. 8 by solid lines and introduced in our small shower generator.

The dependence of the average number and dispersion of electrons and muons on the energy of the primary proton, shown in the figures above, was obtained for vertical showers. Simulations for inclined showers allowed to parametrize all these dependencies as functions of the arriving particle zenith angle.

For the particular size, the radial distribution of shower particles (for electrons and muons separately, with different scale paramerer $r_0$) in the form of Eq. (1) can be determined using the age parameter $s$ adjusted to small shower CORSIKA results, considering the geometry of inclined showers.

These findings are already sufficient to formulate a simple algorithm generating small EAS that will resemble the EAS generated by the CORSIKA program. Examples of showers generated in this manner are presented in Fig. 9.

Zoom In Zoom Out Reset image size

Eventually, the procedure is developed to generate the particular density of shower particles (electrons and muons) in the shower initiated by any cosmic ray primary particle of any energy, pointing to any point on the observation plane at any shower inclination. This allows to compare the results of our small shower generator, not only with the CORSIKA showers that it is based on but also with experimental results.

The shower particles density spectrum was measured since at least the middle of the last century. The form of the spectrum determined agreed with a simple power law formula, for example the one measured by Cocconi, Loverdo, and Tongiorgi in 1946 for densities from approximately 10 to 1000 particles per meter squared [21, 22]

$$N(x) = 700 \times x^{-1.47} , \tag{ 2 }$$

where $x$ represents the particle density (m $^{-2}$) and $N$ is the rate of events of densities higher than $x$ (hour $^{-1}$). Another example is the result obtained by Broadbent et al. in 1950: $\left( 620 \times x^{-1.425}\right)$ for slightly lower densities [23], and the one by Norman in 1956: $\left(540 \times x^{-1.39}\right)$ (for $x \lt 500$ per m $^{2}$) [24], or Greisen in 1960, who presents a similar expression for the density spectrum in the range $1 \lt x \lt 10^4$ per m $^{2}$ with the index of $(-1.3)$ [25].

With our fast small shower generator, which not only reproduces the average shower characteristics determined for a fixed energy of the primary particle (including its mass and angle of arrival) but also considers respective dispersions, the multidimensional Monte Carlo integration can be carried out to obtain the shower particle density observed with a single small detector. The general formula is quite trivial:

$$\begin{aligned}[b] f(\varrho_{e/\mu}) = & \sum_A \: \int \limits_{10^9}^{10^{16}} {\rm d}E \:\: \Phi_A(E)\: \int\limits_0^{90}{\rm d}\phi\: \int {\rm d}{N_{e/\mu}}\:\\ & \times p\left( N_{e/\mu},\langle N_{e/\mu}\rangle\right) (E,A,\phi) \\ & \times \:\int\limits_{\rm -2km}^{\rm 2km}{\rm d}x \:{\rm d}y\: 2\pi\:\sin(\phi)\:\cos(\phi)\:\varrho(r) , \end{aligned} \tag{ 3 }$$

where $A$ represents the mass number of the primary cosmic ray nucleus, $\Phi_A(E)$ is the cosmic ray energy spectrum of the component of mass $A$, $\langle N_{e/\mu}\rangle (E,A,\phi)$ indicates the average electron/muon size of the EAS initiated by the particle of mass $A$, energy $E$ originates from the direction specified by the zenith angle $\phi$, and $p\left( N_{e/\mu},\langle N_{e/\mu}\rangle \right)$ is the probability density distribution that the actual electron/muon shower size is $N_{e/\mu}$ when the expected value is $\langle N_{e/\mu}\rangle$.

The charged particle density is $\varrho = \Big[ \rho_e \left( r(x, y, \phi)\right)+$ $\rho_\mu \left(r(x, y, \phi) \right) \Big]$, where $\rho_{e/\mu}$ represents the electron/ muon density described by the NKG-type formula in Eq. (1) with adjusted parameters $r_0$ and $s$ , which again depend on primary particle energy, mass, and angle. The radius $r$ is the distance to the shower axis, whereas $x$ and $y$ specify the position of the shower core in the observation plane.

The upper limit of the primary particle energy spectrum depends on the values of densities we attempt to study. This study is concerned with small densities and small showers; hence, if we limit ourselves to densities not exceeding a few hundred per m $^2$, the 10 $^{16}$-eV limit is quite sufficient, as will be verified in this paper. The $2$ km limit on the $r$ integration is again considered with a surplus that is safe for the problem we are discussing.

Summing over the mass spectrum of cosmic ray particles, for practical reasons, as it is conventionally practiced, was changed to summing over a few groups of particles with similar masses: protons ( $A = 1$), Helium ( $A = 4$), CNO ( $A = 14$), Medium ( $A = 28$), and Iron group ( $A = 56$). Uncertainties associated with such a simplification are negligible, considering the statistical scattering and smearing of EAS parameters by the process of particle transport in the atmosphere. In the energy region of interest, a large amount of data exist on mass composition obtained from direct satellite, balloon, and extensive air shower experiments at high altitudes and sea levels. However, the question of the distribution of cosmic ray nuclei masses remains unresolved. There are essentially two models competing in the literature, both with solid theoretical justification, and both describing several of the observed shower parameters. The first one, called heavy dominant (HD), may be related to a supernova acceleration and rigidity dependent propagation model, such that the proton component is assumed to bend at an energy of approximately 10 $^{14}$eV. The second, a proton dominant model (PD), assumes a proton dominant chemical composition over the entire knee energy region of our interest. The fractions of mass groups in two models at the total energy of 10 $^{14}$eV are: protons - 26 (34), Helium - 14 (17), CNO - 19 (19), Medium - 18 (16), Iron - 23 (14), where the first numbers are for the HD model, while those in brackets are for the PD model, respectively [26]. The absolute flux of each composition agrees with that obtained by measurements in the energy region at approximately 10 $^{13}$ – 10 $^{14}$ eV [27].

Certainly, if the total observed flux of primary cosmic radiation solely comprises iron nuclei, the flux of muons, as well as electrons, would be significantly larger than that of a purely proton composition; however, the differences in particle intensities at sea level for mixed mass spectrum compositions, PD and HD, are small, the calculated total muon flux at the sea level is 97 m $^{-2}$ s $^{-1}$ , and the electron flux is 29 m $^{-2}$ s $^{-1}$ for both models. The total flux of primary cosmic rays for both compositions was the same [28].

The integration results in Eq. (3) are depicted in Fig. 10, in comparison with measured results presented above [22-24]. As can be observed, the agreement exhibited is very good.

Zoom In Zoom Out Reset image size

Studies on the particle density spectra allow, to some extent, conclusions to be drawn on the energy spectrum of cosmic ray particles at the top of the atmosphere. In addition, averaging over the position of the shower axis ( $x$, $y$) does not definitely cancel the complex proportionality (averaged) observed on a single detector particle density with the mean energy of primary particle corresponding to this density. Obviously, the simple proportionality between the density $\rho$ and the energy $E$ is smeared out owing to the propagation of the shower in the atmosphere; however, the study on the slope of the density spectrum allows general conclusions to be drawn on the primary energy spectrum. In accordance with this idea, in 1957, Zawadzki proposed, for the first time, the observation of an abrupt change in the density spectrum [29], eventually confirmed later in several EAS experiments and conventionally known today as the "knee."

By integrating the spectra illustrated in Fig. 11 , we can obtain the fractions of corresponding observations and, for example, the rate of registering single muons.

Zoom In Zoom Out Reset image size

If we measure the number of particles observed on a 1 m² detector, each observed value corresponds to a different distribution of the primary particle energy. The results of the calculation are presented in Fig. 11, where the energy spectra of the primary particles leading to observations of 1, 2, 3, 4, and 5 electrons (solid lines) and muons (dashed lines) are illustrated. These results lead to some important conclusions:

$-$ for single muon or electron registrations, the lower energy limit of the primary particle energy is below 10 $^{11}$ eV

$-$ the rate of registration of single muons is several times higher than that of single electrons.

$-$ simultaneous registrations of several electrons (on a 1 m² detector) effectively start from the primary particle energy of 10 $^{12}$ eV for double registrations, up to 10 $^{14}$ eV for the simultaneous registrations of 5 particles.

$-$ cases of simultaneous registrations of more than one muon are significantly rarer than those of more than one electron, and the primary energy required is substantially higher.

By integrating the spectra shown in Fig. 11 , we can obtain the fractions of corresponding observations and, accordingly, the rate of registering single muons.

As already mentioned in Sec. I, the main application of the small shower generator is to assist in the interpretation of data from small shower arrays, either for educational purposes or for its applications in integrated networks of such local stations.

For example, let us assume that such stations would consist of four identical detectors located not far from each other (at a distance of 5 meters). We assume that the detector is ideal and generates no noise, which allows the interpretation of events when only one of the detectors has registered anything, usually one particle. The practical importance of cases where one detector registers several particles according to the results shown in Fig. 11 is negligible. Although these events are very rare, such cases can be analysed via the simulations of the small shower generator.

Let us further assume that every registration of a particle by any detector will trigger an event and that the binary state (hit/no hit) of each detector will be stored.

The small shower generator will help in answering the question on the energy of the primary particle required, or more precisely, the question on the energy distribution that should be associated with a given type of coincidence. Examples of such results are presented in Fig. 12.

Zoom In Zoom Out Reset image size

With the spectra illustrated in Fig. 12 , we can determine the rate of particular coincidence. The obtained results are presented in Table 1.

As mentioned in Sec. I, the size of detectors in local, school arrays is a very important parameter in their design. It seems that the 1 m $\times$1 m size is impractical (and too expensive). Among the planned school experiments, detector layouts of 0.5 m² are being investigated; however, they are not significantly smaller than those of 1 m $\times$1 m. Another approach is to adopt smaller detector sizes of 10 cm $\times$20 cm that can also be used for other measurements/experiments [10].

The effect of the dominance of the electron component for higher rank coincidences, which we have shown in Fig. 12b), is substantially important for the planning of small, school shower arrays and their detailed location. Detectors placed under concrete roofs in physics labs will record shower events significantly less frequently. To obtain a four-fold coincidence once an hour, detectors with a size of 10 cm $\times$20 cm should be positioned, such that they are not shielded by anything from the soft, electron component of the showers. The same detectors shielded and only exposed to the muon, the hard component, will give one quadruple coincidence per year.