Enhancing the cosmic-ray mass sensitivity of air-shower arrays by combining radio and muon detectors

Abstract

The muonic and electromagnetic components of air showers are sensitive to the mass of the primary cosmic particle. The sizes of the components can be measured with particle detectors on ground, and the electromagnetic component in addition indirectly via its radio emission in the atmosphere. The electromagnetic particles do not reach the ground for very inclined showers. On the contrary, the atmosphere is transparent for the radio emission and its footprint on ground increases with the zenith angle. Therefore, the radio technique offers a reliable detection over the full range of zenith angles, and in particular for inclined showers. In this work, the mass sensitivity of a combination of the radio emission with the muons is investigated in a case study for the site of the Pierre Auger Observatory using CORSIKA Monte Carlo simulations of showers in the EeV energy range. It is shown, that the radio-muon combination features superior mass separation power in particular for inclined showers, when compared to established mass observables such as a combination of muons and electrons or the shower maximum $X_{\mathrm{max}}$. Accurate measurements of the energy-dependent mass composition of ultra-high energy cosmic rays are essential to understand their still unknown origin. Thus, the combination of muon and radio detectors can enhance the scientific performance of future air-shower arrays and offers a promising upgrade option for existing arrays.

1 Introduction

More than 100 years after the discovery of cosmic rays, the origin of the ultra-high energy cosmic rays is still under investigations. Recently, the Pierre Auger Collaboration found evidence for an extragalactic origin of cosmic rays above $8 \times 10^{18}$ eV by measuring an anisotropic distribution of the arrival directions [1]. However, the sources accelerating these cosmic rays up to the highest energies are still unknown. To answer open questions about their sources and propagation mechanisms through studies of the arrival directions and the energy spectrum, accurate measurements of the mass composition are essential. Cosmic rays above $10^{14}$ eV are measured indirectly via extensive air showers in the Earth’s atmosphere. Thereby, the primary cosmic ray induces a cascade of secondary particles, of which mainly muons, electrons and photons arrive on ground and can be measured with particle detectors. In addition, mainly the electrons induce fluorescence light, Cherenkov light, and radio emission in the atmosphere, which can be measured on ground [2,3,4]. All three methods have comparable intrinsic sensitivity to the electromagnetic shower component, but optical detectors are restricted to clear and dark nights while radio antennas can be operated at almost any weather conditions. In the last few years radio arrays have achieved a measurement accuracy competitive to the optical techniques. This stimulated this research on possible synergies between radio and particle detectors, i.e., the two techniques allowing for continuous operation.

The development of the muonic and electromagnetic components in a cosmic-ray air shower depends on the mass of the primary particle. Since air showers induced by heavier particles develop faster, the mass of the primary particle can be estimated statistically in two ways. First, the atmospheric depth of the shower maximum, $X_{\mathrm{max}}$, is on average smaller for heavier particles. $X_{\mathrm{max}}$ can be measured best by optical and radio detectors, where the most accurate measurements of energy and $X_{\mathrm{max}}$ are currently provided by fluorescence telescopes [5], because the systematic uncertainties are well understood. Nevertheless, recent radio arrays have already reached similar accuracy [6,7,8]. Second, more muons and less electrons are produced compared to showers of light primary particles. This can be used to estimate the primary mass by measuring the muonic and electromagnetic components of the same air shower separately, since it influences the ratio between the numbers of muons and electrons at the shower maximum as well as on ground.

An established method is to measure the number densities of electrons and muons at a reference distance to the shower axis on ground (e.g. in CASA-MIA [9, 10], AGASA [11], KASCADE-Grande [12], and AugerPrime, the Upgrade of the Pierre Auger Observatory [13]). Except for showers more inclined than 40$^{\circ }$, the muons rarely interact or decay in the atmosphere and their number is approximately constant from the shower maximum to the ground. However, the electrons are partly absorbed in the atmosphere and suffer larger energy losses so that their number depends on the distance to the shower maximum. Especially for showers at large zenith angles, the distance to the shower maximum is long and the number of electrons falls below the detection threshold. On the contrary, the radio emission is produced by the electrons and positrons along the shower and is not absorbed in the atmosphere. Thus, the radio signal provides information about the size of the electromagnetic component for all zenith angles. Furthermore, the width of the radio footprint on ground rises with the zenith angle [14], which enables the economic detection of inclined showers with radio antenna arrays with a large spacing.

