Atmospheric muons at PeV energies in radio neutrino detectors

Experiments seeking to detect radio emission stemming from neutrino interactions will soon reach sensitivities that bring a detection within reach. Since experiments like RNO-G or the future IceCube-Gen2 target more than an order of magnitude more effective volume than existing experiments, the renewed and detailed study of rare backgrounds is needed. In this paper, we study the potential background from energy losses of highly energetic atmospheric muons. Due to both limited experimental measurements and limited modeling in hadronic interaction models, the expected event rate is subject to large uncertainties. Here, we estimate rate predictions and their uncertainties for different models and instrumental parameters. We also study possible routes towards mitigation of the muon background, such as parent air shower detection, and illustrate what is needed to make the first measurement of the prompt muon flux at energies above 10 PeV.


Introduction
The existence of a high-energy astrophysical neutrino flux is by now firmly established, e.g.[1].However, at energies beyond tens of PeV the flux predictions for neutrinos produced directly in sources or when ultra-high energy cosmic rays interact with cosmological photon fields vary widely e.g.[2][3][4][5][6].It is however clear that detectors with larger effective detection volumes than currently exist are necessary to discover EeV neutrinos.Radio neutrino observatories offer a promising approach to this challenge by exploiting the kilometer-scale attenuation length of radio emission in ice and the relatively low cost per detection unit [7].Among these observatories are the Radio Neutrino Observatory Greenland (RNO-G) [8], which is currently under construction, and the planned radio array for the extension of IceCube, IceCube-Gen2 [9].Both of these experiments are discovery-focused, making it essential to have a robust understanding of the signals and backgrounds involved.
The radio signal to be detected is generated by the particle cascade following a neutrino interaction in ice.The build-up of a net negative charge at the shower front leads to the emission of coherent radiation, the Askaryan emission [10].Due to a Cherenkov-like effect, the emission is strongest at the Cherenkov angle (∼56°in ice).The signal amplitude at a given observer distance scales linearly with the shower energy [11] and is typically detectable above the thermal noise at 5 PeV to 10 PeV [12].Due to this emission mechanism, any particle cascade induced in the ice with the necessary energy deposit creates a detectable signal, independent of its parent particle.This means that also high-energy muons stemming from air showers could act as a background in neutrino detectors whenever they initiate a shower [13].
Radio detectors do not need to be installed deep in the glacial ice.The antennas are typically located within 200 m below the surface, which makes them sensitive to potential anthropogenic noise 1 , as well as air shower induced backgrounds.In general, three different types of air shower backgrounds are distinguished: (1) the in-air radio emission of air showers that is refracted into the ice to the antennas; (2) the core of incompletely developed air showers can penetrate into the ice, where it induces a cascade that emits radio signals; and (3) in-ice particle showers following an energy loss of an atmospheric muon.The signatures of (1) and ( 2) have previously been studied and quantified [14][15][16][17].Both signals can be triangulated to close to the surface and therefore provide signatures that can be suppressed on an analysis level.Reflections in the ice may complicate the reconstruction, but this is true also for the neutrino signals itself.For both direct air shower backgrounds, a reasonable estimate of the background rate is possible, because the distribution of shower maxima as function of energy is relatively well-known.
The number of muon-induced background events, however, has been studied less.It has in principle been shown that muons are a non-negligible background to radio neutrino detectors in ice [13].However, the predicted event rate depends on the muon flux, which in turn strongly depends on the hadronic interaction model, and the cosmic ray composition, all of which are less-well determined, in particular at the highest energies.Furthermore, instrumental parameters, foremost the triggering system, determine the observable rate.We present a comprehensive study of the muon-induced background in this article to guide future searches for neutrinos beyond PeV energies.For energies up to 10 5 GeV, the influence of hadronic interaction models and cosmic ray spectrum on the atmospheric muon flux is studied in [18].

Predictions of muons at PeV energies and beyond
Atmospheric muons are produced in extensive air showers, which occur when high-energy cosmic rays penetrate the Earth's atmosphere.The cosmic ray nucleon interacts with an air nucleus and produces short-lived intermediate particles, mostly pions (the lightest known meson) and a few heavier particles with shorter life times, such as kaons, D-mesons, etc. Their decay gives rise to an atmospheric lepton flux, including muons.The energy range in which atmospheric muons are detectable with radio neutrino experiments is limited by the minimum muon energy that is required to produce an in-ice particle cascade with a measurable radio signal (around 10 PeV).At these energies, the flux of parent cosmic rays is low, which results in a very small muon flux.Nonetheless, this muon rate is likely comparable to the expected neutrino rate at these energies, making radio neutrino detectors the first experiments where atmospheric muons of EeV energies become relevant.While the much-discussed Muon Puzzle [19] describes a discrepancy between predicted and observed muon production in air showers for muons with energies around 1 GeV (more muons are measured than predicted by Monte Carlo simulations), the situation for muons above PeV energy is different: These muons are usually produced within the first three interactions of an air shower, rather than continuously throughout the shower development [20].The energy of a parent particle is distributed among  GeV proton induced vertical air shower according to [21].The bottom contribution is assumed to be 10% of the charm contribution [22].The photo-conversion includes muons from γ → µµ and from γ → ρ, J/Ψ → µµ.The uncertainties are a conservative estimate (constant factor in energy), taking into account experimental limitations and comparison between different event generators, they are meant rather as illustration of the current state of the field than as firm estimate, compare [23][24][25][26].
its children, which leads to lower energy particles with each further interaction in the cascade.Consequently, one has to concentrate on the highest energy interactions to study the relevant muon background.Unfortunately, these interactions are far outside of the energy regime currently observable at accelerators, which makes far-reaching extrapolations necessary.

