Characteristics of extensive air showers around the energy threshold for ground-particle-based $\gamma$-ray observatories

Abstract

Very-high-energy $\gamma$-ray astronomy based on the measurement of air shower particles at ground-level has only recently been established as a viable approach, complementing the well established air Cherenkov technique. This approach requires high (mountain) altitudes and very high surface coverage particle detectors. While in general the properties of air showers are well established for many decades, the extreme situation of ground-level detection of very small showers from low energy primaries has not yet been well characterised for the purposes of $\gamma$-ray astronomy. Here we attempt such a characterisation, with the aim of supporting the optimisation of next-generation $\gamma$-ray observatories based on this technique. We address all of the key ground level observables and provide parameterisations for use in detector optimisation for shower energies around 1 TeV. We emphasise two primary aspects: the need for large area detectors to effectively measure low-energy showers, and the importance of muon identification for the purpose of background rejection.

1 Introduction

In the past two decades very-high-energy $\gamma$-ray astronomy has become established as a rich astronomical discipline, based primarily on observations with Imaging Atmospheric Cherenkov Telescope (IACT) arrays. However, given the limitations of IACTs in terms of Field-of-View (FoV, typically a few degrees diameter) and duty cycle (limited to typically 10–15%), new approaches are being developed to provide very large field of view and continuous observations. The direct detection of extensive air shower (EAS) particles at ground level is the most obvious complementary approach. Air Shower Particle Detectors (ASPDs), can observe a large fraction of the overhead sky at any given time and, with their close to 100% up-time, can survey roughly 2/3 of the sky on a daily basis. However, since they measure only the ground-level properties of air showers, the accuracy with which they reconstruct the properties of the primary $\gamma$-ray is typically significantly worse than is the case for IACTs. Typical ASPD (IACT) performance values around 1 TeV are angular resolution $0.4^{\circ }$ (CF $0.06^{\circ }$), hadron rejection efficiency 99% (CF 99.9%) and energy resolution of 50% (CF 10%) – see [1,2,3,4,5].

During their propagation through the atmosphere, the number of particles in the air shower grows due to production of secondary particles until an atmospheric depth $X_\text {max}$. After this maximum, the energy content of the shower declines as ionisation losses dominate over new particle production. ASPDs are typically located at high altitude sites to be as close to $X_\text {max}$ as possible and maximise the particle count per shower, improving the detection efficiency and accuracy. Ground-based $\gamma$-ray astronomy relies on discrimination of $\gamma$-rays from the dominant background of cosmic ray protons and nuclei on the basis of air shower characteristics. The most obvious differentiating characteristic of hadronic versus purely electromagnetic (EM) cascades is the presence of pions, and subsequent production of muons in pion decay as well as EM sub-showers. For a recent detailed review of the properties of air showers we refer the reader to [6].

Several instrumentation approaches have been used for ASPDs at different observatories around the world: the Tibet AS-$\gamma$ experiment [7] used scintillator-based detectors, ARGO-YBJ is based on Resistive Plate Chambers (RPCs) [8], but up to now, the most successful technique in terms of sensitivity, angular and energy resolution is the water Cherenkov technique introduced in MILAGRO [9] and currently implemented in the High Altitude Water Cherenkov (HAWC) $\gamma$-ray observatory [10]. HAWC has reached the critical sensitivity required to detect $\gamma$-ray sources in significant numbers [11]. The Large High Altitude Air Shower Observatory [12] is a facility which is currently under construction and is expected to become operational in the coming years. In addition, there is a significant push towards an ASPD in the Southern Hemisphere, with a range of concepts being explored [13,14,15,16,17,18]. Here we study the properties of $\gamma$-ray and background showers around the threshold where the ASPD approach becomes possible for a high altitude detector, focussing on aspects that influence the design and utilisation of such a detector.

2 Simulations and atmospheric propagation

