1 Introduction

The hadronic cascades associated to charged cosmic ray primary particles typically produce large numbers of muons, primarily from the decay of charged pions. The potential to use these muons to discriminate between hadronic and electromagnetic cascades and hence do gamma-ray astronomy, has long been recognised (see for example [1]).

Fig. 1
figure 1

Example simulated muon images in a 28 m telescope at different energies and impact distances. All simulated muons are produced at 11 km above sea level and are observed at 1835 m with a zenith angle of 20\(^\circ \)

Full size image

The performance of the LHAASO array [2, 3], demonstrates the power of this approach for very-high and ultra-high energy gamma-ray astronomy (from tens to hundreds of TeV). An important factor contributing to the success of the LHAASO array is the very large total area of muon detectors, a factor of nearly 20 larger than the CASA-MIA [4] array, resulting in a more than hundredfold improvement in background rejection.

Ground level muons become a useful separation tool for showers above \(\sim 1\) TeV at high altitude [5]. However, excellent hadron rejection power, that is, over a factor \(10^{4}\) reduction, is possible only at tens of TeV [3].

Imaging Atmospheric Cherenkov Telescope (IACT) arrays have superior rejection power to other ground-based arrays in the domain around 1 TeV, exploiting primarily the differences in shower width and substructure between electromagnetic and hadronic showers in this energy range [6,7,8]. However, IACTs have not so far demonstrated excellent rejection power at tens of TeV. Traditional separation approaches are limited by the fact that, beyond small impact distances, large events are typically not fully contained by the camera image, which significantly affects the estimation of the necessary shower parameters. The current background rejection power attained by the traditional separation methods at energies above a few tens of TeVs reaches levels between \(10^{-2}\) and \(10^{-3}\) [6, 7, 9]. This loss of performance at high energies can also be seen, albeit indirectly, in the expected background rate after background rejection cuts for CTA.Footnote 1 The surviving background rate falls by a factor 500 between 0.1 and 1 TeV but less than a factor of 10 between 10 TeV and a 100 TeV – while the proton flux falls by a factor of 50 per decade.

The ring-like images produced when ground-level muons pass through IACTs have long been used as a means of calibration [10], see [11] for a recent review. Recently, [12] suggested that the identification of a much higher fraction of muons produced in extensive air showers is possible with IACTs.

Large telescopes such as the central telescope of H.E.S.S. [13] and the Large-Sized Telescopes (LSTs) of CTA [14] enable the detection of individual muons out to large impact distance. This has traditionally been seen as a problem due to their apparent similarity to gamma-like events [15], but can also be seen as an opportunity for improvement of the background rejection power at the highest energies if characteristic differences of muons to gamma rays are identifiable.

In this paper we explore the potential for muon measurement with IACTs as a tool for background rejection by characterizing the number of muons that are detectable by a large Cherenkov telescope in proton- and gamma-initiated showers of different energies. In our simulations we adopt a hybrid approach to allow exploration of the very high background rejection power regime at acceptable cost in terms of computing time. This approach is introduced and motivated in Sect. 2. Section 3 discusses the muon content of showers from the perspective of air-Cherenkov detection and details the criteria used to label muons as detectable by a telescope. Section 4 presents the result of said criteria applied to showers initiated by both protons and gamma rays, both for large- and medium-dish telescopes. Finally, Sect. 5 discusses the implications and potential for muons as a means to improving IACT background rejection.

Table 1 Comparison of muon treatment between the CORSIKA+sim_telarray approach and the simplified muon model
Full size table

2 Cherenkov light from muons