Muon production in air showers
Atmospheric muons are produced in the hadronic cascade of an air shower mainly through the decay of short-lived mesons, namely charged pions and kaons (conventional component) [27].At very high energies the Lorentz time dilatation increases the decay length of pions and kaons to a multiple of their interaction length (ℓ int ) in air, making it more likely that they will interact and lose energy before they can decay.The contribution of particles with a shorter lifetime τ then becomes dominant as shown in Figure 1.Due to their almost immediate decay the contribution of short-lived hadrons with cτ ≪ ℓ int is called prompt flux and dominates above 10 6 GeV.Charmed hadrons (with D 0 , D + , D + s , Λ + c , Ω 0 c and their antiparticles) have large (∼10%) branching ratios into semi-leptonic modes and a lifetime τ ∼ 10 −12 s, implying a prompt decay with a probability of order 1 up to energies around 10 7 GeV [26].In principle, also bottom hadrons (B 0 , B + , B + s , B + c , Λ 0 b , Ξ 0 b , Ξ + b ), which have similar lifetimes and semileptonic decays, contribute to the prompt muon flux.B-mesons are less frequently produced by cosmic rays in the atmosphere, but their decay length is smaller, yielding a contribution to the 10 7 GeV muon flux which is 10% of the one from charm hadron decays [22].Additional significant contributions are: unflavored mesons (η, η ′ , ρ 0 , ω, ϕ) [22], photo-conversion into a muon pair (γZ → µ + µ − Z) e.g.Bethe-Heitler process, Drell-Yan processes [22] and photon conversion into a vector meson decaying into muons.These dominate the muon flux above ∼3 × 10 8 GeV [21,28].A sketch of the contributions according to [21] is shown in Figure 1.The uncertainties are a rough estimate considering experimental limitations and differences between event generators [24][25][26]29].
Taking into account these different sources, the atmospheric muon flux can then be expressed as the sum of five components: The high energy muon flux is mainly driven by the outcome of the first interaction of an air shower.The relativistic hadron-ion collisions under low momentum transfer are in the non-perturbative regime of quantum chromodynamics (QCD) [19,30], where hadron production cannot be calculated directly from first principles.Instead, effective theories and phenomenology are used.An explicit prediction of the muon flux from perturbative QCD (pQCD) is not present, but judging by the variation in different first-principle pQCD calculations predictions for the muon neutrino flux [23,24,[31][32][33][34][35] there is a sizable uncertainty.To simulate the hadron production, different (phenomenological) hadronic interaction models are used.In air shower simulations, hadronic interactions are the largest source of uncertainties, because the center-of-mass energy in the first interactions significantly exceeds the maximum energy studied at the LHC and interactions in the forward direction, i.e. high pseudorapidities are not well covered [27,36].When extrapolating to higher energies, the model predictions thus diverge even further, also with respect to the pQCD calculations.A detailed discussion of post-LHC hadronic interaction models follows in Section 2.3.
Next to the particle physics processes in the air shower, the atmospheric muon flux is determined by the cosmic ray composition.The type of the primary particle entering the atmosphere and its number of nucleons has an influence on the number of muons produced.The muon number grows less-than-linear with the primary energy of an air shower [20].This is a consequence of the energy fraction f given to charged pions in each interaction f ∼ (2/3) n , after n generations.For nuclear primaries, a nucleus with atomic number A can be treated as the sum of A separate proton air showers all starting at the same point, each with 1/A of the primary energy [20].The lower energy nucleons which initiate the shower generate fewer interaction generations, and so lose less energy to electromagnetic components [20].Therefore the number of muons is larger for heavy primaries than for showers initiated by light nuclei of the same energy.For very high-energy muons, which are created within the first interactions, this picture changes: since a proton contains the kinetic energy in one nucleon, it can produce higher energy particles than an iron primary with the same energy.Therefore, a 3 × 10 10 GeV proton shower can produce muons up to 10 10 GeV, while an iron-induced shower with the same energy and arrival direction only produces muons up to 2 × 10 8 GeV [21], as shown in Figure 2.

Muon flux simulations
For this article, the atmospheric muon flux is calculated using Matrix Cascade Equations (MCEq) [25], which describe the evolution of particle densities as they propagate through the atmosphere, using the CORSIKA parametrizations [37] as atmospheric model.The 1-dimensional cascade equations neglect the lateral versus the longitudinal development of the shower, which is important at lower energies, where the transverse momentum of the particles may be relatively important and imply a larger lateral displacement.Since this paper focuses on energies > 1 PeV this seems an acceptable limitation.Compared to computational

N per primary
Number of muons for air showers with different primaries (proton, helium, and iron) and energies (10 7 GeV and 10 10 GeV) with 60°zenith arrival direction at 3200 m (Summit Station, Greenland).Calculations were performed with Sibyll-2.3c.
extensive Monte Carlo codes like AIRES [38] and CORSIKA [37], MCEq provides a way to estimate the relative importance of a given parameter, for which accurate studies with full shower simulations would require very large statistics.

Dependence on hadronic interaction models
Several theoretical approximations describing particle production are available for different energy ranges and kinematic regimes.Different approaches have to be combined to model all hadronic interactions in air showers.For this paper, the post-LHC hadronic interaction models EPOS-LHC [39], QGSJet-II.04[40], and Sibyll-2.3c[41] are considered.While EPOS-LHC has a more general focus on minimum-bias proton-proton and heavy ion collisions, the latter two are focused on air shower simulation.
A theoretical prediction of the muon flux above PeV energies should include at least four components (see Equation 2.1), which are, however, not all taken into account in the same way in the hadronic interaction models.While the conventional flux is implemented in all models considered, only Sibyll-2.3cincludes charm production (D + , D 0 , D s , Λ c ) through a parametrization; forward charm production is intrinsically included in the nucleon PDF [19].Sibyll-2.3calso includes muons from unflavored mesons and J/Ψ.EPOS-LHC does not include charm, its prompt component arises from the decay of unflavored mesons.QGSJet-II.04only considers η decay as a production channel for prompt muons [42].The calculated muon fluxes, therefore, start to vary widely at 1 PeV where the prompt flux dominates, see Figure 3 left.EPOS-LHC and QGSJet-II.04yield the lowest muon flux because they neglect the charm component and in the latter case most unflavored mesons.The photoproduction of muon pairs that becomes relevant at PeV energies is not implemented in any model [43].Given that only Sibyll-2.3cincludes charm and unflavored mesons, it is currently the most complete model to predict the muon flux above PeV energies.However, even  Sibyll-2.3c[41], and EPOS-LHC [39] for the Global Spline Fit (GSF) [44] cosmic ray composition.Dashed lines represent the prompt muon flux and dash-dotted is conventional muon flux.The right shows different cosmic ray composition models using Sibyll-2.3c:GSF [44], HillasGaisser (H4a) [45], T+UFA, [46,47], T+H [4,46] and HillasGaisser (H3a) [45].
Sibyll-2.3cstill under-predicts the flux of muons at the highest energies due to missing production channels from photo-conversion and B-mesons, which should be addressed in future theoretical work.Given our current understanding, these components only become dominant above 3 × 10 8 GeV, so they should be a minor contribution to the background in radio detectors.The main theoretical uncertainties arise from charm cross-section calculations.The theoretical calculations are limited by the uncertainties in the scale, charm mass, and the nuclear PDFs [19].A non-perturbative intrinsic charm component may also contribute [24].This short overview illustrates the difficulty of predictions beyond LHC energies.On the one hand, theoretical uncertainties are present due to the non-availability of measurements, and on the other hand, known processes have not been implemented in all codes, since the priorities have been weighted differently for existing hadronic interaction models.We therefore use the spread between the three hadronic interaction models as indication of the current uncertainties, while keeping in mind that they do not provide the full range of possible systematic uncertainties at this point.

Dependence on cosmic ray composition
The cosmic ray spectrum covers several decades of energy up to 10 11 GeV, including particles from galactic and extra-galactic origin.Just below the so-called ankle at 8 × 10 9 GeV, the transition region from galactic to extra-galactic cosmic rays is expected [48,49], with detailed explanations still varying.
Measurements of the cosmic ray composition above a few 10 5 GeV suffer from the uncertainties in the hadronic interaction models, since the composition has to be inferred from shower parameters such as the position of the shower maximum X max , which provides composition models with much room for interpretation.Since the ultra-high energy muon flux directly depends on the cosmic ray composition, different models have been investigated to study the uncertainty stemming from this aspect.We also combine models to study the influence of galactic and extra-galactic components.This is done to show the spread in models, rather than choosing one over the other for correctness.
The well-known Hillas Gaisser models are theoretical simplifications for extreme scenarios: a heavy composition after the ankle (H3a) [45] and a proton-rich composition (H4a) [45].This is contrasted by the Global Spline Fit (GSF) [44], a data-driven parameterization that considers measurements of more than ten experiments and provides uncertainties at each energy.GSF is agnostic to theoretical models explaining the derived composition in terms of sources and propagation.
Thoudam et al. [46] published different theory-driven cosmic ray spectra up to EeV energies.In the following, their prediction for cosmic rays stemming from Supernova remnants (SNR-CR) and Wolf-Rayet stars (WR-CR) is used as a galactic component, labeled T, and are combined with different extra-galactic components.
The model by Unger, Farrar and Anchordoqui (UFA) [47] predicts a strong pure-proton component concentrated at only about one order of magnitude in energy below the ankle.For our combination into the T+UFA model the results are optimized for a pure nitrogen galactic composition, which matches the predicted composition for WR-CR [46].
The extra-galactic component of Heinze et al. [4] is based on a framework in which an ensemble of generalized ultra high energy cosmic ray accelerators is characterized by an universal spectral index (equal for all injection species), a maximal rigidity, and the normalizations for five nuclear element groups.The source evolution is included as an additional free parameter.This allows for a parameter scan with a best fit result.The composition used in this paper is obtained by a fit to the Auger data from 2019.The resulting muon flux is shown in Figure 3 for Sibyll-2.3cas hadronic interaction model.
As discussed in Section 2.1, the relevant quantity to produce high-energy muons from different primaries is the energy per nucleon.For hydrogen as a primary (with A = 1) the nucleon energy is equal to the primary energy, for heavier elements the energy scales with 1/A where A is the atomic number.On the left-hand side of Figure 4 the proton flux for the chosen models is shown, while on the right-hand side the proton fraction taking into account all nuclei, relative to their neutron number is shown.For a pure proton flux, the fraction would be 1, given that hydrogen consists only of one proton and one electron.For pure iron (with 26 protons and 30 neutrons), the fraction would be ∼ 0.46.The models start to deviate at 10 7 GeV, close to the transition region from galactic cosmic rays to extra-galactic cosmic rays.Here, the theory-based models (T+UFA, T+H) have a dip in the proton flux.The proton fraction of the GSF flux only decreases in the transition region to a fraction of 0.9 and is significantly higher than the theoretical models (around 0.5).The GSF model therefore also predicts the highest muon flux, with the exception of the proton-only scenario at energies >3 × 10 17 eV.We will therefore use mainly GSF to estimate the muon numbers going forward to remain conservative, keeping in mind that it is just one realization of the uncertainty stemming from the cosmic ray composition.

Signatures of muons in radio instruments
When a muon travels through the ice it initiates showers along its track.At PeV energies and above, the relevant shower production mechanisms are hard bremsstrahlung and pair creation events, referred to as catastrophic energy losses.As a rule of thumb, the energy of the parent particle inducing the cosmic ray air shower is roughly one decade higher than the subsequent muon.The in-ice particle cascade on the other hand typically has a shower energy one decade lower than the initiating muon.The in-ice shower energy is the important quantity for the radio emission and hence the one which determines if a muon triggers the in-ice radio detector.The Monte Carlo framework NuRadioMC [50] with its extension to simulate secondary interactions [13,51] is used to simulate the muon interaction in-ice, the subsequent Askaryan radio emission, the propagation of the radio signal to the detector, and finally the detector response to the electric field.In order to also track secondary losses of all types of leptons, the lepton propagation code PROPOSAL [52] has been included in NuRadioMC and is used for our simulations.We study the dependence on the instrument details (Section 3.1), on the muon flux itself (Section 3.2), as well as strategies to mitigate this muon background for neutrino detection (Section 4.1).
For the purpose of this work, we simulate a detector array of 35 stations, which is similar to the Radio Neutrino Observatory Greenland (RNO-G).Each station is comprised of a dipole antenna (Vpol) located at a depth of 100 m in the ice (deep component), and three log-periodic dipole antennas (LPDA) pointing straight down located at the surface (shallow component).The stations are arranged in a square grid with a spacing of 1.25 km.
Simulations are performed for several triggers to study the dependence on instrument details.The assumed noise temperature is 300 K, in both deep and shallow component.At a depth of 100 m signal-only trigger with a simple threshold of 1.5σ noise and 2.5σ noise in the band of 96 MHz to 220 MHz is evaluated.For the shallow component a high-low threshold trigger of 2.5σ noise and a two-out-of-three coincidence in the band of 80 MHz to 180 MHz is applied.The triggers for the deep component are a simplification of the phased array trigger that is the current state of the art in radio neutrino detection [53].Simulating a true phased array using a fixed trigger rate would be the best approximation of a real instrument, as done in e.g.[54].To save computing time, we chose to use the simplified trigger of a single dipole.While the 2.5σ case is likely close to the current implementation for RNO-G [12], a 1.5σ trigger is used as a proxy for potential future optimizations.A true phased array implementation will likely affect the absolute event numbers (e.g.[13]), but should not affect the relative scaling of different effects.The shallow trigger represents an optimistic performance of the current RNO-G trigger.
In order to express the detector performance, the effective area is calculated.This is done by simulating muon interactions within an ice volume containing the detector array, the initial muon position is on the air-ice planar interface.Since only the projection of the detector is perpendicular to the direction of the flux, the simulated area has to be corrected with cos(θ).The effective area (A eff ) is the projected surface area multiplied with the trigger efficiency: The expected event rate is obtained combining the effective area with an incident muon flux integrated over energy and the solid angle element of the flux, which in spherical coordinates yields an additional factor of sin(θ): (3.2)

Dependence on instrumental details
As shown in Figure 5, the shallow antennas detect the fewest muons.This is expected, as the LPDAs also have a comparatively smaller neutrino effective volume, due to their location close to the surface.As a consequence of the ice profile, in which the index of refraction increases with depth, signals propagate less often to the surface, but are bent instead towards the denser ice.The shallow LPDAs detect mostly horizontal muons above 65°zenith angle, because of the geometry constraint by the Cherenkov cone, while the deep antennas have a broad detection range with a peak around 55°zenith angle.A lower detection threshold increases the number of muons from 0.07 per year and 35 stations to 0.16 per year and 35 stations.The higher muon yield can mostly be attributed to muons in the range of 10 7 GeV to 10 8 GeV.The uncertainties shown in Figure 5 are statistical uncertainties only based on the Feldman Cousins confidence belts [55], which provide upper limits for null results and twosided confidence intervals for non-null results, which converge to a Poisson error.At lower energies, only a few geometries allow the antenna to register a signal, hence the statistics are small, and uncertainties increase due to the comparatively high muon flux at low energy.Most events (97%) are only seen in one station, regardless of the trigger configuration.

Dependence on hadronic interaction models and cosmic ray composition
The differences in the flux predictions due to the hadronic interaction models propagate almost directly into the muon rate of an in-ice radio detector.Figure 6 left shows the number of muons predicted for three different hadronic interaction models per year and 35 stations and the same cosmic ray composition.As discussed in Section 2.3, Sibyll-2.3cincludes the most production mechanisms, which explains the larger flux.In Figure 6 right, the expected muon rate for the same hadronic interaction model, but for three different cosmic ray compositions are shown.The GSF model yields the highest muon rate with a maximum between 10 7 GeV to 10 8 GeV muon energy, which is expected due to the higher proton content.

Relation to parent air shower
While the projected number of muons is relatively small, muons can still pose a problem if the neutrino rate is comparatively low.One possibility to distinguish between an atmospheric muon and a neutrino is to detect the air shower from which the muon originates.This would identify the muon and provide a veto mechanism on muon events, as also discussed in [7].

Detectability of the parent air shower
To calculate the veto efficiency, it is essential to have information about the energy and arrival direction of the air shower, as well as the distance to the nearest detector station.As highenergy muons are boosted along the air shower's axis, the cosmic ray arrival direction can . Distribution of the probability for a muon with a given energy (color code) to stem from a cosmic ray as a function of the cosmic ray energy.The x-axis indicates the energy of the inducing cosmic ray, the y-axis shows the probability that a muon with the energy indicated in the label originates from a certain cosmic ray.The relation was obtained using Sibyll-2.3cas hadronic interaction model and GSF as cosmic ray composition.While differing in detailed shape, other combinations of hadronic interaction models and cosmic ray composition provide qualitatively similar distributions.
assumed to be the same as the muon arrival direction.The location of the air shower core can be determined by projecting the muon vertex position along the arrival direction until it intersects the boundary between the ice and air.
To establish a relationship between a muon and the corresponding cosmic ray energy, Bayes' theorem can be applied.By solving the Matrix Cascade Equation with Sibyll-2.3cas hadronic interaction model for different types of primary cosmic rays (pr) -namely proton, helium, carbon, and iron -and over a range of cosmic ray energies (10 bins between 10 6 GeV to 10 11 GeV) the muon flux at ground-level can be calculated.Once the muon flux for a specific cosmic ray induced shower is known, it has to be folded with the actual flux of the primary to obtain the muon flux for all cosmic rays.Here, the number of the different primaries is drawn from the GSF cosmic ray spectrum.The probability to produce a muon with a certain energy given a cosmic ray energy p(E CR |E µ ) is calculated by The number of muons N µ is calculated for each shower, therefore it has to be summed over all possible primaries, pr.The number of cosmic rays N CR is calculated from the cosmic ray flux and also needs to be summed over all primaries.This sum is normalized by summing over all possible cosmic ray energies the muon can stem from.The distribution for different muon energies stemming from a cosmic ray with a certain energy is shown in Figure 7 for Sibyll-2.3cand GSF.The plot shows, that a muon with a given energy can stem from a variety of cosmic ray energies, most likely from a cosmic ray with an energy ∼ 10× higher than the muon energy.Since few cosmic rays have been measured above 10 11 GeV the predictions for the rate are very uncertain.Correspondingly the shape of the highest energy muon distributions vary significantly between flux models.Muons with E µ = 3×10 8 GeV (dark blue) most likely originate from proton and helium induced showers.Their probability distribution features a second peak around 2 × 10 10 GeV, which is likely associated to the cosmic ray spectrum.The GSF model provides a steep falling proton flux around 1 × 10 10 GeV and a rising helium flux with a peak around 2 × 10 10 GeV.The relation between muon and cosmic ray energy depends on the choice of hadronic interaction model and cosmic ray composition.
To calculate the veto efficiency in a RNO-G like array, for each muon event, an air shower is selected according to muon arrival direction and placed inside the array as previously described.The resulting radio signal is simulated with CORSIKA [37] and the radio extension CoREAS [56] and then folded with the detector response using NuRadioReco [57].Since the amplitude of the air shower signal scales linearly with the cosmic ray energy [58] it can now be calculated, which air shower energy is necessary to exceed a simple 2.5σ noise trigger threshold in an upward-pointing shallow LPDA antenna and hence veto the muon event.In the last step, the probability that a muon event stems from an air shower with an energy higher than the trigger threshold energy is calculated and assigned to that muon.Combined with the predicted muon flux, the number of muons that can be vetoed by detecting the parent cosmic ray can be calculated.Figure 8 shows a veto efficiency close to 100% for muon energies > 10 9 GeV.Muons originating from inclined air showers are more likely to be vetoed since the radio signal covers a larger area but becomes fainter at the same time.Therefore the veto efficiency increases with higher zenith angles only for higher energies.

Timing of air shower and muon
While the muon and the air shower stem from the same cosmic ray, the signal arrival time at the detector and subsequently at the data acquisition unit (DAQ) differs.The air shower propagates through the atmosphere with a zenith angle θ.The position where the air shower axis intersects with the ice surface is called core position, with t = 0. Here, the radio emission from the air shower is assumed to be a plane wave at the shower front, traveling the distance from the axis to the shallow antenna according to its arrival direction θ and the velocity of light in air, see Equation 4.2.The muon travels along the arrival direction of the air shower and continues into the ice until it creates a shower.From there the radio emission propagates through the ice on a bent path to the antennas.Once received by a deep antenna, the signal travels along the cable to the DAQ at the surface, see Equation 4.3.The time difference as registered in the DAQ of the radio signal stemming directly from the air shower and the subsequent muon is the difference of t µ→DAQ and t CR→DAQ with the following definitions: 3) The cable delay for 100 m coaxial cable is ∼500 ns, with c coax = 2/3c.The cable from the shallow antennas is typically 10 m, which provides a lower bound of the time difference at ∼450 ns.The full distribution is shown in Figure 9.The muon can travel up to 4 km in the ice which increases the possible travel time up to several microseconds, moreover the propagation velocity in ice is slower than in air according to the refractive index.Any air shower veto would need to take this travel time into account, by either allowing for readout with no trigger dead-time (i.e.double-buffering) or sufficiently long record lengths.Selftriggering on the air shower is challenging due to the potentially small signals and the resulting high trigger rate.A longer record length would allow a post-processing search, which simplifies background identification, however, its implementation into a low-power DAQ system is not easily possible.