We studied the combination of radio and muon detection for the site of the Pierre Auger Observatory [15] in Malargüe, Argentina, where the combination of muon and radio detectors is already realized. Its main detector is dedicated to measure cosmic rays above $10^{18.5}$ eV and comprises a 3000 km$^2$ large surface array of water-Cherenkov detectors overlooked by fluorescence detectors at four sites [16]. In the enhancements area of the Observatory, the AMIGA (Auger Muons and Infill for the Ground Array) [17] and AERA (Auger Engineering Radio Array) [18] detectors measure showers in coincidence down to energies of $10^{17.5}$ eV. AMIGA consists of water-Cherenkov detectors on an area of 23 km$^2$, arranged on a dense grid of 750 m spacing compared to 1500 m in the standard array, to lower the energy threshold. Below the water-Cherenkov detectors, underground scintillators (AMIGA Muon Detector) are being installed at 2.3 m depth as part of the AugerPrime Upgrade. The water-Cherenkov detectors are sensitive to all shower particles whereas the underground scintillators are shielded from the electromagnetic component of the shower. Thus, the underground detectors solely measure muons with an energy threshold of about 1 GeV. AERA comprises antenna stations distributed on an area of 17 km$^2$ inside AMIGA to detect the radio emission of cosmic-ray air showers in the frequency band of 30–80 MHz. By measuring the radio emission, AERA is mainly sensitive to the charged electromagnetic component of the shower. Motivated, by the existing AMIGA and AERA arrays, we study for the same simulated air showers muons above 1 GeV and the radiation energy in the band of 30 – 80 MHz.

In this work, the capability of the radio emission is evaluated to serve as a mass estimator in combination with the muons reaching ground. To study the pure shower physics, air-shower simulations mostly independent of detector properties are used. Apart from the energy threshold for muons and the radio frequency band, we do not simulate any specific detector response or array spacing, but instead study particle numbers and densities and the total radiation energies. To facilitate the application to the real detectors, we use the geomagnetic field and height above sea level of the Auger site of 1552 m a.s.l. Different parameters such as the particle numbers and radiation energy are studied for proton- and iron-induced air showers. The mass separation power is investigated depending on the zenith angle and the primary energy. It shows that in particular for inclined showers using the radio emission is superior compared to classical detection methods using solely particles on ground.

2 Air-shower simulations

The air-shower simulations used in this work are calculated with the simulation code CORSIKA [19], using the hadronic interaction model QGSJETII-04 [20]. The true particle distributions calculated from CORSIKA are used to investigate the muon density at different distances to the shower axis and the radiation energy. To facilitate the foreseen application, the simulations used in this work were prepared in the energy range of the AMIGA Muon Detector and AERA. Thereby two simulation sets are used:

• 1.

  • fixed zenith angle of $38^{\circ }$

  • energy range $10^{17.5}{-}10^{19}$ eV, isotropically distributed

  • azimuth 0$^{\circ }$–360$^{\circ }$, isotropically distributed

  • 1000 simulations

• 2. distribution of the direction and energy according to AERA measurements:

  • energy range of $2 \times 10^{16}{-}4 \times 10^{19}$ eV

  • full zenith and azimuth angle range

  • 10,604 simulations.

Both sets contain an equal number of proton and iron primaries. The air showers are simulated until an observation level of 1552 m a.s.l., corresponding to 870 g cm$^{-2}$ atmospheric overburden and according to the average altitude of the AMIGA and AERA detectors. The Earth’s magnetic field is simulated according to the Auger site. For all results, only muons above an energy of 1 GeV are considered, corresponding to the energy threshold of the AMIGA Muon Detector at 2.3 m underground, i.e. an additional overburden of 540 g cm$^{-2}$ due to the soil [17, 21]. In order to study the general potential of combining muon and radio detection, we do not apply the specific detector responses or array layouts of AMIGA or AERA, but work with the true observables from the CORSIKA simulations.

3 Definitions and methods

In this work, the mass separation power is evaluated using the figure of merit (FOM). The figure of merit quantifies the separation of two classes for an observable that has different mean values for each class, e.g, the separation of proton and iron showers by the muon number. The FOM is defined as the following relation between the mean value and the standard deviation of the observable for the two classes:

$$\begin{aligned} FOM = \frac{ | \mu _{Fe} - \mu _p | }{ \sqrt{ \sigma _p^2 + \sigma _{Fe}^2 } }, \end{aligned}$$

(1)

where $\mu _{\text {i}}$ are the mean values and $\sigma _i$ the standard deviations of the proton (p) and iron (Fe) classes. The figure of merit is a valid statistical estimator for Gaussian distributions. For example, a $FOM = 1$ indicates that the means of the two classes are separated by one standard deviation of the difference between these observables. In particular, we use the FOM to quantify the separation between proton and iron showers based on the true observables derived from the CORSIKA simulations.

To calculate the mean and the standard deviation of the investigated observables in this study, they are normalized to an energy of $10^{18}$ eV for all simulated air showers. The energy dependencies are derived from power-law fits to the complete set of simulations:

$$\begin{aligned} a \left( E \right) = a_0 \cdot \left( \frac{E}{10^{18} ~\text {eV}} \right) ^{\gamma }, \end{aligned}$$

(2)

where a(E) is any energy-dependent observable, $a_0$ its value at $10^{18}$ eV, E the primary energy used in the simulation, and $\gamma$ the index of the power law used to approximate the energy dependence.

Fig. 1