The properties of the Cherenkov emission of a single atmospheric muon are very well determined and straight-forward to calculate in comparison to air-showers in general. The suppressed bremsstrahlung cross-section of muons with respect to electrons allows a majority of muons to reach ground-level with only ionisation losses. Similarly, the reduced multiple scattering of muons with respect to electrons means that the assumption of a linear trajectory is reasonable in most cases. For these reasons the simulation of muons in full detail may not be necessary to capture the essential characteristics of Cherenkov light from muons in showers. We implement a simplified muon model (SMM) which, starting from basic muon properties, approximates the Cherenkov light production and telescope simulation with an analytical treatment, described below. Table 1 compares the SMM approach to the combination of the CORSIKA [19] package for shower and Cherenkov light simulation and the sim_telarray [17] package for the telescope response and camera simulation. In order to verify the predictions of the SMM, we produced a small set of full CORSIKA+sim_telarray muon simulations with energies between 5 and 100 GeV for different starting heights. Note that the Cherenkov threshold for muons at 1835 m above sea level is slightly above 5 GeV.

The key parameters affecting the Cherenkov image properties of individual muons are the initial energy and production height in the atmosphere. Muons that reach ground level and land close to, or intersect, a telescope dish, produce a ring-shaped image in the telescope camera, with a full circle for muons hitting the dish and reduced sections of arc as the impact distance becomes larger. The surface brightness of the images, however, remains mostly constant, which allows muons to trigger out to large impact distances, even when the ring section captured by the camera is small enough to no longer resemble an arc, but rather a small cluster of pixels. Figure 1 illustrates the evolution of muon image properties with impact distance, as imaged in the 28 m telescope of H.E.S.S., based on CORSIKA and sim_telarray simulations.

The first step to determine whether an incoming muon can be detected by a telescope located at a certain distance from its ground impact point, is to compute the amount of Cherenkov photons collected by the telescope camera as a function of said distance. This distribution is calculated assuming a straight trajectory of the muon through the atmosphere from the production height h\(_{\mathrm {prod}}\). The atmospheric density \(\rho (h)\) and refractive index n(h) profiles are described by the same model used by sim_telarray at the H.E.S.S. location [17]. The wavelength dependence of the refractive index is ignored. For a muon with incoming zenith angle \(\theta _z\), the actual path through the atmosphere is then described by \(l = h/\cos \theta _z\) for \(h \in [h_{\mathrm {ground}}, h_{\mathrm {prod}}]\), where \(h_{\mathrm {ground}}\) is taken to be 1835 m above sea level. The emitted photons are subject to wavelength-dependent atmospheric absorption A\((\lambda ,h)\), which is integrated along the photon path, assumed here for simplicity to be the same as the muon path starting at the point where the photon is produced. The photons produced at a height h then arrive at the ground at a distance R(h) from the point where the muon hits the ground

$$\begin{aligned} R(h) = (h-h_{\mathrm {ground}})\cdot \frac{\sin \theta _c}{\cos \theta _z}, \end{aligned}$$
(1)

where \(\theta _c = \arccos ((n\cdot \beta )^{-1})\) is the Cherenkov angle of a muon traveling with velocity \(\beta = v/c\) in a medium with refractive index n. Note that \(\theta _c\) varies with height, because of energy losses and the refractive index profile.

The number of Cherenkov photons \(N_{\gamma }\) initially produced by the muon per path length dl between wavelengths \(\lambda _1\) and \(\lambda _2\), is described by

$$\begin{aligned} \frac{dN_{\gamma }}{dl} = 2\pi \alpha \left( 1-(\beta n)^{-2}\right) \int _{\lambda _1}^{\lambda _2} \frac{d\lambda }{\lambda ^2}, \end{aligned}$$
(2)

where \(\alpha \) is the fine-structure constant. This quantity needs to be convolved with the telescope response, which consists of many different elements. For this, we used the response of the telescopes in the H.E.S.S. array. The wavelength-dependent telescope response W\(_T(\lambda )\) is the combination of the mirror reflectivity, quantum efficiency of the camera and plexiglas transmittance of the camera window. Additionally, there are wavelength-independent corrections for the telescope area projection and the camera and Winston cones shadowing, which combine to a factor \(f \sim 0.6\) for a 28 m telescope. Combining the number of photons initially produced with the telescope response results in the number of Cherenkov photons produced by a muon per unit path length that are detected by the telescope:

$$\begin{aligned} \frac{dN_{\gamma , T}}{dl} = 2f\pi \alpha \left( 1{-}(\beta n)^{-2}\right) \int _{\lambda _1}^{\lambda _2} A(\lambda , h) W_T(\lambda )\frac{d\lambda }{\lambda ^2}. \end{aligned}$$
(3)
Fig. 2
figure 2

Top: Comparison of the amount of photoelectrons predicted by the analytical model for a 28 m dish (black line) and the resulting distribution for simulated muon showers for energies above 20 GeV, scattering distance smaller than 5 m and starting height \(\sim 10\) km. Bottom: Effective area comparison between the simulations and the simplified model for different starting heights. Includes all events, even those that undergo significant scattering and bremsstrahlung

Full size image

The photons distribute radially on the ground from the impact position of the muon, which defines the origin of the coordinate system. Using Eqs. 1 and 3 we can compute the ground density of detected photons

$$\begin{aligned} \rho _{\gamma , T}(R) = \frac{dN_{\gamma , T}}{dA} = \frac{dN_{\gamma , T}}{2\pi R dR} = \frac{dN_{\gamma , T}}{dl}\frac{dl}{2\pi R dR}, \end{aligned}$$
(4)

where R is the distance to the muon impact position. Placing a test telescope with a circular mirror of diameter D\(_T\) in the ground at position \((x = 0, y = R)\), the amount of photoelectrons (P.E.) collected by the dish as a function of R is then given by

$$\begin{aligned} {\mathrm {N}}_{{\mathrm {P.E.}}}(R) = \int _{R-R_T}^{R+R_T}dx\int _{-\alpha }^{\alpha }dy \cdot \rho _{\gamma , T}({\scriptstyle \sqrt{x^2 + y^2}}), \end{aligned}$$
(5)

where \(\alpha = \sqrt{R_T^2-(x-R)^2}\).

The number of photoelectrons predicted by the SMM is compared in the top panel of Fig. 2 with that resulting from the CORSIKA+sim_telarray simulations. Note that if a muon has undergone significant scattering while traveling through the atmosphere, the impact position on the ground is no longer a meaningful parameter to describe its trajectory. For this reason, the comparison shown in Fig. 2 is done with a subset of the full simulations with events selected for modest scattering. As expected, the agreement is very good. A further check utilizing the entire set of simulations is described below.

Once the properties of the photoelectrons produced by a muon and captured by the camera are known, the next step is to determine if they would activate the camera trigger. For this we base our trigger definition on the one used by FlashCam, the camera installed in the largest H.E.S.S. telescope since October 2019 [20]. The criterion is passed when an image has a group of nine neighboring pixels with a total of more than 68 photoelectrons. For each combination of muon starting height, initial energy and telescope distance from muon position, the corresponding muon image in the telescope is generated from the Cherenkov light distribution in the SMM and tested against this criterion. This allows for a trigger decision dependent on those three parameters only, which we will refer to as the simplified trigger criterion.

To verify that this is an accurate description of the real trigger conditions we once again use the small sample of CORSIKA+sim_telarray simulations. For this comparison, we include muon events undergoing significant scattering, a process which, together with bremsstrahlung, is ignored by the SMM (see Table 1). To do this, we compare, for different muon starting heights and as a function of energy and impact distance, the total number of muons arriving at the ground in the CORSIKA sample with the number of muons detected by the telescope after processing the sample with sim_telarray. Integrating the ratio of these numbers radially results in an effective area for a given energy and starting height. Similarly, for the same set of input CORSIKA events, we use their energy, production height and angular direction to compute our simplified trigger criterion. We then compute the muon effective area in the same way. These quantities are shown in the bottom panel of Fig. 2. This comparison includes muon events that undergo scattering and bremsstrahlung, confirming that the simplified trigger criterion is accurate to well within 10% across all energies.