Consequences for experiments
After having studied dependencies of the muon flux predictions on hadronic interaction models, composition, and instrumental details to set the stage of the uncertainties in the flux predictions, we now discuss the experimental consequences and mitigation strategies.
We will investigate whether neutrino and muon flux predictions can be treated as independent (Section 5.1), whether neutrinos and muons can be distinguished based on their experimental signature in terms of expected rates, energy or zenith distribution (Section 5.2), and whether radio detectors can be used to measure the prompt muon flux at 100 PeV energies and above (Section 5.3).

Possible connection between muon flux and neutrino flux
In Section 2.4, we established that the muon flux strongly depends on the cosmic ray composition at Earth, specifically the proton fraction, which is in turn related to the cosmic ray composition at the sources.The production of cosmogenic neutrinos is also influenced by the cosmic ray composition, as ultra-high energy cosmic rays interact with the cosmic microwave background and the extra-galactic background light [5].Moreover, the proton component plays a significant role in the generation of neutrinos, since protons produce more neutrinos than heavier nuclei when propagating through the Universe [59,60].This raises the question whether background and signal can be treated as independent from each other.
In the following analysis, we assume different cosmic ray compositions consistent with the Auger published data from 2019 and evaluate the resulting neutrino and muon events for an in-ice radio neutrino detector.We combine the galactic component by Thoudam (denoted T) with three extra-galactic components by Heinze et al. [4]: the best fit (H best fit ) with a maximal rigidity R = 1.58 × 10 9 GeV, a source evolution parameter m = 4.0, and spectral index γ = -0.7; a fit with a flat source evolution (H flat evol : R = 2.81 × 10 9 GeV; m = 0.0; γ = 0.75) and a fit with a high maximal rigidity (H high Rmax : R = 4.46 × 10 9 GeV; m = -5.6;γ = 1.6).As the fits are supposed to resemble the measured cosmic ray composition on Earth, we expect the resulting muon flux to be similar, but the large measurement uncertainties still leave room to accommodate different interpretations.The neutrino flux is calculated using the method described in [4], while the muon flux is calculated using the Matrix Cascade Equations, as detailed in Section 2.2.The result is shown in Figure 10.While the different models alter the numbers of detected muons by only a factor of two, the variation in the number of detected neutrinos changes by about a factor of ten.Since the galactic component stays unchanged, this means that only a small influence of the extra-galactic component is visible in the muon flux.The flux differs strongest at muon energies above 10 7 GeV, which is in agreement with the expected transition region from galactic to extra-galactic cosmic rays.
In other words, most muons at the relevant energies are generated by cosmic rays of 10 8 GeV to 10 9.5 GeV while cosmogenic neutrinos relevant for radio neutrino detectors stem -14 -
from cosmic rays of above 10 10 GeV.This can also be seen in the fact that the change in muon number is significantly smaller than in the neutrino number from the same models.This means that the muon background expectation can in general be treated independently from the neutrino production models.Of course, keeping in mind that some model-dependent cases are imaginable, where background and signal need to be considered together, in particular when including new physics.

