Study of longitudinal development of air showers in the knee energy range

Cosmic ray are energetic particles from space. The most important features of cosmic ray spectrum are the so-called “knee” around 3 × $10^{15}$ eV and the “ankle” around 5 × $10^{18}$ eV [1]. The origin of the “knee” of cosmic ray energy spectra is still under debate, however, it is widely believed to be a strong constraint for acceleration and propagation models [2, 3]. One method to understand the origin of the “knee” of cosmic ray spectrum is directly measuring the spectrum and isotropy of cosmic ray, especially the energy spectra and the “knee” of individual species [2, 3, 4].
For particles around the “knee” energy, currently they are only detected by ground-based experiments through the extensive air shower (denoted as EAS hereafter) of the primary particles. There are several secondary particles produced in EAS, most of which are gamma rays, electrons/positrons, muons and neutrons/anti-neutrons. Meanwhile, Cherenkov light and fluorescence light are produced during the propagation of charged secondary particles (mostly electrons and positrons) in the atmosphere. The density of secondary particles at the observation level versus the perpendicular distance to the shower axis (namely, the lateral distribution) has been well studied in the literature [5, 6, 7, 8]. This was used to reconstruct the core position, energy and the type information of primary particles. But, more so, it is an important way for detecting the cosmic ray by EAS. The number of secondary particles versus the slant depth of atmosphere traversed (namely the longitudinal development) also contains energy and type information of the primary particles, which is another important way for detecting cosmic ray by EAS.
The longitudinal development of EAS can be studied through Cherenkov light, fluorescence light and muon, because they interact very rarely with the atmosphere during propagation, while the other secondary components interact with the atmosphere very frequently before reaching the observation level. The longitudinal developments of Cherenkov light and fluorescence light were normally measured by telescope, which converged the Cherenkov/fluorescence light on the focal plane to have an image of the EAS. Since Cherenkov light and fluorescence light were produced mostly by electrons along the trajectory, the longitudinal developments of Cherenkov light and fluorescence light were regarded as a measurement of the longitudinal development of electrons. The longitudinal development of muon can be reconstructed from the arrival time of muon at the observation level as documented in the literature [9, 10, 11]
Basically, ground-based experiments detect the primary particles by measuring the longitudinal development and lateral distribution of EAS. For example, HESS, MAGIC and VERITAS experiments detect gamma rays by measuring the longitudinal development of Cherenkov light, the reconstructed angle resolution within $0.1^{\circ}$, and the energy resolution around 15% [12, 13, 14]. WFCTA of LHAASO experiment detect cosmic ray nuclei also by measuring the image of Cherenkov light, the reconstructed energy resolution for proton around 10—15% in the “knee” energy range, the length, width and angular distance variable (the angular distance between the center of the image and the direction of the primary particles) from the image, were used to identify the proton from other nuclei [15]. Pierre Auger experiment, Telescope Array experiment and HiRes experiment constructed several telescopes to measure the longitudinal development of fluorescence light from cosmic ray nuclei around and above the “ankle” energy, measuring the all-particle spectrum and composition of cosmic ray  [16, 17, 18]. The experiments that measure the lateral distribution of the secondary particles from EAS, e.g., KASCADE experiment, HAWC experiment and many others were also reported [19, 20, 21, 22]. In the last decades, there had been more hybrid experiments that measured both the longitudinal development and lateral distribution, e.g. the LHAASO experiment, the Pierre Auger experiment, the Telescope Array experiment and others, measuring the EAS comprehensively [23, 16, 17]. However, some details about them remain unclear. It is therefore important to understand the advantages and disadvantages of the two detection methods, and combine them to have a better performance. In this paper, the differences of longitudinal developments between electron and Cherenkov light were investigated. The energy and composition measurements based on longitudinal development of EAS were compared with that based on lateral distribution. The uncertainties of energy and composition measurements based on longitudinal development and lateral distribution were both composed of two parts, with the first from the shower-to-shower fluctuations, and the other from the measurement uncertainties due to the detector. The measurement uncertainties due to detector were not studied in this paper, since it is different for different detectors, and can be reduced with better designed detectors. So the performances derived in this paper can be regarded as an upper limit for the variables studied, and they can still provide important information for the selection of detector type, detector design and data analysis for physical measurements.The paper is organized as below, the simulation procedures and parameters are introduced in section 2. The properties of longitudinal development, the differences between electron and Cherenkov light are presented in section 3. The performances of energy measurement and composition discrimination from longitudinal development of secondary particles, and their comparison with those from lateral distribution are studied in section 4. The differences between $EPOS-LHC$ and $QGSJet-II-04$ models are presented in section 5, and finally section 6 presents the summary.