Fig. 3
figure 3

Distribution of the arrival time of Cherenkov photons for a muon of energy 20 GeV and starting height 11 km as seen from several impact distances. The quantity in the y-axis is normalised for equal area at each impact distance

Full size image

The distribution of arrival times of the Cherenkov photons produced by the muon is straightforward to compute with the setup described above. Defining the instant in which a muon produced at height \(h_{\mu }\) arrives at the ground as \(t_{\mu }\), we can compute the relative delay in arrival time of a photon produced by the muon at height h, \(t_{\gamma }(h)\) as

$$\begin{aligned} \Delta t(h) = t_{\gamma }(h) - t_{\mu } = \frac{h\cdot c}{\cos \theta _c(h)} - t_{\mu }, \end{aligned}$$
(6)

where \(\theta _c(h)\) is the Cherenkov angle. Using Eq. 1 we compute \(\Delta t(R)\), that is, the relative photon arrival time as a function of the ground distance to the muon impact point.

Now, combining this quantity with Eq. 5 allows us to compute the number of photons that arrive to the dish of a telescope \(\text {N}_{{\mathrm{P.E.}}}\) for each value of \(\Delta t(R)\). An example of such distributions for a muon of energy 20 GeV can be seen in Fig. 3. Note that the distribution for each telescope location is normalised for equal area.

Fig. 4
figure 4

Characteristics of muons produced in proton- and gamma-ray-initiated showers. The top left panel shows the average number of muons present as a function of primary energy, with the shaded area representing the standard deviation. The remaining panels show the distributions of angular distance to the shower, energy and production height for muons produced by showers in two energy ranges, normalised for equal area. Solid lines correspond to “low-energy” (LE, 10–20 TeV) showers, while dashed lines represent muons in“high-energy” (HE, 100–120 TeV) showers

Full size image

3 Detectable muons in showers

Considering shower development only down to the Cherenkov threshold for muons allows extremely rapid simulations with CORSIKA, even up to primary energies of many hundreds of TeV. For both gamma-ray- and proton-initiated showers, over \(10^7\) showers were produced with energies between \(10^1\) and \(10^{2.5}\) TeV, distributed as \(\propto E^{-1}\) to ensure enough events at the highest energies. For all showers, the primary particle initial direction was \(\theta _z=20^\circ \). From these showers, the production height, energy, ground level direction and impact point were extracted for all the muons present. Figure 4 summarises the basic properties of these muons.

The very high statistics allows us to probe the characteristics of the very rare most muon-poor proton showers, which form the irreducible background of the muon-tagging approach. We define a muon as detectable if they fulfill three separate conditions, with the first being the simplified trigger criterion described in Sect. 2. The remaining conditions refer to the muon trajectory and are described below. We assume for simplicity that the telescope is always pointed towards the shower axis and located at the shower ground impact position. The second condition is then given by the fact that the angular distance between the muon trajectory and the axis must be smaller than half of the telescope field of view (FoV), taken here to be that of the central H.E.S.S. telescope, 3.5\(^\circ \).

Beyond simply being able to detect light from muons, in order to use their presence as a background rejection criterion, it is crucial to accurately identify them as such. Muon identification can be carried out via different techniques, which exploit the properties of the muon Cherenkov signal described in Sect. 2. A thorough study of the different muon identification techniques and their performance is beyond the scope of this paper (see [12, 21, 22] for some recent efforts). Current background rejection methods implemented in IACTs rely on the properties of the time-integrated Cherenkov light images. In this case, a complete overlap in the camera image between the main shower component and the light coming from muons makes this identification very difficult. In order to take this into account, and also to explore the effect of different muon identification efficiencies, we impose a third and final condition for detectability: a requirement on the minimum angular distance between the trajectories of the muon and the primary particle. We choose a reference value for this distance of 0.3\(^\circ \) which corresponds to the survival of, as can be seen in Fig. 5 around half of the muons that pass the first and second detectability conditions described above.