Observational signatures
We now consider the practical implications for neutrino observations and analyses with a radio neutrino telescope.The observational signature for in-ice radio neutrino detectors is an electric field whose amplitude is proportional to the shower energy.The signal strength depends on the fractional energy which is deposited in the shower, so the shower energy rather than the muon or neutrino energy is the relevant observational quantity.The shower energy (which requires a reconstructed vertex distance and viewing angle, see e.g.[7,12] for details), together with the arrival direction are likely the only two reconstructed quantities that can be used to distinguish signal from background, unless a veto from air shower tagging or multiple station/pulses coincidences is possible.
The detected arrival directions of muon and neutrinos differ only slightly, as shown in Figure 11, because they are dominated by the detector geometry, which is also illustrated by the different shape of the distributions for shallow and deep component.This, however, prohibits a distinction between muons and neutrinos on an event-by-event basis and complicates it even using the whole distributions at low statistics.The only unique signature of neutrinos is an arrival direction > 90°zenith, since muons get absorbed in the Earth.However, this is only a very small fraction of the expected events.
To summarize, Figure 12 combines the most conservative and optimistic models for muon and neutrino predictions for an RNO-G like detector in terms of shower energy.In the most conservative case, an RNO-G like detector will detect 0.07 muons a year (0.16 muons with a 1.5σ trigger), and in the most optimistic case, only 0.002 muons (0.01 muons with a 1.5σ trigger).While there are thus differences in the extreme case of O(30) between the muon predictions, current neutrino flux predictions in contract vary by more than a factor of O(150).
The combination of Sibyll-2.3cas hadronic interaction model and the Global Spline Fit (GSF) yields the highest muon rate, the theoretically driven model T+H is approximately a factor two lower.QGSJet-II.04combined with T+H and the proton-poor cosmic ray composition of H3a together yields almost no muons.Recall that QGSJet-II.04does not include charm which results in an underestimation of the muon flux above PeV energies, where the prompt muon component dominates.The differences using Sibyll-2.3cand GSF, and T+H respectively are therefore likely a better estimate for the uncertainties of the muon event rate, reducing the uncertainty budget to a factor of 2, keeping in mind that Sibyll-2.3cstill does not model all components of the muon production.The neutrino flux predictions are influenced by source and propagation modeling, as well as the cosmic ray composition, as indicated by the two predictions for the composition as reported by the Pierre Auger Collaboration and the Telescope Array (TA).Without additional experimental evidence the entire neutrino parameter space has to be considered equally likely for discovery experiments.
The maxima of the muon distributions predicted in all considered scenarios are at around 10 7 GeV and fall steeply towards higher energies.Above 10 8 GeV all shown neutrino predictions are higher than the muon expectation, which provides an avenue towards a possible analysis cut at high energies.A recent study of the discovery potential for the diffuse flux of ultra-high energy cosmic neutrinos also showed the usefulness of using the reconstructed shower energy as a discriminator for the atmospheric muon background [61].
In addition, it should be noted that all showers with an energy < 10 6 GeV have their vertex position within 20 m radius of the deep antenna.While the community is pushing towards lowering the energy threshold of detectors to gain an overlap to existing (optical) experiments, the current simulations make a number of approximations which are no longer completely valid in these cases, e.g.observing the far field of the radio emission, the separation of emission and propagation, and a constant index of refraction in the emission zone.The predictions of event rates at low energies therefore carry additional uncertainties.However, Figure 12 also shows that the background problem likely becomes larger at low energies, Figure 12.Left: Expected muon event rate for a 2.5σ trigger in the deep component evaluating four extreme scenarios in combining hadronic interaction model and cosmic ray composition as stated in the label.Right: Various predictions for an expected neutrino event rate with a 2.5σ trigger in the deep component, including cosmogenic neutrinos (TA [5,62,63], Auger [5], T+H [4,46]) and neutrinos from sources (Fang and Murase [64], Rodrigues et al. [6]).
in particular since the muon flux rises much more steeply towards lower energies than the neutrino flux predictions.This is shown in a different way in Figure 13, which illustrated potential minimum energy cuts that could be imposed to gain a cleaner neutrino sample.For instance, cutting at a shower energy of 10 7.5 GeV would retain 80% or more of all expected neutrinos, but improve the signal-to-background ratio with a factor of 5−10 depending on the model.This in turn, however, raises the question how successful an extension of the detector sensitivities to lower neutrino energies can be, given the increasing muon background.