The cascade processes of primary particles in the atmosphere were simulated with CORSIKA software (v77410) [24]. For high-energy hadronic interaction processes, the $EPOS-LHC$ and $QGSJet-II-04$ models are employed, the main results are based on $EPOS-LHC$ model. The $Fluka$ and $EGS4$ models were employed for low-energy hadronic interaction processes and electromagnetic interaction processes, respectively. There are five mass groups, namely proton, helium, CNO group (the mass number, noted as A, is 14), MgAlSi group (A is 27) and iron. All of which were simulated with fixed energy at $log_{10}(E/GeV)$=5.1, 5.3, 5.5, 5.7, 5.9, 6.1, 6.5, 6.9, 7.3, 7.7, 8.1, 8.5 and 8.9, respectively, corresponding to the energy range from below the “knee” to approaching the “ankle”. The altitude of detection plane was set to 4400 m a.s.l., and the horizontal and vertical components of the earth magnetic field at the observation site were set to be 34.618 $\mu T$ and 36.13 $\mu T$, respectively. Because Cherenkov light has different properties of longitudinal development for different zenith angles (denoted as $\theta$ below), to reduce the amount of simulation, three zenith angles ($\theta$=$0^{\circ}$, $28^{\circ}$ and $45^{\circ}$) were simulated, which corresponds to a maximum atmosphere depth of $\sim 850$ $g/cm^{2}$ with the choosing altitude. The azimuth angle of primary particles is uniformly distributed in the range of $[0^{\circ},360^{\circ}]$. The secondary components produced in EAS include Cherenkov light, electron, gamma rays, muon, neutron and so on [5]. The threshold energies for hadrons, muons, electrons and photons were set to 0.1, 0.1, 0.001 and 0.001 GeV, respectively. The longitudinal development of secondary particles was recorded in steps of 10 $g/cm^{2}$. Cherenkov light was produced and propagated to the observation level. To reduce the size of Cherenkov data in Cherenkov light simulation, only the energy of primary particles less than $10^{7.4}$ GeV were simulated, and only Cherenkov photons inside several regions at observation level were recorded. Each region corresponds to one telescope, and all photon bunches reaching spheres with 3m radius around the center of the region were detected and their position, and arrival direction were stored. In the production of secondary particles and Cherenkov light, the atmospheric refractive index and its evolution with height are taken into account, and the density of the atmosphere along the altitude is described by the US standard atmosphere model. The wavelength distribution of Cherenkov light is also simulated. The perpendicular distance between the telescope and shower axis (denoted as $R_{p}$ below) was fixed at $R_{p}$=50, 100, 150, 200, 300 and 400 m, respectively. In reality, Cherenkov light was absorbed inside the ozone region and scattered by atmosphere and aerosol. Telescopes were built to converge the Cherenkov light on the focal plane to have an image of the EAS, however, they were dependent on the atmosphere models, the aerosol models and the parameters of the telescope. This study focused on the physical properties of longitudinal development in CORSIKA data therefore, no telescope simulation was performed.