Zenith angle dependence of the mean true number of muons at observation level ($\hbox {E}_{\mu } > 1$ GeV). The number of muons is normalized to an energy of $10^{18}$ eV. For zenith angles larger than 40$^{\circ }$ the traveled distance becomes large and a fraction of the muons decays during the shower development before reaching the observation level of 1552 m a.s.l. As also in the following figures, the lines and their surrounding bands denote the mean values and the standard deviations

Fig. 2

True muon density at different distances to the shower axis above a muon energy of 1 GeV for showers induced by a proton (a) and an iron nucleus (b) with a zenith angle of 38$^{\circ }$. The muon density is directly correlated to the total number of muons on ground and shows approximately the same energy dependence

4 The muonic component

Muons interact much less with matter, i.e. the atmosphere, than electrons and have only negligible energy losses from bremsstrahlung due to their larger mass. The muon is an unstable particle, but its decay is only mediated by the weak interaction and therefore slow with a mean lifetime of 2.2 $\upmu$s. Hence, this decay only becomes important for inclined showers, where the distance from $X_{\mathrm{max}}$ to the observation level is of the same order of the distance traveled in the (time dilated) lifetime. In Ref. [21] it is shown, that the number of muons at $X_{\mathrm{max}}$ and on ground deviates for zenith angles above 40$^{\circ }$. Moreover, the total number of muons of a shower when reaching the ground depends strongly on the primary energy. A power-law fit to the simulations with the full zenith angle range shows that the muon number increases with the energy with a mean index of $\gamma = 0.93$ (for proton $\gamma = 0.922\, \pm \,0.004$ and for iron $\gamma = 0.934 \pm 0.004$). This mean index will be used in the following to correct for the energy dependence of the number of muons $N_\mu$. For comparison, older versions of the different hadronic models predict similar values between $\gamma = 0.88 {-} 0.92$ [22].

In Fig. 1 the mean true number of muons $N_\mu$ at an observation level of 1552 m a.s.l. is plotted over the zenith angle. For each shower the number of muons is normalized to the corresponding value at an energy of $10^{18}$ eV, using the obtained fit results of the energy dependence. For zenith angles up to around 40$^{\circ }$ the number of muons is stable up to the observation level. For larger zenith angles a growing fraction of the muons decays before reaching the observation level.

The number of muons at the observation level is larger for iron than for proton showers for all zenith angles. The spread (standard deviation) is larger for proton showers due to larger shower-to-shower fluctuations. However, the separation between proton and iron is larger than the spread, which shows the potential of the number of muons on ground for estimation of the primary mass, provided that the primary energy is known.

4.1 Observable: muon density at a reference distance

The particle footprint of a cosmic-ray air shower extends over several square kilometers on ground at the energies investigated here. Therefore, realistic experiments cannot measure the total number of muons directly, but instead locally measure the number of muons at several positions with sparse detector arrays. This corresponds to a measurement of the muon density $\rho _{\mu }$ ($= \hbox {number of muons per area}$) at certain distances from the shower axis. Thus, we calculated the muon density at certain distances from the shower axis from the muon output of the CORSIKA simulations. As shown in Fig. 2, the muon density at a chosen distance is directly correlated to the total number of muons at the observation level and can be used as an observable for the latter. The figure shows the true muon density for proton (Fig. 2a) and iron (Fig. 2b) primaries at distances from 300–1000 m as well as the true total number of muons at an observation level of 1552 m a.s.l. compared to the primary energy for showers with a zenith angle of 38$^{\circ }$. The muon density decreases with the distance to the shower axis. Furthermore, the muon density shows a slightly lower energy dependence (smaller index $\gamma$) than the total number of muons for all distances, e.g., for a distance of 600 m in the simulations including all zenith angles the mean index is $\gamma = 0.90$ (for proton $\gamma = 0.894 \pm 0.002$ and for iron $\gamma = 0.907 \pm 0.002$). As explained below, we have selected 600 m as reference distance because of the high mass separation power of the muon density at 600 m. Hence, in the following the muon density is normalized to a primary energy of $10^{18}$ eV with an index of $\gamma = 0.90$.

Fig. 3

True muon density compared to the distance to the shower axis for showers with a zenith angle of 38$^{\circ }$, normalized to $10^{18}$ eV. The muon density decreases with the distance. Thereby, the relative difference between proton- and iron-induced showers increases, whereas the relative spread is constant

The mean muon density is shown over the distances of 300–1000 m from the shower axis in Fig. 3 for showers with a zenith angle of 38$^{\circ }$, normalized to a primary energy of $10^{18}$ eV. The muon density decreases with the distance to the shower axis for both proton and iron showers. The mass separation power depends on the relative difference and the spread at each distance, which is quantified by the figure or merit shown for different ranges of zenith angles in Fig. 4. In addition, the figure of merit is shown for the zenith angle range of 0$^{\circ }$–55$^{\circ }$, at which AMIGA features full detection efficiency [23]. For zenith angles below 20$^{\circ }$, the figure of merit slightly decreases with the distance. For showers more inclined than 40$^{\circ }$, it increases with the distance and does not reach a maximum until 1000 m. Combining the zenith angles from 0$^{\circ }$ to 55$^{\circ }$ leads to a broad maximum in the figure of merit between 600 and 800 m for a muon energy threshold of 1 GeV. Since we expect real experiments to suffer less from measurement uncertainties at higher values of $\rho _\mu$, we have selected 600 m as reference distance for the present study.