Note that possible muon identification techniques that exploit the arriving times of Cherenkov photons would not be limited by image overlap. Additionally, even techniques that are based in the integrated image can also strongly be affected by the shape of the separable muon component rather than its size in pixels. Hence this requirement on the minimum angular distance is only a proxy of the fact that only a fraction of the detected muons will realistically be tagged as such. In any case, as we show later in Sect. 4, the background rejection power attainable with the muon tagging strategy would be competitive up to a muon identification efficiency of a few percent.

Fig. 5
figure 5

Fraction of surviving muons as a function of the minimum shower separation threshold for several energies. Both the trigger and FoV criteria have already been applied. The dashed line indicates the chosen reference value of 0.3\(^\circ \) which translates to roughly 50% survival due to this criterion

Full size image
Fig. 6
figure 6

Top panel: Logarithm of the number of detectable muons in showers vs primary energy. The distribution is normalised at each primary energy bin. The type of particle initiating the showers is indicated on the top left corner. Bottom panel: Probability of 0 muons in a shower vs primary energy. Normalised together with the upper panel at each energy

Full size image
Fig. 7
figure 7

Cumulative distribution of the number of detectable muons in proton-initiated showers (solid lines) and gamma-ray-initiated showers (dashed lines) for several primary energy ranges. A round marker is placed at the x-position defined by a 70% signal efficiency. The corresponding telescope size is indicated with D\(_T\) in each panel

Full size image

4 Results

Figure 6 shows the logarithmic distribution of the number of muons classified as detectable, \(\log _{10}N_{\mu , {\mathrm {det}}}\), for showers with both protons and gamma rays as the primary respectively. The energy of the primary particle energies ranges between 10 and 250 TeV. The differences between both distributions are striking, revealing that a large fraction of the initially high (see Fig. 4) number of muons in proton showers can be in principle detectable by a large telescope. Conversely, for gamma-ray showers, the distribution is shifted towards a much lower number of detectable muons. The bottom panel of both figures shows the probability that a shower does not contain any detectable muons. Above 100 TeV this number is of the order of \(10^{-5}\) for proton showers, a number which assumes a 50% muon identification efficiency as described in Sect. 3. This number represents the irreducible background of this approach, and is orders of magnitude below the current rejection power reached at those energies [6, 7, 23], even assuming a less effective muon identification strategy.

Let us refer to the probability distributions of detectable muons shown in Fig. 6 as \(f(N_{\mu , {\mathrm {det}}} | E_{\mathrm {pri}})\), where \(E_{\mathrm {pri}}\) is the energy of the primary particle. From this quantity, we can compute the fraction of showers expected to contain less then a certain number N of detectable muons as:

$$\begin{aligned} {\mathrm {F}}(N_{\mu , {\mathrm {det}}} < N | E_{\mathrm {pri}} ) = \int _{0}^{N} f(N_{\mu , {\mathrm {det}}} | E_{\mathrm {pri}})dN_{\mu , {\mathrm {det}}}. \end{aligned}$$
(7)

Note that this corresponds to an integral as a function of the y-axis of each panel in Fig. 6. We show this quantity corresponding to several ranges of primary particle energies for both gamma-ray and proton showers in the left panel of Fig. 7. Beyond 30 TeV substantial separation power would clearly result from effective identification of single muons within shower images. For energies higher than 80 TeV the separation power could reach \(10^{-5}\), well beyond that achieved so far for IACT arrays.

4.1 Smaller telescopes