During the cascade process of the EAS, the number of secondary particles first increased as the EAS developed. It was then decreased due to energy loss. Moreover, two important parameters can be extracted from the longitudinal development. For electron and muon, one is the slant atmosphere depth at which the number of secondary particles per atmospheric depth bin reached maximum (denoted as $X^{max}_{e}$ for electron and $X^{max}_{\mu}$ for muon), another one is the number of secondary particles per atmospheric depth bin at $X^{max}_{e}$/$X^{max}_{\mu}$ (denoted as $N^{max}_{e}$ for electron and $N^{max}_{\mu}$ for muon). Meanwhile, for Cherenkov light, one is the slant atmosphere depth at which the number of Cherenkov photons per atmospheric depth bin observed by a telescope located $R_{p}$ away from the shower axis reached maximum (denoted as $X^{max}_{Cer}$), another one is the number of Cherenkov light at $X^{max}_{Cer}$ (denoted as $N^{max}_{Cer}$. As can be seen below, $X^{max}_{Cer}$ is not an intrinsic properties of air shower, it is different for telescopes located at different distance to the shower axis, so a $R_{p}$ value is attached when $X^{max}_{Cer}$ is giving. In general, $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$ can be used to distinguish between different composition. $N^{max}_{e}$, $N^{max}_{\mu}$ and $N^{max}_{Cer}$ are good energy estimators. In this section, the longitudinal developments of electrons, muons and Cherenkov light will be discussed.

For ground cosmic ray experiment, there are two kinds of measurements. One is measuring the secondary components through their lateral distribution at a fixed altitude, the properties of which have been studied in detail by  [5]. Another is measuring the longitudinal developments of secondary components along the trajectory (the properties of longitudinal development have been studied in section 3.) Energy measurement and composition discrimination ability are two of the most important aspects for cosmic ray measurements. The comparison of those two measurements between variables from longitudinal development and variables from lateral distribution will be discussed below.

Electron and Cherenkov light from both longitudinal development and lateral distribution are widely used for energy reconstruction. For electron, it was pointed out that the density of electron at 200 m away from the shower axis is less affected by the shower-to-shower fluctuations (see Fig. 11 in reference [5] for details), and it is a better energy estimator compared to total electron number. The comparison between the sigma value of the Gaussian function fitted to $ln(\rho_{e})$ (electron density) distribution for different zenith angles ($\theta=0^{\circ}$, $28^{\circ}$ and $45^{\circ}$) and the sigma value of the Gaussian function fitted to $ln(N^{max}_{e})$ distribution is shown in Fig. 12. Since the energy and zenith angle of primary particles are both fixed, the sigma value represents the shower-to-shower fluctuations, which is related to the reconstructed energy resolution.

The left and right panels of Fig. 12 correspond to primary proton and iron respectively. As illustrated in this figure, $N^{max}_{e}$ has a smaller fluctuation, and it is quite flat with energy. The shower-to-shower fluctuations of $\rho_{e}$ has a clear dependence on energy, which is due to the fluctuations of secondary particles having reached minimum at shower maximum. The measurement before and after shower maximum will have larger fluctuations. However, by choosing appropriate zenith angle for the desired energy, the shower-to-shower fluctuations of $\rho_{e}$ is very close to that of $N^{max}_{e}$. For example, for proton-induced shower below 10 PeV, the shower maximum is smaller than the slant atmospheric depth for all zenith angles, and the $\theta=0^{\circ}$ is closer to the shower maximum, having a smaller fluctuation of $\rho_{e}$ than $\theta=45^{\circ}$. For the proton-induced shower above 50 PeV, the shower maximum is larger than the slant atmospheric depth for $\theta=28^{\circ}$, and is smaller than that for $\theta=45^{\circ}$, so $\theta=45^{\circ}$ has a smaller fluctuation of $\rho_{e}$ than other zenith angles. For iron below 100 PeV, the situation is different, the shower maximum is always smaller than atmospheric depth for all zenith angles, so that $\theta=0^{\circ}$ has a smaller fluctuation of $\rho_{e}$ than other zenith angles.

Cherenkov light is another widely used method for energy reconstruction. There are two variables for Cherenkov light. The first is the total number of Cherenkov photons at observation level (denoted as $N_{Cer}$ below). Another one is the number of Cherenkov photons per atmospheric depth bin at its shower maximum (or $N^{max}_{Cer}$). The comparison of the shower-to-shower fluctuations between $N^{max}_{e}$, $N^{max}_{Cer}$ and $N_{Cer}$ with different $R_{p}$ values ($R_{p}$=150, 200, 300, 400 m) and different composition for $\theta=0^{\circ}$ and $\theta=45^{\circ}$ are shown in Fig. 13. For proton with $\theta=0^{\circ}$, at higher energy, the longitudinal development is more fluctuated, since there are more fraction of events with the shower maximum larger than the vertical atmospheric depth. As clearly seen, $N^{max}_{e}$ from electron longitudinal development has the best resolution in the whole energy range. The Cherenkov light measurement at fixed observation level (or $N_{Cer}$) has better resolution compared to $N^{max}_{Cer}$. This implies that accumulative measurement of Cherenkov light on observation level is a better choice for energy measurement.

For ground cosmic ray experiment, the composition discrimination power is limited compared to space cosmic ray experiment. Therefore, it is important to find powerful variables for composition discrimination. The variables from lateral distribution for composition discrimination have been studied in detail  [5]. It was extensively discussed that the density of muon at 250 m away from the shower axis (denoted as $\rho_{\mu}$ below) is the most powerful variable for nuclei discrimination. The shower-to-shower fluctuations of $\rho_{\mu}$ is important for the composition discrimination ability. The comparison of shower-to-shower fluctuations between $\rho_{\mu}$ and $N^{max}_{\mu}$ at different zenith angles for proton and iron is shown in Fig. 14.

Compared to electron in Fig. 12, the difference in the shower-to-shower fluctuations between the $\rho_{\mu}$ and $N^{max}_{\mu}$ is much smaller in Fig.14. The shower-to-shower fluctuations of $\rho_{\mu}$ is closer to that of $N^{max}_{\mu}$. This is due to the fact that muon almost do not interact with atmosphere, and the fluctuations of muon number after shower maximum is similar to the one at shower maximum. Since $N^{max}_{\mu}$ is more difficult to measure compared to $\rho_{\mu}$, the discrimination ability of $\rho_{\mu}$ will be compared to shower maximum from the longitudinal development point of view.

The variables $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$ from longitudinal development were widely used as a composition estimators. The distribution of $X^{max}_{e}$, $X^{max}_{\mu}$ from longitudinal development and $\rho_{\mu}$ from lateral distribution for proton and iron at different energies are shown in Fig. 15. As it can be seen, $X^{max}_{e}$ has similar ability with $X^{max}_{\mu}$, but $\rho_{\mu}$ has a much better composition discrimination power than $X^{max}_{e}$ and $X^{max}_{\mu}$. So $\rho_{\mu}$ from lateral distribution is a very powerful composition discrimination variable for cosmic ray physics.

The shower maximum measured from Cherenkov light (or $X^{max}_{Cer}$) is an estimator of $X^{max}_{e}$, which can also be used for composition discrimination. The distribution of $X^{max}_{Cer}$ for proton and iron at different $R_{p}$ values and different energies are shown in Fig. 16. As seen in the figure, the composition separation power of $X^{max}_{Cer}$ is worse than that of $X^{max}_{e}$ at around 100TeV, while they have similar discrimination power at around 1PeV.

Energy spectrum of cosmic ray measured by ground-based experiments depends on hadronic interaction models. For example, the proton energy spectrum measured by KASCADE experiment between $QGSJet-01$ and $Sibyll-2.1$ models are different, with the flux difference as nearly twice [19]. According to [5], the differences of lateral distribution of secondary particles between $EPOS-LHC$ and $QGSJet-II-04$ hadronic models have been studied. The effects of hadronic interaction models on longitudinal development will be studied in this section and compared with the effects of hadronic interaction models on lateral distribution.

The difference in percentage (defined as $(QGSJet-EPOS)/EPOS\ \times 100\%$) of $N^{max}_{e}$, $N^{max}_{\mu}$, $N^{max}_{Cer}$ , $N_{Cer}$, $\rho_{e}$ and $\rho_{\mu}$ between the two models are shown in Fig.17. As depicted in figure 17, the differences of all variables except $N^{max}_{\mu}$ is within $\pm 6\%$, while for $N^{max}_{\mu}$, the difference is larger for higher energy, around 9% for $log_{10}(E/GeV)$=8.1. The difference of $\rho_{\mu}$ which is smaller than $N^{max}_{\mu}$ at high energy could be due to the difference of muon number which is mainly in the forward region of the shower (see Fig. 8 of reference [5] for details), and the $\rho_{\mu}$ measured far away from the shower core. It is also noticed that the differences between $EPOS-LHC$ and $QGSJet-II-04$ hadronic models are very similar for different type of nuclei for $log_{10}(E/GeV)$=8.1.

The location of the shower maximum of an air shower is another measurement that characterizes the properties of longitudinal development. The differences of $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$ between hadronic models are shown in Fig.18. As illustrated, when $log_{10}(E/GeV)<=$ 6.9, the differences of $X^{max}_{e}$ and $X^{max}_{\mu}$ is within 6 $g/cm^{2}$, except $X^{max}_{\mu}$, which is within around 10 $g/cm^{2}$.

The shower-to-shower fluctuations of shower size will affect the composition separation power and energy resolution. The upper panel of Fig.19 shows the sigma value of a Gaussian function fitted to the distributions of $N^{max}_{e}$, $N^{max}_{\mu}$, $N^{max}_{Cer}$, $N_{Cer}$, $\rho_{e}$ and $\rho_{\mu}$ in percentage, while the lower panel shows the sigma value of a Gaussian function fitted to $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$. The samples are proton and iron nuclei with energy $log_{10}(E/GeV)=$ 6.9 simulated with $EPOS-LHC$ and $QGSJet-II-04$ hadronic models, respectively. As clearly observed, the differences of the sigma value for shower size is within $2.5\%$, while the differences of the sigma value for shower maximum is within $2g/cm^{2}$. Similar conclusions were found for other energies.

In general, for the difference of nuclei between $EPOS-LHC$ and $QGSJet-II-04$ models, the differences of the number of electron and Cherenkov light in longitudinal development are close to the differences of the number from lateral distribution at observation level. The muon density at $R_{p}$=250 m (or $\rho_{\mu}$) has a smaller hadronic interaction model dependence than $N^{max}_{\mu}$. When $log_{10}(E/GeV)<$ 6.9, the differences of $X^{max}_{e}$ and $X^{max}_{Cer}$ are within 6 $g/cm^{2}$, while the difference of $X^{max}_{\mu}$ is lager, with a maximum of around 10 $g/cm^{2}$. The shower-to-shower fluctuations both for the shower size and shower maximum are small between $EPOS-LHC$ and $QGSJet-II-04$ models.

In this work, the properties of the longitudinal development of electron, muon and Cherenkov light produced in EAS in the knee energy region were studied, with a CORSIKA simulated samples for different energy, composition and zenith angles. The shower maximum ($X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$) and the number of secondary particles at shower maximum ($N^{max}_{e}$, $N^{max}_{\mu}$ and $N^{max}_{Cer}$) were derived from the longitudinal development of electron, muon and Cherenkov light by fitting them with a “Gaussian-in-age” function. The relationship between $X^{max}_{Cer}$ and $X^{max}_{e}$ is almost composition independent, however, the relationship is dependent on $R_{p}$ and zenith angle, therefore, $X^{max}_{Cer}$ is kind of a complicated composition independent $X^{max}_{e}$ estimator. The uncertainty of $X^{max}_{e}$ reconstruction from $X^{max}_{Cer}$ is about 10-15 $g/cm^{2}$ for all compositions and $R_{p}$ values, when energy is larger than$\sim$1 PeV. From 100 TeV to 1 PeV, the uncertainty is from 10-30 $g/cm^{2}$ for $R_{p}$=150 m, and 10-45 $g/cm^{2}$ for $R_{p}$=400 m.

The performances of energy measurement and composition discrimination ability from $N^{max}_{e}$, $N^{max}_{\mu}$, $N^{max}_{Cer}$, $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$ from longitudinal development were studied and compared with the variable $\rho_{e}$, $\rho_{\mu}$ from lateral distribution. It was found that $N^{max}_{e}$ has smallest shower-to-shower fluctuations compared to number of Cherenkov light and the density of electron (or $\rho_{e}$) at observation level. The smallest shower-to-shower fluctuations of $N^{max}_{e}$ is around 5% for nuclei. However, by choosing the appropriate zenith angle at different energy, the fluctuations of $\rho_{e}$ is very close to the fluctuations of $N^{max}_{e}$. Meanwhile, the combination of electron number and muon number at their shower maximum, namely $N^{max}_{e}$+25$N^{max}_{\mu}$, is a very good composition independent energy estimator. $X_{max}$ from longitudinal development can be used for composition discrimination, which is linear with lnA. The shower-to-shower fluctuations of $X^{max}_{e}$ is 50—55 $g/cm^{2}$ for proton, and 20—25 $g/cm^{2}$ for iron. For proton/iron separation, the composition separation power of $X^{max}_{Cer}$ is worse than $X^{max}_{e}$ at around 100 TeV, while they have similar discrimination power at around 1PeV. However, the density of muon (or $\rho_{\mu}$) at observation level is much better than $X^{max}_{e}$, $X^{max}_{\mu}$ and $X^{max}_{Cer}$.

The differences of the variables both from longitudinal development and lateral distribution between $EPOS-LHC$ and $QGSJet-II-04$ hadronic model were studied. For the differences between $EPOS-LHC$ and $QGSJet-II-04$ models, the differences of the number of electron and Cherenkov light in longitudinal development are close to the differences of the variables from lateral distribution at observation level, while the $N^{max}_{\mu}$ from longitudinal development has a larger model dependencies than $\rho_{\mu}$ from lateral distribution. The differences of $X^{max}_{e}$ were also similar with $X^{max}_{Cer}$, but the difference of $X^{max}_{\mu}$ is larger than both of them. The differences of shower-to-shower fluctuations between $EPOS-LHC$ and $QGSJet-II-04$ models are negligible.

Therefore, our results provide important information for the selection of detector type, detector design and data analysis for physical measurements.

This work is supported by the Science and Technology Department of Sichuan Province (Grant number 2021YFSY0031), the National Natural Science Foundation of China (Grant number 12205244), the National Key R&D Program of China (Grant number 2018YFA0404201). The numerical calculations in this paper have been done on Hefei advanced computing center.