Fig. 4

Figure of merit compared to the distance to the shower axis for different zenith angle ranges and for combining showers at all zenith angles. The muon density is normalized to a primary energy of $10^{18}$ eV for all showers. The figure of merit decreases with the zenith angle due to larger shower-to-shower fluctuations. It shows a maximum between 600–800 m for the range of 0$^{\circ }$–55$^{\circ }$

Fig. 5

Zenith angle dependence of the muon density at 600 m above a muon energy of 1 GeV. The muon density increases slightly for zenith angles up to 50$^{\circ }$, which indicates that up to these zenith angles the production of high energy muons dominates. For larger zenith angles it decreases due to muon decay

The dependence of the muon density $\rho _{\mu }^{600}$ on the zenith angle is shown in Fig. 5. It slightly increases for zenith angles up to around 50$^{\circ }$ and decreases again at higher zenith angles due to muon decay.

5 The electromagnetic component and the radio emission

The number of electrons $N_e$ in an air shower can be measured by particle detectors on ground in a similar way as the muons. However, the electrons are partly absorbed in the atmosphere on their way to the ground and suffer much larger energy losses, e.g. by bremsstrahlung, than the much heavier muons. Therefore, their number in the shower strongly depends on the distance in atmospheric depth to $X_{\mathrm{max}}$, which increases with the zenith angle $\theta$ roughly with $1/\cos {\theta }$. The number of electrons at an observation level of 1552 m a.s.l. and at the electromagnetic $X_{\mathrm{max}}$ are plotted over the zenith angle in Fig. 6. As expected, the number at the observation level decreases with the zenith angle and is about three orders of magnitude smaller at an angle of 80$^{\circ }$ compared to vertical showers. This dependence on the zenith angle has to be taken into account when using the electron number for mass separation measurements, which leads to additional uncertainties from the measurement of the arrival direction of the shower. Depending on the size and type of the particle detectors and on the observation level, the number of electrons falls below the detection threshold and only the muons are detected for very inclined showers.

Fig. 6

Zenith angle dependence of the true number of electrons above an energy of 250 keV. The number of electrons at 1552 m a.s.l. decreases by about three orders of magnitude over the plotted zenith angle range due to absorption in the atmosphere. Proton showers contain more electrons than iron showers except for very inclined showers, where the electrons are mainly products of muon decays. On the contrary, the number of electrons at $X_{\mathrm{max}}$ is nearly independent of the zenith angle and higher in proton showers

Moreover, it becomes apparent that the difference between proton and iron showers becomes smaller and finally flips the direction, so that for $\theta>$ 65$^{\circ }$ iron showers contain more electrons at the observation level than proton showers do. The shower-to-shower fluctuations are increased around the zenith angle of the flip. In addition, the slope becomes smaller towards higher zenith angles. These observations are explained by the increasing number of muon which decay into electrons (cf. Fig. 1). Hence, these electrons are mostly created by the muonic component. Thus, for inclined showers, the number of electrons is correlated with the number of muons in the shower, which is larger for heavier primary particles. Due to this flip in the proton-iron separation and the strong decrease, the number of electrons at the observation level does not provide a reliable mass estimator for inclined air showers.

On the contrary, the number of electrons at $X_{\mathrm{max}}$ does not depend on the zenith angle. Independent of the shower inclination, it is larger for proton than for iron-induced air showers and, thus, provides information about the mass of the primary particle. However, the number of electrons at $X_{\mathrm{max}}$ cannot be measured directly by air-shower arrays on ground, but indirectly by the electromagnetic energy deposited in the atmosphere. This electromagnetic shower energy is slightly larger for proton than for iron showers, e.g., by 4.5% for a primary energy of $10^{18}$ eV and 3% at $10^{19}$ eV [24]. The electromagnetic energy of an air shower can be measured by its radio emission and in dark clear nights additionally by the fluorescence light and the Cherenkov light produced in the atmosphere. Hence, in contrast to direct measurements of the number of electrons on ground, indirect measurements of the number of electrons at $X_{\mathrm{max}}$ by radio (or optical) detectors provide a useful observable for the combination with muon measurements over all zenith angles including very inclined showers.

5.1 Observable: the radiation energy of the radio emission

The radio emission of an air-shower is induced by the electromagnetic particles in the shower [3, 4]. Hence, the energy contained in the radio emission – the radiation energy $E_{\text {rad}}$ – provides a calorimetric measurement of the electromagnetic shower energy $E_{\text {em}}$. This correlation between the radiation energy and the electromagnetic shower energy was observed at the Pierre Auger Observatory [25]. It was modeled and corrected for various dependencies on the arrival direction in Ref. [26], which is used here to calculate a corrected radiation energy $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ from the shower simulations (see Appendix A). Thereby, the radio emission in the frequency band of 30–80 MHz is considered.