Measuring the muon flux
Finally, one can invert the approach taken above and ask whether radio detectors can be used to measure the prompt muon flux above PeV energies.As shown in Figure 8, across all energies and arrival directions, roughly 50% of the detected muons can be related to an air shower that is also detected by the same instrument, meaning that a clear identification of muon events will be possible.In the case of RNO-G, ∼ 0.3 tagged muon events are expected in 10 years at a trigger threshold of 2.5σ based on Sibyll-2.3cand GSF.Hence, the array will be too small to make a probable detection of a muon over the planned operation time.Even with an optimized trigger to 1.5σ noiseless signal equivalent, the largest flux predictions (Sibyll-2.3cand GSF) still predict < 1 tagged events in 10 years of operations.In addition, all muon signals will be very close to the threshold and thus the yet unknown analysis efficiency, as well as unstudied properties of the near-surface ice will have to be considered to solidify this number.
However, future radio detectors are already being planned, in particular, IceCube-Gen2 [9].While the precise expected event numbers will depend on the details of the detector such as the exact hardware implementation of the trigger, bandwidth of the system, and analysis efficiencies, at this point an estimate is already possible.Using the detector configuration and trigger settings as foreseen for IceCube-Gen2 [54] which includes a full simulation of . Left: Shown is the ratio of signal to background: the cumulative sum of neutrino events above the indicated shower energy in the x-axis is divided by the cumulative sum of expected muon events above the same shower energy.The color code represents the same neutrino flux models as in Figure 12.The assumed muon rate is calculated with Sibyll-2.3cand GSF as cosmic ray composition.Right: Fraction of neutrino events with shower energies higher than indicated in the x-axis.
the phased array trigger system, we simulated the atmospheric muon rates and find that IceCube-Gen2 will observe ∼ 1.9 tagged muon events in 10 years for the currently highest flux expectations of Sibyll-2.3cand GSF (see Figure 14) and ∼ 0.1 for QGSJet-II.04and H3a.With an optimized trigger, one can envision improving on these numbers to reach an expectation significantly > 0. This would allow the in-ice radio array of IceCube-Gen2 to provide the first measurements of the prompt muon flux at 10 PeV.We publish the expected muon background as a function of shower energy and incident direction for all cosmic ray composition and interaction models discussed in this paper as supplemental material (see Appendix A), so that this forecast can be incorporated in future analyses such as [61].

