Muons as a tool for background rejection in Imaging Atmospheric Cherenkov Telescope arrays

The presence of muons in air-showers initiated by cosmic ray protons and nuclei is well established as a powerful tool to separate such showers from those initiated by gamma rays. However, so far this approach has been fully exploited only for ground level particle detecting arrays. We explore the feasibility of using Cherenkov light from muons as a background rejection tool for imaging atmospheric Cherenkov telescope arrays at the highest energies. We adopt an analytical model of the Cherenkov light from individual muons to allow rapid simulation of a large number of showers in a hybrid mode. This allows us to explore the very high background rejection power regime at acceptable cost in terms of computing time. We show that for very large ($\gtrsim$20 m mirror diameter) telescopes, efficient identification of muon light can potentially lead to background rejection levels up to 10$^{-5}$ whilst retaining high efficiency for gamma rays. While many challenges remain in the effective exploitation of the muon Cherenkov light in the data analysis for imaging Cherenkov telescope arrays, our study indicates that for arrays containing at least one large telescope, this is a very worthwhile endeavor.


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]). a e-mail: Laura.Olivera-Nieto@mpi-hd.mpg.de 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 ∼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 groundbased 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 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 E>20 GeV 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 Section 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 tele-scope. 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, Section 5 discusses the implications and potential for muons as a means to improving IACT background rejection.

Cherenkov light from muons
The properties of the Cherenkov emission of a single atmospheric muon are very well determined and straightforward 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. of the SMM, we produced a small set of full COR-SIKA+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 prod . The atmospheric density ρ(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 θ z , the actual path through the atmosphere is then described by l = h/ cos θ z for h ∈ [h ground , h prod ], where h ground is taken to be 1835 m above sea level. The emitted photons are subject to wavelength-dependent atmospheric absorption A(λ, 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 where θ c = arccos((n · β) −1 ) is the Cherenkov angle of a muon traveling with velocity β = v/c in a medium with refractive index n. Note that θ c varies with height, because of energy losses and the refractive index profile. The number of Cherenkov photons N γ initially produced by the muon per path length dl between wavelengths λ 1 and λ 2 , is described by where α 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 (λ) is the combination of the mirror reflectivity, quantum efficiency of the camera and plexiglas transmittance of the camera window. Additionally, there are wavelengthindependent corrections for the telescope area projection and the camera and Winston cones shadowing, which combine to a factor f ∼ 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: The photons distribute radially on the ground from the impact position of the muon, which defines the origin of the coordinate system. Using Equations 1 and 3 we can compute the ground density of detected photons 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 The number of photoelectrons predicted by the SMM is compared in the top panel of Figure 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 Figure 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 Figure 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.
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 µ arrives at the ground as t µ , we can compute the relative delay in arrival time of a photon produced by the muon at height h, t γ (h) as where θ c (h) is the Cherenkov angle. Using Equation 1 we compute ∆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 Equation 5 allows us to compute the number of photons that arrive to the dish of a telescope N P.E. for each value of ∆t(R). An example of such distributions for a muon of energy 20 GeV can be seen in Figure 4. Note that the distribution for each telescope location is normalised for equal area.

Detectable muons in showers
Considering shower development only down to the Cherenkov threshold for muons allows extremely rapid sim- ulations 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 ∝ E −1 to ensure enough events at the highest energies. For all showers, the primary particle initial direction was θ z = 20 • . From these showers, the production height, energy, ground level direction and impact point were extracted for all the muons present. Figure 3 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 muontagging approach. We define a muon as detectable if they fulfill three separate conditions, with the first being the simplified trigger criterion described in Section 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. 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 • .
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 Section 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 • which corresponds to the survival of, as can be seen in Figure 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 Section 4, the background rejection power attainable with the muon tagging strategy would be competitive up to a muon identification efficiency of a few percent.  Figure 6 shows the logarithmic distribution of the number of muons classified as detectable, log 10 N µ,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 Figure 3) number of muons in proton showers can be in principle detectable by a large telescope. Conversely, for gammaray 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 Section 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 Figure 6 as f (N µ,det |E pri ), where E 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:

Results
Note that this corresponds to an integral as a function of the y-axis of each panel in Figure 6. We show this quantity corresponding to several ranges of primary  Cumulative distribution of the number of detectable muons in proton-initiated showers (solid lines) and gamma-rayinitiated 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.
particle energies for both gamma-ray and proton showers in the left panel of Figure 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.

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 Figures 6 and the left panel of Figure 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 Figure 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.

Discussion
It is clear from Figure 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 Figure 3. 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 Figure 4, 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 Figure 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 Figure 3). 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 Section 4 were produced for showers arriving from a distance of 20 • from zenith. We also produced smaller samples of showers from 0 • and 40 • 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 Equation 1).
As can be seen in Figure 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.

Declarations
Funding This research was supported by the Max Planck Society and ETH Zurich. AMWM is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project Number 452934793, MI 2787/1-1.

Conflicts of interest Not applicable.
Data Availability Statement Data sharing not applicable to this article as no datasets were generated or analysed during the current study. Code availability Not applicable.