Fig. 7

True radiation energy compared to the zenith angle. The radiation energy is corrected for the angle between the geomagnetic field and the shower axis as well as for the air density at the mean $X_{\mathrm{max}}$ at the corresponding zenith angle. Furthermore, the radiation energy is calculated for the observation level of 1552 m a.s.l. The radiation energy increases with the zenith angle up to 50$^{\circ }$, whereas this effect is larger for proton showers. This leads to an enhanced mass sensitivity at this zenith angle

The corrected radiation energy $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ is plotted over the zenith angle in Fig. 7 after normalization to $10^{18}$ eV by assuming $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ $\propto E^2$ according to Ref. [25]. $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ grows with the zenith angle, as inclined showers extend over a larger geometric distance on which the radiation energy is released. In addition, the shower-to-shower fluctuations decrease with increasing zenith angle. Proton showers release more radiation energy than iron showers due to the larger amount of electrons and positrons at $X_{\mathrm{max}}$. In fact, the difference between proton and iron showers grows with the zenith angle. This is expected, since the mean free path of the shower particles grows and, thus, the difference between the total path lengths of all electrons and positrons. In addition, the full development of an iron shower is shorter in atmospheric depth than of a proton shower with the same primary energy. This becomes apparent, when the iron shower ($\hbox {A}=56$) is described as 56 parallel developing “proton” sub-showers with each a primary energy of E₀/56 as in Ref. [27]. Each of these sub-showers develops faster than the proton shower with the primary energy E₀ and hence the whole iron shower does. Therefore, the proton shower travels a longer geometric distance on which more radiation energy is released. This effect becomes larger for more inclined showers, where the ratio between the atmospheric depth and the geometric distance becomes larger.

6 Mass estimation by combining observables

Combining the two observables, the muon density and the radiation energy, leads to a mass sensitive parameter introduced in this section. Summarizing the results for the electromagnetic component, the radiation energy shows an enlarged mass sensitivity at higher zenith angles, reaching a plateau at around 50$^{\circ }$. In contrast, the number of electrons looses its mass sensitivity at angles above 60$^{\circ }$ due to the fact that mainly electrons originating from muon decay reach the observation level. In addition, the uncertainties due to shower-to-shower fluctuations are particularly large for inclined showers. The muonic component shows a different behavior. The difference between proton and iron showers for the muon density at 600 m only decreases slowly for large zenith angles.

Fig. 8

Correlation between the muon density and the radiation energy and a power-law fit

Fig. 9

The correlation between the muon density and the radiation energy (as in Fig. 8) using different hadronic models. The average index $\gamma$ of a power-law fit to the proton and iron distributions shows only minor differences. Therefore, only QGSJETII-04 is used for all following investigations

The quantities of the electromagnetic component feature higher values for proton showers ($N_e$ only up to 60$^{\circ }$ zenith angle), the muon density is higher for iron. Therefore, it is expected that the mass sensitivity is enhanced by combining these complementary observables. The muon density is plotted over the radiation energy in Fig. 8 for QGSJETII-04 simulations with a zenith angle of 38$^{\circ }$. The simulations are repeated using different hadronic interaction models, i.e. Sibyll 2.3 [28] and EPOS-LHC [29], for comparison in Fig. 9 (shown as power-law fits to the simulations). Sibyll 2.3 predicts more muons for proton showers at small radiation energies and EPOS-LHC at higher radiation energies. The models mainly differ in the absolute scale of the predicted muon density. The absolute scale is relevant for the interpretation of a mass estimator in terms of atomic mass numbers, but not for the separation between light and heavy primaries investigated here. The difference in the indices of the power law used to normalize the energy dependence (cf. Sect. 3) is statistically significant, but small in size. Therefore, we have not examined the detailed impact of using different hadronic interaction models in the context of this conceptual study, and performed the full analysis using QGSJETII-04.

Fig. 10

Ratio between the muon density and the square root of the radiation energy compared to the zenith angle. The separation of proton and iron showers exceeds the spread originating from shower-to-shower fluctuations. The ratio decreases for larger zenith angles, since the muon density decreases (cf. Fig. 1)

In Fig. 10 the ratio between the muon density at 600 m axis distance and the square root of the radiation energy is plotted over the zenith angle. Fitting the ratio over the primary energy with a power law for the simulations over all zenith angles leads to a correlation of $\rho _{\mu }^{600} / \sqrt{{\mathrm{S}}_{{\mathrm{RD}}}^{\rho _{\theta }}}$ $\propto$ E$^{0.058}$, which is used here to normalize the ratio to $10^{18}$ eV. Since both, the muon density at 600 m and the radiation energy, slightly increase with zenith angle until about 50$^{\circ }$, the ratio is nearly independent of the shower inclination until 50$^{\circ }$ zenith angle. It decreases at zenith angles above 50$^{\circ }$, since the muon density decreases (see Fig. 1). The separation between proton and iron showers is larger than the spreads (standard deviations) of both distributions for all investigated zenith angles.