Conclusion and outlook
We presented a study of the background of atmospheric muons at PeV energies and beyond for radio neutrino detectors in ice.The ultra-high energy muon flux is highly dependent on hadronic interaction models and the proton fraction of cosmic ray composition.Sibyll-2.3ccurrently provides the most complete hadronic interaction model for these high energies, since it considers the conventional component, the contribution from charmed hadrons and muons from unflavored mesons, neglecting only the subdominant contribution from B-mesons and photo-conversion into muon pairs.An explicit muon flux prediction from pQCD above 10 PeV is lacking, leaving experiments with large uncertainties on their prediction of the high-energy muon rate.The main uncertainties arise from the unknown charm cross-section, which is not accessible in current particle colliders.Judging by the variations in pQCD prediction for the muon neutrino flux [32,65] there is a sizable uncertainty and more work from the theoretical side would help to establish better predictions for the muon background.The cosmic ray composition influences the muon rate mostly through the parameter of the proton fraction.Changing from a proton-rich to a proton-poor model, yields a difference of a factor of two in flux prediction.
The total observed flux is very sensitive to instrument geometry and in particular trigger settings.An RNO-G like detector will, at full completion, observe about 0.07 muons per year, using the Sibyll-2.3cprediction and a 2.5σ-threshold.At a trigger of 1.5σ this number would rise to 0.16 muons per year.These numbers should be compared to the very uncertain flux predictions for neutrinos, which are ranging from 2.7 neutrinos to 0.01 neutrinos per year in RNO-G.
Since both the neutrino and muon fluxes depend on the proton fraction of the cosmic ray composition, it was studied whether they are correlated.It could be shown, that muon and neutrino flux predictions mostly decouple.Most ultra-high energy muons stem from cosmic rays at energies lower than those that cause the cosmogenic neutrino flux.On the downside, one can therefore not reduce uncertainties through a combined treatment of signal and background.
A possible mitigation strategy is to detect cosmic rays and thereby identify muon events: if the parent air shower of the muon can be detected, it provides a signature unique to muon events.In a detector with shallow antennas, such an air shower tagging is possible directly, using the same system.The efficiency of this mechanism is energy and arrival direction dependent with good efficiency for showers with zenith angles lower than 55°and muon energies above 10 9 GeV.One could consider adding a more closely spaced array in shallowonly stations for RNO-G, which will likely improve the veto efficiency for less inclined showers.However, for high efficiency, such an in-fill array would have to have a spacing of O(100)m, making it too dense to be feasibly installed.
A discrimination between muon and neutrino signals only based on the arrival direction is unlikely, as the distributions follow mostly the detector acceptance.It is, however, likely that neutrinos and muons show a different energy spectrum.The muon flux will likely not be measurable above 10 9 GeV shower energy, already being smaller than most neutrino fluxes at 10 8 GeV shower energy.The obtainable resolution of the shower energy of radio neutrino detectors is expected to be better than a factor of two [12], which seems sufficient to assign a significant signalness probability for high energy events.Combined with an air shower veto, which is most efficient at high energies, this should allow for a relatively background-free neutrino shower detection above 10 8 GeV.
An RNO-G-like detector is likely too small to make a first measurement of the prompt muon flux at energies above 10 PeV.This could be done by using those muons that are identified as stemming from an air shower, but the expected number of these kinds of events is < 1 in 10 years.However, a much larger detector like the planned radio array of IceCube-Gen2 has the potential for the first muon measurements at these energies.Given the incompleteness of current hadronic models and lacking pQCD predictions for muons at these energies, a direct measurement will provide important additional handles, in particular on forward charm production in QCD.Measuring the muon background will also contribute to our understanding of cosmic ray composition.