$\gamma$-ray and proton induced air showers are simulated using the CORSIKA (COsmic Ray SImulations for KAscade) 7.4005 package [19]. Within the CORSIKA package, the hadronic interaction event generator FLUKA [20, 21] is used below 80 GeV and above this energy the QGSJet-II model [22] is used. The choice of hadronic interaction models can influence the gross development of the simulated shower. However a detailed study of these effects is beyond the scope of the work presented here. For a recent systematic study on the impact of the choice of hadronic interaction models on air shower simulation in the energy range of our interest we refer to [23]. The electromagnetic interactions are handled by the EGS4 model [24]. The low kinetic energy cut-off values are chosen to be 50 MeV for hadrons and muons and 0.3 MeV for electrons and photons.

With this configuration a library of air shower simulations is generated for proton and $\gamma$ ray induced air showers with discrete values in primary energy, zenith angle, and altitude of ground-level. For each set of parameters $10^{4}$ air showers are simulated, and the values of the used parameters will be explicitly mentioned at each of the figures showing analysis results. The discrete values of zenith angles and ground-levels used in the simulations are indicated on Fig. 14 in Appendix A and are chosen such that they cover the plausible altitude and zenith range for a future ASPD observatory. For the charged background of cosmic nuclei we considered only protons, as they are the dominant source of $\gamma$-ray like background events in ASPD $\gamma$-ray observatories in this energy range. The parameter that drives the number of interactions and energy losses in the shower development is the amount of material that the particles in the EAS traverse during their propagation through the atmosphere, the slant depth, given by the line integral over the atmospheric density $\rho _{\text {atm}}$ from a location of altitude z above sea level in the direction of shower origin: $X(z,\theta ) = \int _\infty ^{z}\rho _{\mathrm {atm}}(z',\theta )\cos \theta dz'$, where $\theta$ is the zenith angle. Figure 14 in Appendix A provides a reference for the relationship of slant depth to altitude (z) and zenith angle $\theta$ for the US standard atmosphere density profile, for use in interpreting the results given here in terms of slant depth. At all feasible altitudes for an observatory in this energy range, shower maximum typically occurs well before the shower reaches the ground. The dependency of the average value of $X_{\text {max}}$ on energy, usually called the elongation rate, can be parametrised as $\langle X_{\text {max}} \rangle = a + b \log _{10}(E/\text {GeV})$ [25, 26], where E is the energy of the primary particle. For simulated $\gamma$-ray showers we find $a = 98\hbox { g~cm}{^{-2}}$ and $b=83\hbox { g~cm}{^{-2}}$ and for proton primaries $a = 111\hbox { g~cm}{^{-2}}$ and $b =74\hbox { g~cm}{^{-2}}$. This means that for an observatory at an altitude of 5 km a 100 GeV (1 TeV) $\gamma$-ray induced shower from zenith is on average 7.7 (5.5) radiation lengths beyond the atmospheric depth at which it reached $X_{\text {max}}$, implying very few particles reaching ground-level and large fluctuations.

Fig. 1

Distribution of the fraction $x_{\text {gr}} = E_i/E_{\text {gr}}$ of the total energy that reaches ground $E_{\text {gr}}$ per particle type $i=(e^{\pm }, \mu ^{\pm }, \gamma$ or other) for vertical EASs with a ground altitude of 5 km above sea level. The left panels show the distributions for $\gamma$-ray primaries while proton primaries are shown in the right panels. Top and bottom panels correspond to 100 GeV and 1 TeV primary energy respectively

3 Ground-level particles

The vast majority of the energy in EASs is carried to the ground by photons, electrons and positrons (hereafter simply electrons), and muons ($\mu ^{+}$ and $\mu ^{-}$). Figure 1 shows for each of these particle types the distribution of the fraction of the total energy arriving at the ground ($x_{\text {gr}} = E_{i}/E_{\text {gr}}$) for 100 GeV and 1 TeV proton and $\gamma$-ray primaries. For $\gamma$-ray showers the well-known dominance of photons is apparent, as well as the occasional production of muons (and even very occasionally the energetic dominance at the ground of muons). For proton initiated showers the majority of the arriving energy is typically in the form of muons. For both primary particle types huge fluctuations are evident at these energies, as expected due to the small number of particles reaching the ground. We note that the category “other” is dominated by protons and neutrons, which can carry a significant fraction of the energy in rare cases.

The distributions of individual particle energies at ground-level are very different for muons with respect to electrons and photons, as illustrated in Fig. 2. The peak of the energy distribution in terms of number per log energy interval($dN/d\log E$) is around $\sim$ 6 MeV for photons, $\sim$ 20 MeV for electrons and 2–3 GeV for muons. In terms of total energy per log interval ($EdN/d\log E$) the peak lies at $\sim$ 150 MeV for photons, $\sim$ 600 MeV for electrons and 30–40 GeV for muons. The shape of these distributions are similar for showers initiated by primary particles with different energies.

Fig. 2

Number density (top) and energy density distributions (bottom) of ground-level particles at 5 km altitude for vertical showers initiated by primary protons of three energies, solid lines correspond to 1 TeV primary energy, dashed lines to 10 TeV, and dashed-dotted lines to 100 GeV. The distributions of the 10 TeV case have been scaled down by a factor of ten, while the 100 GeV curves have been multiplied by a factor of ten

Fig. 3

Distributions of the number of muons within 100 m of the impact point as a function of ground-level electromagnetic energy for vertical showers observed at 5 km altitude. Proton (blue) and $\gamma$-ray (red) initiated showers of different energies are shown, with the area of the squares scaling linearly with density

Many ASPDs make use of calorimetric detectors that are sensitive to the total ground-level electromagnetic energy ($E_{\text {em}}$), rather than the number of electrons. Muons, in contrast, rarely loose a significant amount of their energy in typical detectors and hence their number ($N_{\mu }$) is the relevant quantity. Figure 3 shows the distributions of muon number $N_{\mu }$ versus electromagnetic energy $E_{\text {em}}$ for proton and $\gamma$-ray initiated showers of several energies. A large difference in $N_\mu$ at fixed $E_{\text {em}}$ between $\gamma$-ray and proton induced EASs is apparent, with useful separation power appearing at energies above about 1 TeV provided both of these quantities can be adequately measured.

In this paper only proton simulations are considered for charged cosmic rays. The cosmic ray flux is dominated in the energy range of interest by protons and helium, with helium having on average a slightly higher ratio of $N_\mu /E_{\text {em}}$ which makes them easier to discriminate from $\gamma$-ray showers than protons. The contribution of heavier elements (like Iron for example) are significantly reduced below a TeV primary energy as very few particles reach ground level. At higher energy they would appear as air showers with a relatively high ratio of $N_\mu /E_{\text {em}}$ and are therefore easier to identify as background than proton or helium initiated showers.

4 Parameterisation of electromagnetic energy at ground

Fig. 4

Distributions of the electromagnetic energy that reaches the ground at 5 km altitude for vertical showers, expressed as a fraction of the primary particle energy $x_{\text {em}} = E_{\text {em}} / E$. $\gamma$-ray initiated showers are shown in the top panel and protons in the bottom panel

Figure 4 shows the distribution of the fraction of the primary energy which reaches the ground in the electromagnetic part of EASs ($x_{\text {em}}$) for different energies and primaries. It is clear that for proton-initiated showers not only does less electromagnetic energy reach the ground, but that fluctuations from shower to shower are a lot larger than for $\gamma$-ray showers. This is especially the case for the lowest energy showers considered, where just a few muons may carry the majority of the total energy at ground. In addition, we observe that the quantity $\log _{10} (x_{\mathrm{{em}}})$ follows an approximately Gaussian distribution. Therefore, the distributions in $\log _{10} (x_{\mathrm{{em}}})$ can be parametrised to first order by the root mean square (rms) $\sigma (\log _{10} (x_{\mathrm{{em}}}))$ and the mean $\langle \log _{10} (x_{\mathrm{{em}}})\rangle$ of the distributions.

Fig. 5

Dependence of fractional electromagnetic energy arriving at ground $x_{\mathrm{{em}}} = E_{\text {em}} / E_{\mathrm{\gamma }}$ on slant depth, for primary $\gamma$-rays of different energies. The top panel shows the dependence of the mean of the $\log _{10} (x_{\mathrm{{em}}})$ distribution, while the bottom panel shows the behaviour of the rms ($\sigma$) of the same distribution. Each data point corresponds to the slant depth of a combination of zenith angle and ground-level from Fig. 14, filled markers correspond to simulation sets that have a zenith angle larger than $30^{\circ }$. Statistical uncertainties are of similar size as the markers or smaller

Fig. 6

Dependence of fractional electromagnetic energy arriving at ground $x_{\mathrm{{em}}} = E_{\text {em}} / E_{\mathrm{p}}$ on slant depth, for primary protons of different energies. The top panel shows the dependence of the mean of the $\log _{10} (x_{\mathrm{{em}}})$ distribution, while the bottom panel shows the behaviour of the rms ($\sigma$) of the same distribution. Each data point corresponds to the slant depth of a combination of zenith angle and ground-level from Fig. 14, filled markers correspond to simulation sets that have a zenith angle larger than $30^{\circ }$. Error bars are of similar size as the markers or smaller

The resulting characterisation of the distributions of shower $x_{\text {em}}$ as a function of the slant depth is shown in Fig. 5 for $\gamma$-ray, and in Fig. 6 for proton, initiated showers. As expected, the mean $E_{\mathrm{{em}}}$ decreases rapidly with increasing slant depth. For the lowest energy showers this results in a fraction of simulated showers with no particles arriving at the ground which cannot be used in the estimation of mean and rms in log-space. If this fraction is more than 10% of the total sample, we omit the sample completely from Figs. 5 and 6. The dependencies of the mean and the width of the distributions on slant depth are reasonably well described by first or second order polynomials as illustrated by the fit functions represented by the lines in Figs. 5 and 6. The fit parameters are given in Appendix B.

Fig. 7

The fraction of showers that have $E_{\text {em}}>10$ GeV as a function of slant depth. Top: $\gamma$-ray induced showers. Bottom: proton induced showers. The markers are obtained from the full simulation set, while the lines are obtained using the parameterisation based on the fitted behaviour from Figs. 5 and 6. Each data point corresponds to the slant depth of a combination of zenith angle and ground-level from Fig. 14

Well above detection threshold effects, the collection area of an ASPD for well-measured showers can be considered equal to the projected footprint of the array. Around the threshold however, the collection area is strongly zenith angle and primary energy dependent, with showers fluctuating deep in the atmosphere providing some collection area at the lowest accessible energies. To first approximation the detectability of a $\gamma$-ray shower at ground level depends on the arriving EM energy $E_{\text {em}}$, and the collection area on the fraction of showers that have $E_{\text {em}}$ above a threshold value for a given primary energy and slant depth. In Fig. 7, this fraction is shown as a function of slant depth, with an assumed threshold $E_{\text {em}}>$10 GeV. This threshold has been chosen as an example of an optimistic detection threshold for a future observatory. Figure 7 compares the results of the full Monte–Carlo shower simulations (markers) with the parameterisations given in Appendix B (lines, see Figs. 5 and 6). The agreement between the full Monte–Carlo and simple parameterisation is within a few percent for gamma-ray and within 10% for protons, illustrating the usefulness of the parameterisations for first order detection threshold approximations as presented here. The parameterisation for proton showers overestimates the fraction of showers with $E_{\text {em}} > 10$ GeV by $\sim$ 5%, which is likely due to a deviation from a pure gaussian behaviour. For first order approximations the size of these deviations are not very relevant. However when estimating detection efficiency far in the tails of these distributions it is advised not to rely on simple parameterisations, but rather deploy a full Monte Carlo and accurate detector simulation. As for Figs. 5 and 6, the dependency on slant depth is typically stronger for the $\gamma$-ray induced air showers, which is expected for pure electromagnetic cascades.

5 Lateral extent

In addition to the total $E_{\text {em}}$ that reaches the ground, the area over which this energy is spread is a crucial parameter in shower detectability for a given array design. We adopt the radius $r_{50}$, in which 50% of the total electromagnetic energy that reaches the ground is contained, as an indicator of the shower extent. $r_{50}$ is calculated in the plane perpendicular to the direction of propagation of the primary particle. The choice of 50% rather than a larger containment fraction is motivated by the large values obtained (even for half containment) for low energy showers. In addition, for air shower reconstruction a substantial fraction of the footprint needs to be measured, but not necessary all the way to sparsely populated tail of the lateral particle distribution. Therefore, $r_{50}$ relates more closely to the desirable minimum size of a particle array than for example $r_{90}$. We note that the Molière radius contains approximately 90% of the energy in the full electromagnetic cascade [27]. In air at 5 km above sea level the Molière radius $R_{m}$ is approximately 130 m.

Fig. 8

Relationship of the lateral extent parameter $r_{50}$ to other shower parameters for vertical EAS observed at 5 km altitude. Results are shown for $\gamma$-ray (left panels) and proton (right panels) primary particles simulated at two different energies. From top to bottom the following EAS parameters are shown: a The average number of particles at ground with energy of more than 1 MeV. b Average EM energy at ground. In addition, the bands on this panel indicate the range that contains 90% of $r_{50}$ values. c The maximum energy of a single EM-particle ($e^\pm ,\gamma$). d Depth of shower maximum, obtained from the fit routine within CORSIKA. The parameters shown in a–c are scaled with the energy of the primary particle

Figure 8 illustrates the relationship of the parameter $r_{50}$ to several other related EAS characteristics. For a fixed primary energy, $r_{50}$ depends on the stage of development of the EAS, which can be clearly observed from the correlations with the EAS parameters shown in Fig. 8. From the behaviour of $r_{50}$ with respect to the number of particles at ground N (panel a) it is clear (at least for 5 TeV showers) that the smallest values of $r_{50}$ are provided by showers which reach ground-level before reaching their maximum development.¹ However, this effect is not so clear from the relationship to shower maximum $X_{\text {max}}$ (panel d) as obtained by fitting the longitudinal shower development with a Gaisser-Hillas function. The reason for this is presumably the difficulty of extrapolating the longitudinal profile to a below-ground $X_{\text {max}}$, in particular for small showers with significant fluctuations. For showers that reach $X_{\mathrm{max}}$ well before arriving to the ground $r_{50}$ does show the expected correlation with shower maximum.

A better handle on the stage of EAS development might be obtained by the energy $E_{\text {max}}$ (panel c) of the most energetic single electromagnetic particle in the shower at ground, which correlates to the number of interactions in the EM-cascade before the ground is reached. Whilst $E_{\text {max}}$ is not a straight-forward parameter to measure, Fig. 8 indicates that both it and $r_{50}$ may be very useful observables as indicators of the stage of shower development. This idea is reinforced by the fact that the behaviour of $r_{50}$ with the fraction of energy arriving at the ground (Fig. 8b) is very similar for the two primary energies illustrated.

Fig. 9

Lateral shower extent parameter $r_{50}$ as a function of slant depth. The points correspond to the median value of $r_{50}$. The top panel shows the dependency for $\gamma$-ray showers while the behaviour for protons is shown in the bottom panel. The results of fitting a second order polynomial are shown as solid lines, fit parameters are given in Appendix B. Each data point corresponds to the slant depth of a combination of zenith angle and ground-level from Fig. 14. Statistical uncertainties are of similar size as the markers or smaller

Figure 9 shows the dependence of the median $r_{50}$ on slant depth, indicating again the very large typical extent of showers well beyond shower maximum. Again the relationship is fit with a second order polynomial and the best fit parameters provided in Appendix B for both gammas and protons. The typical few percent scatter around these lines give an indication of the accuracy of this parameterisation.

For the lowest energy showers, fluctuations are very large and the trends are less clear, but clearly the typical extent is very large in comparison to the traditional regime of air shower measurements. For higher energy showers closer to $X_{\mathrm{max}}$ we find that the 90% energy containment radius $r_{90}$ is close to the Molière radius as expected.

The deviation from this behaviour as the shower peters out is consistent with the disappearance of particles at or above the critical energy from which the Molière radius is defined, with subsequent rapid multiple scattering and large displacements from the shower axis for the remaining low energy electrons. As an example, considering 100 GeV $\gamma$ rays with a typical $r_{50}$ of 150 m, it is apparent from Fig. 8 that such showers have typically 100 particles sharing about 2 GeV of kinetic energy and no particles with energy $>200$ MeV. Another consequence of the low particle number and large scattering is that the lateral distribution of particles becomes very flat and with very large fluctuations. In general we find that $r_{90}\sim 6 r_{50}$, but again with large fluctuations for the low energy cases. To give an indication of the typical appearance of showers of these energies at ground-level, Fig. 10 compares the lateral distribution of EM energy for examples of individual $\gamma$-ray initiated showers of different energies. For the two 1 TeV examples the distributions are clearly centrally peaked and rather similar, but at 100 GeV fluctuations are very large, shower-to-shower difference are very large, and no clear shower core may be present.

Fig. 10

Electromagnetic energy density $\rho$ as a function of distance to the shower axis for four individual vertical $\gamma$-ray induced showers of two different energies observed at an altitude of 5 km

6 Muon number

As is clear already from Fig. 3, the number of muons at ground relative to electromagnetic energy provides a powerful discriminant between $\gamma$-ray and proton showers. As well as the number of arriving muons, the lateral extent of muons is a key consideration. Muons are often present at rather large distances to the shower axis as they gain significant transverse moment from both the parent pions and the decay process. For showers in the energy range considered here, the shower axis is often poorly defined based on ground-level observables. Here we consider the distance of muons from the barycentre of the electromagnetic particles in the shower, as the key factor for discrimination is the association of a given muon with the EM-component of the shower. Given the rather high rate of background muons unassociated to any shower at ground, it is unlikely that the association of a single muon can be used for discrimination purposes. In the following we focus on the fraction of showers in which at least two muons reach ground level, within radii, 50 m, 100 m and 200 m, of the EM centroid. Figure 11 illustrates these fractions as a function of both primary energy and electromagnetic energy at ground for proton initiated showers. In addition, we provide here for reference the equivalent $\gamma$-ray energies which result in the same mean $E_{\text {em}}$ for this zenith angle and altitude.

Fig. 11

Top: fraction of proton induced showers that have more than two muons within 50/100/200 m distance of the shower EM barycentre as a function of primary particle energy. Bottom: the same fraction but as a function of the electromagnetic energy brought to the ground ($E_{\text {em}}$). The upper label for the abscissa indicates the average energy of a vertical incoming $\gamma$-ray primary which results in the same $E_{\text {em}}$ as observed at 5 km altitude. The 4.1 km above sea level corresponds to the altitude of the HAWC observatory

As showers are observed at progressively larger slant depths, the beam of muons broadens and an increasing number of muons decay before reaching the ground, both effects leading to a lower muon count per unit area. The upper panel of Fig. 11 illustrates this reduction in muon number per unit area with decreased observation altitude, but as is clear in the lower panel of this figure, the attenuation effect is much stronger for the electromagnetic component. It is apparent that a detector linear scale of $\gg$50 m is needed to exploit muons for background rejection for $\gamma$-rays below 1 TeV.

7 Zenith angle dependence

In the previous sections we have shown mainly results for vertical showers or as a function of the slant depth of the observation altitude. The slant depth determines the total amount of target material and therefore relates directly to the number of interactions within the particle cascade. In this section we will illustrate features of the ground-particle distributions at a fixed slant depth of $\hbox {800 g cm}^{-2}$, while changing the zenith angle and observation altitude as indicated in Fig. 14. Figure 12 shows the mean and the width of the distributions that describe the amount of energy that reaches ground in electromagnetic particles (equivalent to Figs. 5 and 6) as a function of zenith angle. There is no significant dependence of these distributions on the zenith angle, justifying to parameterise these distributions only as a function of slant depth.

Fig. 12

The dependence of the fractional electromagnetic energy arriving at ground $x_{\mathrm{em}} = E_{\text {em}} / E_{\mathrm{p}}$ on zenith angle ($\theta$), for primary protons and $\gamma$ rays of 1 TeV. The top panel shows the dependence of the mean of the $\log _{10} (x_{\mathrm{{em}}})$ distribution, while the bottom panel shows the behaviour of the rms ($\sigma$) of the same distribution. Each data point corresponds to a slant depth of $\hbox {800 g cm}^{-2}$ and a combination of zenith angle and ground-level from Fig. 14, ground altitudes are between 2.1 and 5.6 km and the zenith angle range is from $0^{\circ }$ and $55^{\circ }$. Error bars indicating the statistical uncertainty are typically of the size of the marker or smaller

Fig. 13

Top: dependence on the zenith angle $\theta$ of shower size parameter $r_{50}$ for a fixed slant depth of $\hbox {800 g cm}^{-2}$. Bottom: Same as top panel, but for the average number of muons $N_{\mu }$ in proton showers. Each data point corresponds to a slant depth of $\hbox {800 g cm}^{-2}$ and a combination of zenith angle and ground-level from Fig. 14, ground altitudes are between 2.1 and 5.6 km and the zenith angle range is from $0^{\circ }$ and $55^{\circ }$. Error bars indicating the statistical uncertainty are typically of the size of the marker or smaller

Figure 13 shows two parameters that have a clear dependence on zenith angle (or altitude) for fixed slant depth. At fixed slant depth the particle cascade develops in a lower density atmosphere with increasing zenith angle (or altitude) which effectively increased the physical length of the longitudinal shower development. As a consequence, particles have more time to diffuse away from the shower axis as is clearly illustrated by the linear increase of $r_{50}$ with $1- \cos (\theta )$ in the top panel of Fig.  13. This effect leads to discrepancies ranging from a few percent up to 10 percent from the parameterisations derived in Sect. 5 for zenith angles above 30${^\circ }$. Similar, a roughly linear decrease is observed for the number of muons per area. With increasing zenith angle muons are on average produced further away from the detection point which leads to larger lateral distance at ground and higher probability to decay before reaching ground level.

8 Implications for array design

For altitudes much above 5 km it becomes extremely difficult to find a suitable site for a $\gamma$-ray observatory. It therefore seems likely that a future observatory must be designed to deal with being $\sim$ 6 radiation lengths below $X_{\mathrm{max}}$ even for vertical showers. For primary $\gamma$-rays below 1 TeV this translates to low typical ground-level particle energies and hence as large, and strongly fluctuating, shower footprint (see Figs. 8 and 10). Given that the nominal 10 GeV ground-level EM energy threshold assumed for Fig. 7 is already extremely challenging, an array on a scale significantly larger than $r_{50}$ (i.e. $\sim$ 100 m), and with fill factor approaching 100%, is needed for $\gamma$-ray astronomy at these energies. For hadron rejection based on muon identification the relevant scale is even larger as can be seen from Fig. 11. A detector on the scale of HAWC (array area $\sim ~\hbox {20,000~m}^{2}$) will typically detect only a modest fraction of the arriving EM energy of a sub-TeV shower and a small fraction of the muons in a background TeV proton shower. Given the absence of a clear core for many far-beyond $X_{\mathrm{max}}$ low energy showers (see e.g. Fig. 10) morphology-based hadron rejection will be extremely challenging and clear muon identification appears to be a promising approach for most of the energy range considered here. Figure 11 indicates that muon counting can be exploited for $\gamma$-ray energies as low as 300 GeV with a sufficiently large array. The Large High Altitude Air Shower Observatory [12] currently under construction in China will have a large water-Cherenkov detector array of $\hbox {78,000~m}^{2}$ with a nearly full ground coverage. The sensitivity, especially at the low-energy threshold, is expected to scale faster than the square root of the instrumented area, as a larger fraction of the events will be well contained within the array (see Fig. 8) and more muons will be detected to be useful in background rejection (see Fig. 11). When taking Fig. 11 as a proxy for hadronic background rejection rates, LHAASO is expected to reach similar efficiencies to HAWC at a quarter decade lower gamma-ray energy. There is currently another serious effort ongoing to plan a next-generation facility in the Southern Hemisphere [18], hence there seems to be a bright future for the application of ASPDs in gamma-ray astronomy.

9 Conclusions

The parameterisations of proton and $\gamma$-ray showers at ground level presented here should prove useful in the design of next generation $\gamma$-ray observatories based on ground-particle detection. In the shower tail, far beyond $X_{\mathrm{max}}$ where sub-TeV measurements are possible, the typical extent of showers becomes very large and large array footprint as well as close to 100% fill factor is required. The tagging of individual muons offers an opportunity for background rejection even in the sub-TeV regime, where the very large fluctuations in $\gamma$-ray showers at ground will make rejection based on the lateral distribution of particles very difficult.

Appendices

Appendix A

Figure 14 provides a reference for slant depth $X(z,\theta ) = \int _\infty ^{z}\rho _{\mathrm {atm}}(z',\theta )\cos \theta dz'$ versus zenith angle and observation altitude for a US standard atmosphere. The markers indicate the discrete combinations of zenith angles and ground levels for which simulation are generated.

Fig. 14

Slant depth at ground level as a function of altitude and zenith angle for a US standard atmosphere. Each contour-line indicates a $\hbox {50~g~cm}^{-2}$ change. The cross markers indicate the discrete values at which simulations sets are generated in the relevant altitude and zenith angle range and the full circles indicate the sets used in Sect. 7

Appendix B

Table 1 Fit results to parameterise the distributions of $E_{\text {em}}$ as a function slant depth. Listed are the parameters from a polynomial fit f$(X_{\text {sl}}) = p_0 + p_1 X_{\text {sl}} + p_2 X_{\text {sl}}^2$, where $X_{\text {sl}}$ is the slant depth, and f is given in the first column of the header of each sub-table

In Table 1 we include the fit parameters from Figs. 5 and 6. These parameters can be used to derive the mean and standard distribution of the distributions at a given energy for gamma and proton primaries and the selected slant depth.

In Tables 2 and 3 we include the fit parameters from Fig. 9. These parameters can be used to estimate shower sizes as a function of slant depth. Especially for low energy results, it is not wise to extrapolate these functions outside the range indicated in Fig. 9.

Table 2 Best fit parameters for the function $r_{50}(X_{\text {sl}}) = p_0 + p_1 X_{\text {sl}} + p_2 X_{\text {sl}}^2$ for $\gamma$-ray initiated air showers. $X_{\text {sl}}$ is the slant depth

Table 3 Best fit parameters for the function $r_{50}(X_{\text {sl}}) = p_0 + p_1 X_{\text {sl}} + p_2 X_{\text {sl}}^2$ for proton initiated air showers. $X_{\text {sl}}$ is the slant depth