Fig. 11

Figure of merit of the different shower observables. The different mass estimators are corrected for their dependence on the true primary energy, which is known as input parameter of each simulation. In addition, the uncorrected ratios are shown providing a more realistic estimate of the potential for real air-shower arrays. The bands depict the uncertainties due to shower-to-shower fluctuations. All observables are true values derived from CORSIKA simulations. They do neither include the effects of a specific layout of an air-shower array nor detector specific uncertainties

The mass separation power represented by the figure of merit is shown in Fig. 11 for the ratio of the muon density and the radiation energy, the ratio of the muon density and the number of electrons at 1552 m a.s.l., the muon density, and $X_{\mathrm{max}}$. The figure of merit of $X_{\mathrm{max}}$ serves as a reference, since $X_{\mathrm{max}}$ is a widely used observable for measurements of the mass composition [2]. The $X_{\mathrm{max}}$ values are directly obtained from the same CORSIKA simulations as the other observables, not including any assumptions on a specific detection technique. All observables are shown with and without normalization for their dependencies on the true energy of the primary particle, which is known exactly in the simulations. The muon density alone has only mass separation power if the primary energy is known, and therefore has no counterpart without normalization. In a realistic experiment, the primary energy might not be reconstructed accurately enough for the purpose of normalization. Therefore, we also show the figure of merit for the unnormalized observables, so the potential improvement by normalization for the primary energy can be assessed by comparing both sets of curves. Whether or not normalizing for the primary energy, the classical method of the muon-electron ratio looses mass sensitivity with increasing shower inclination. In contrast, both $X_{\mathrm{max}}$ and the new combination of the muon density and the radiation energy can be used for the purpose of mass separation for the full range of investigated zenith angles.

7 Discussion

The results shown here compare the well established mass estimator using the particles numbers with the novel method combining the muons with the radio emission. The first proof of principle conducted here illustrates the potential of radio-emission measurements to enhance mass estimation in particular for inclined showers. Potentially, the mass sensitivity can be improved further by investigating other radio observables, such as the radio amplitude or energy fluence at a reference distance instead of the integral $\sqrt{S_{{\mathrm{RD}}}^{\rho _{\theta }}}$ over the whole footprint, and by tuning the reference distance for the muon density as a function of zenith angle instead of using a fixed distance of 600 m. Furthermore, the method of mass estimation via the ratio of the muon density and the radiation energy can be combined with the independent mass estimator $X_{\mathrm{max}}$, measurable by radio as well as optical detectors, to reduce the overall uncertainties on the mass.

Figure 11 shows the intrinsic mass sensitivities of various observables not including detector effects and measurement uncertainties. In addition to the uncertainties of the individual observables, the uncertainty on the reconstructed energy of the primary particle will impact the total accuracy for the mass. Therefore, it is an advantage of the new radio-muon mass estimator that it only weakly depends on the energy of the primary particle. Compared to the electron-muon ratio, the normalization for the energy has a relatively small influence on the figures of merit for the radio-muon combination and for $X_{\mathrm{max}}$. In particular $X_{\mathrm{max}}$ and the energy content of the electromagnetic shower component can be measured not only by radio arrays but also by other techniques. Due to their similarities in the sensitivity to the electromagnetic shower component, qualitatively similarly results can be expected for the combination of muon detection with fluorescence or air-Cherenkov light. However, these techniques suffer from their limited duty cycle and atmospheric light absorption. The latter hampers the air-Cherenkov measurement particularly for inclined showers. Hence, only the combination of muon detectors with either fluorescence or radio detectors is expected to provide high mass sensitivity for large zenith angles. Coincident events with the fluorescence, muon, and radio detectors of the upgraded Pierre Auger Observatory will enable an independent cross-check of the new mass estimator.

The influences of realistic detector responses, Poissionian fluctuations due to limited detector sizes, measurement uncertainties, and background is investigated for the combination of the AMIGA Muon Detector and AERA of the Pierre Auger Observatory in Refs. [21, 30]. As expected, these effects slightly degrade the mass-separation power, but generally the high potential for mass-composition studies is confirmed. Dedicated simulation studies need to be done to estimate the full potential of the radio-muon combination for showers more inclined than 55$^{\circ }$, since the reference distance in this study was optimized for the range of zenith angles until 55$^{\circ }$ accessible by AMIGA.

As for other methods for the estimation of the mass of the primary particle based on air-shower observable, the use of hadronic interaction models comes with unavoidable systematic uncertainties. As studied by several experiments [31], all available hadronic interactions models have a deficit in the prediction of muons, which is largest at the highest energies [32]. We do not expect that this discrepancy in the muon number will have a significant influence on the mass separation power investigated in this work, which relies on relative differences in the muon content of proton and iron showers. However, the absolute scale of the muon density is shifted, that has to be taken into account when comparing simulations to measured data and when interpreting data based on simulations.