Figure 1 .
Figure1.Contributions to the muon count stemming from a 10 10.5 GeV proton induced vertical air shower according to[21].The bottom contribution is assumed to be 10% of the charm contribution[22].The photo-conversion includes muons from γ → µµ and from γ → ρ, J/Ψ → µµ.The uncertainties are a conservative estimate (constant factor in energy), taking into account experimental limitations and comparison between different event generators, they are meant rather as illustration of the current state of the field than as firm estimate, compare[23][24][25][26].

Figure 4 .
Figure 4. Contribution of protons in terms of nucleon energy for the same composition models as in Figure3and indicated in the label.Left: Absolute hydrogen (proton) flux, where the nucleon energy equals the particle energy.Right: Evaluating all types of nuclei, shown is the fraction of protons relative to neutrons; The nucleon energy is 1/A of the particle energy, with A being the atomic number of the particle.

5 Figure 5 .
Figure 5. Number of triggered muons in one year and an array of 35 stations.The Σ µ denotes the sum over all bins shown.The color code indicates the trigger and its setting.Blue is a deep simple threshold trigger at 100 m with a 1.5σ threshold, violet is a 2.5σ simple threshold trigger.Yellow is a high-low threshold 2.5σ trigger with a two-out-of-three coincidence just below the surface (2 m).The assumed flux is calculated with Sibyll-2.3cand GSF as cosmic ray composition model.The shaded region indicates the 68% CL derived from the effective area calculations of the detector.Left depicts the distribution for different energies, right for different zenith arrival directions.