The large collection area of the telescope is critical to the number of muons whose light is collected by the dish. The simplified model results shown in Fig. 6 and the left panel of Fig. 7 correspond to a telescope with a diameter of 28 m. Repeating the same process for a smaller telescope, with a diameter of 12 m yields very different results. As can be seen in the right panel of Fig. 7, the background rejection power achievable with the muon-tagging approach worsens dramatically by several orders of magnitude. This is due to the fact that large dishes translate to a larger light collection area. Hence more light can be collected by the PMTs, allowing the detection of a higher number of muons, which are faint emitters compared to showers.

5 Discussion

It is clear from Fig. 7 that very significant potential exists for improving the background rejection of IACT arrays containing at least one large telescope at the highest energies provided that individual muon arcs can be identified with reasonable efficiency in the presence of bright shower images. A very large fraction of the muons initially present in the shower can in principle be detected by the telescope, thanks in part to the fact that the typical angular displacement from the shower is small enough to remain inside the FoV, as can be seen in Fig. 4. In order to exploit this potential, effective muon identification is required. There are several promising strategies, which use different characteristics of the muon Cherenkov light. As shown in Fig. 3, light produced by muons arrives to the telescope in short, concentrated bursts, unlike the shower photons, which span a range of tens of nanoseconds [24]. This difference could be exploited to identify muons if the time distribution of the photons, rather than only the time-integrated images, is available [25].

On the other hand, for time-integrated images, it has become very common and effective to employ template fitting approaches [26, 27] during the reconstruction. However, very energetic events produce very bright and large shower images, which are often not entirely contained in the camera field of view. The very bright pixels from the shower image likely dominate the fit and obscure the presence of dim, constant surface-brightness muon images. Simply masking out the main shower image and analyzing the residual emission will likely result in an increased rate of muon identification, provided that the image cleaning procedure is not excessively harsh. More sophisticated approaches, making use of pixel-based deep learning techniques have already shown some promising results [22]. The recent development of dedicated packages such as CTLearn [28] or GammaLearn [29] will likely facilitate the application of such techniques to the muon identification issue.

Large Cherenkov telescopes are widely seen as a useful asset for the low-energy range of the very-high energy spectrum due to their reduced threshold, but of limited use when considering the highest energies. However, as can be seen by comparing both panels of Fig. 7, this same reduced threshold is key when it comes to detecting the faint light coming from muons, which are overwhelmingly more common in proton-induced showers (see Fig. 4). Exploiting this potential of large telescopes through efficient muon identification algorithms could provide significant improvements in background rejection above several tens of TeVs, and in turn, improve the instrument sensitivity at the highest energies.

Another option is, of course, to build ground level muon detectors. However, these are not planned for either the existing or upcoming IACT arrays. Additionally, an improved muon detection technique could be applied retroactively to the entire data archive from an observatory, in the case of existing IACT arrays. This would effectively, and without increased observation time or hardware improvements, increase the detection capability at higher energies.

The results in Sect. 4 were produced for showers arriving from a distance of 20\(^\circ \) from zenith. We also produced smaller samples of showers from 0\(^\circ \) and 40\(^\circ \) to explore the effect of the zenith angle on the result. With increasing zenith angle, the number of detected muons goes up slightly, since the distance out to which muons are able to trigger goes up (see Eq. 1).

As can be seen in Fig. 6, muons are often produced, albeit in low numbers, by the highest energy gamma-ray showers. This is because at those energies, the number of interactions is so large that rare processes, such as muon pair production, become relevant. This indicates that perhaps the most useful separation parameter in terms of gamma-ray efficiency might not solely be the identification of individual muons, but rather a measure of the muonic content of an event.

The hadronic background is in fact not solely made up by protons, but also heavier nuclei [30]. The muon content in these showers is higher, which translates into a higher number of expected detectable muons. In a sense, the proton case is the worst-case scenario, being the hadronic background species most likely to initiate muon-poor showers. However, protons are most likely to masquerade as a gamma ray in traditional background rejection approaches [15], and therefore the room for improvement is greater.