Finally, the novel technique for mass estimation can be applied to other experiments. While dedicated simulation studies will be needed for each experiment, the general findings of this study should be transferable, since neither high-energy muons nor the radio signal suffer from significant absorption in the atmosphere. Only for sites at higher altitudes the effect of partly clipping the air shower at the observation level needs to be investigated more carefully, in particular for vertical and mildly inclined showers.

A variety of activities is already ongoing. The Pierre Auger Observatory is currently being upgraded with scintillators on top of the water-Cherenkov detectors to disentangle the muonic and electromagnetic components of showers up to zenith angles of 60$^{\circ }$. Investigations are ongoing to equip each surface detector station in addition with a radio antenna, which will allow to measure the electromagnetic component as well for inclined shower [33, 34]. Possibilities to lower the energy threshold of the radio detection technique by the right choice of the frequency band were investigated in Ref. [35]. This can be applied to search for air showers induced by PeV photons using the radio-muon combination for gamma-hadron separation. Therefore, activities are ongoing to enhance the IceTop particle detector array at the South Pole with radio antennas [36]. Furthermore, the TAIGA facility comprises muon detectors and radio antennas (Tunka-Rex), at which the technique could directly be applied [37]. The Giant Radio Array for Neutrino Detection (GRAND) is a huge antenna array planned in China focusing on inclined showers [38]. GRANDproto300, its next stage prototype, will be additionally equipped with Auger-like particle detectors for the purpose of applying the radio-muon method to cosmic-ray air showers.

8 Conclusion

In this work, a novel technique is developed to estimate the mass of ultra-high energy cosmic rays by combining the muon signal and the radio emission of air showers. The muonic and electromagnetic components of air showers induced by proton and iron primaries were analyzed based on air-shower simulations. The size of the muonic component can be observed by the muon density $\rho _{\mu }(\text {r}_{\text {ref}})$ at a reference distance to the shower axis. 600 m was found to be a distance with a strong mass sensitivity for showers with zenith angles below 55$^{\circ }$, and the mass sensitivity may further be enhanced in particular for more inclined showers by varying the reference distance as a function of zenith angle. The radiation energy $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ , i.e. the energy contained in the radio emission, is correlated to the size of the electromagnetic component. The ratio $\rho _{\mu }^{600} / \sqrt{{\mathrm{S}}_{{\mathrm{RD}}}^{\rho _{\theta }}}$ represents a mass-sensitive parameter that is larger for iron than for proton showers. The mass sensitivity of the radio-muon combination was investigated and compared to established methods using solely particle measurements, or using the shower maximum $X_{\mathrm{max}}$. With the presented approach, the radio-muon combination features a mass separation power slightly larger than that of $X_{\mathrm{max}}$ for all zenith angles. Only if the energy of the primary particle is known accurately enough, the traditional observable of the electron-muon ratio provides better mass separation for near-vertical showers with zenith angles smaller than 35$^{\circ }$. Otherwise, and generally for more inclined showers, the new method of the radio-muon combination features superior mass estimation.

This emphasizes the potential for this new mass-sensitive parameter in particular for inclined showers. At large zenith angles, the radio emission spans over a large area on ground, which makes sparse detection array and thus large-scale applications feasible. The results show that the novel technique provides additional accuracy for mass measurements of cosmic rays on a per-event level. This is essential for further progress in answering open questions about the origin of ultra-high energy cosmic rays, since not only the total flux, but also the mass composition is expected to feature an anisotropy in the arrival directions [13]. Therefore, the scientific potential of existing particle-detector arrays can easily be enhanced by adding radio antennas, and future air-shower arrays can be planned to feature both, muon and radio detectors.

Appendix A: Calculation of the radiation energy from the electromagnetic shower energy

The radiation energy is correlated to the electromagnetic shower energy as shown in Ref. [25]. The correlation was modeled in [26] based on CoREAS simulations [39] for showers with zenith angles up to 80$^{\circ }$, including various corrections for dependencies on the arrival direction. The model is used in this work to calculate the radiation energy from the shower energy extracted from the CORSIKA simulations.

The radiation energy slightly depends on the arrival direction and the angle $\alpha$ to the geomagnetic field. The geomagnetic fraction of the radiation energy is influenced by the magnitude of the geomagnetic field $B_{\text {Earth}}$ as well as the angle $\alpha$ between the shower axis and $B_{\text {Earth}}$. The radiation energy scales with $\sin ^2{\alpha }$ because of the coherent nature of the radio emission. The charge excess fraction a of the radiation energy grows with the atmospheric density $\rho _{X_{\mathrm{max}}}$ at $X_{\mathrm{max}}$. The atmospheric density decreases with altitude. Thus, $\rho _{X_{\mathrm{max}}}$ depends on the zenith angle and the altitude of the shower maximum $X_{\mathrm{max}}$, which is generally higher in the atmosphere for showers induced by heavier particles. Furthermore, there is a second order dependence on $\rho _{X_{\mathrm{max}}}$. Whereas the radio emission depends on the geometric distance (in m), the air shower develops according to the atmospheric depth (in g/cm$^2$). The ratio between the geometric distance and the atmospheric depth is higher for regions of lower atmospheric density. This leads to a slightly larger radio emission, if $X_{\mathrm{max}}$ is higher in the atmosphere. Therefore, the radiation energy $E_{\text {rad}}$ is corrected for these dependencies by