Figure 6 .
Figure 6.Number of triggered muons per year and an array of 35 stations with a 2.5σ noise trigger in deep component.Left: for three different hadronic interaction models[39][40][41] and GSF as cosmic ray composition model.Right: For selected cosmic ray compositions of Figure3as indicated in the label[4,44,46,47] and Sibyll-2.3cas hadronic interaction model.

Figure 8 .
Figure 8. Number of triggered muons per year and an array of 35 stations for the GSF cosmic and Sibyll-2.3cas hadronic interaction model.The trigger threshold is at 2.5σ noise in the deep (purple) and the shallow (yellow) trigger.The solid line shows the trigger muon without veto, the dashed line shows the muon number when a muon event is vetoed by the parent air shower.The air shower radio signal has to exceed a simple 2.5σ noise trigger threshold without coincidence in the shallow antennas to count as a veto.Left: with respect to the muon energy, right, with respect to the zenith arrival direction.The two lower plots show the veto efficiency for events that triggered the deep component for different zenith bands (left) and energies (right).

Figure 9 .
Figure 9. Travel time distributions for cosmic ray signals (dark purple) and muon events (light purple) traveling from the air shower core position at the ice/air interface to the data acquisition unit (DAQ) of the detector.The blue line shows the absolute time difference as registered in the DAQ of an air shower signal and the subsequent muon signal.The dotted line is the shallow cable delay relevant for cosmic rays, the solid line indicated the cable delay for muon signals, assuming 100 m coaxial cable.