$$\begin{aligned} S_{\mathrm{RD}}^{\rho _{X_{\mathrm{max}}}}= & {} \frac{E_{\text {rad}}}{a(\rho _{X_{\mathrm{max}}})^2 + \left( 1 - a(\rho _{X_{\mathrm{max}}})^2 \right) \cdot \sin ^2{\alpha }} \nonumber \\&\cdot \frac{1}{\left( 1 - p_0 + p_0 \cdot \exp { \left[ p_1 \cdot \left( \rho _{X_{\mathrm{max}}} - \langle \rho \rangle \right) \right] } \right) ^2 },\nonumber \\ \end{aligned}$$

(A.1)

where $p_0 = 0.251 \pm 0.006$ and $p_1 = -2.95 \pm 0.06$ m$^3$/kg, and $\langle \rho \rangle = 0.65$ kg/m$^3$ is the atmospheric density at the average $\langle X _\mathrm {max} \rangle = 669\hbox { g/cm}^2$ for an average zenith angle of 45$^{\circ }$ and a primary energy of $10^{18}$ eV for a 50%-proton / 50%-iron composition.

The correlation of the true radiation energy S_RD, after normalization according to the arrival direction, with the electromagnetic shower energy $E_{\text {em}}$ is modeled in [26] by

$$\begin{aligned} S_{\text {RD}}^{\rho _{X_{\mathrm{max}}}} = A \cdot 10^7~\text {eV} \cdot \left( \frac{E_{\text {em}}}{10^{18}~\text {eV}} \right) ^B \end{aligned}$$

(A.2)

with $A = 1.683 \pm 0.004$ and $B = 2.006 \pm 0.001$. The radiation energy used for this model was later found to be underestimated by about 11% in the simulations used for this parametrization due to settings in the CoREAS simulations [40]. Since this only affects the absolute scale, but not the relative difference between proton and iron showers, this is not considered here.

$X_{\mathrm{max}}$ , which is used in the correlation, is often not accessible in an experiment and in case it is measured, there are additional measurement uncertainties. Therefore, a correction dependent on the zenith angle, for which the measurement uncertainties are much smaller, is formulated based on the mean $X_{\mathrm{max}}$ at the respective zenith angles:

$$\begin{aligned} S_{\text {RD}}^{\rho _{\theta } \prime }= & {} \frac{E_{\text {rad}}}{a(\rho _{\theta })^2 + \left( 1 - a(\rho _{\theta })^2 \right) \cdot \sin ^2{\alpha }} \nonumber \\&\cdot \frac{1}{\left( 1 - p_0 + p_0 \cdot \exp { \left[ p_1 \cdot \left( \rho _{\theta } - \langle \rho \rangle \right) \right] } \right) ^2 }. \end{aligned}$$

(A.3)

In addition, another zenith angle dependent effect has to be taken into account. Depending on the observation level of a detector, the shower might not be fully developed at the altitude of detection. Hence, a part of the shower is clipped before the radio emission of this part is released. The magnitude of this clipping effect depends on the distance between the observer and $X_{\mathrm{max}}$. The radiation energy investigated in the following is corrected for this effect by

$$\begin{aligned}&S_{\text {RD}}^{\rho _{\theta }} \nonumber \\&\quad = \frac{S_{\text {RD}}^{\rho _{\theta } \prime }}{1 - \exp { \left( -8.7~\text {cm}^2 / \text {kg} \left( D_{X_{\mathrm{max}}} + 0.29~\text {kg} / \text {cm}^{2} \right) ^{1.89} \right) }} \nonumber \\ \end{aligned}$$

(A.4)

where $D_{X_{\mathrm{max}}}$ is the distance between the observer and $X_{\mathrm{max}}$ in kg/cm$^2$. The effect of clipping is small for the present study. For the simulations used in this work, the size of the correction is at most 10% of the total radiation energy, and for the majority of the simulations it is smaller than 2%. Clipping will be more relevant for higher energies or for observatories at higher altitude.

In summary, in this work $S_{\text {RD}}^{\rho _{X_{\mathrm{max}}}}$ is calculated from the electromagnetic shower energy of the full shower by Eq. A.2. Then, $E_{\text {rad}}$ is calculated using Eq. A.1, and the radiation energy $S_{\text {RD}}^{\rho _{\theta } \prime }$ is corrected for the atmospheric density depending on the zenith angle by Eq. A.3. Finally, the radiation energy is clipped according to the observation level by Eq. A.4 to gain $S_{{\mathrm{RD}}}^{\rho _{\theta }}$ used as observable in the present work.