A new air-shower observable to constrain hadronic interaction models

TL;DR

This paper introduces a novel observable based on muon energy spectra in air showers to better constrain hadronic interaction models, especially in the critical energy range of 10^{17.5} to 10^{18} eV.

Contribution

A new observable ( mu) derived from muon ground energy spectra is proposed to effectively constrain hadronic interaction models in air shower simulations.

Findings

mu can distinguish between different hadronic models.

The observable is effective within current experimental resolutions.

The study focuses on the energy range 10^{17.5} to 10^{18} eV.

Abstract

The energy spectrum of muons at ground level in air showers are studied and a new observable is proposed to constrain hadronic interaction models used in air shower simulations. An asymmetric Gaussian function is proposed to describe the muon ground energy spectrum and its parameters are studied regarding primary particle, energy and hadronic interaction models. Based on two realistic measurements of the muon density at a given distance from the shower axis, a new observable (\rmu) is defined. Considering realistic values of detector resolutions and number of measured events, it is also shown \rmu can be successfully used to constrain low and high hadronic interaction models. The study is focused in the energy range between $10^{17.5}$ and $10^{18.0}$ eV because of the importance of this interval for particle physics and astrophysical models. The constraining power of the new observable is shown to be large within current experimental capabilities.

Equations

$$DX=\frac{X_{\mathrm{gr},\mathrm{vert}}}{\cos \theta }-X_{\max },$$

$$\frac{\text{d}{N}_{\mu }}{\text{d}x}=\{\begin{array}{ll}{N}_0\exp \left[-\frac{1}{2}{\left(\frac{x-\eta }{\sigma }\right)}^2\right],& \text{if x<η}\\ & \\ N_0\exp \left[-\frac{1}{2}{\left(\frac{x-\eta }{\alpha \sigma }\right)}^2\right],& \text{if x>η}\end{array}$$

$$r_{\mu }={\sec }^{\beta }\theta \frac{S_{\mu }^{\mathrm{bur}}}{S_{\mu }^{\mathrm{sur}}},$$

$$\text{Merit Factor}=\frac{\mid ⟨r_{\mu }{⟩}_a-⟨r_{\mu }{⟩}_b\mid }{\sqrt{{\sigma }_a^2+{\sigma }_b^2}},$$

Full text

††journal: Astroparticle Physics

1 Introduction

The interaction of ultra high energy cosmic rays with atmospheric nuclei allows us to access hadronic interactions at energies much beyond the reach of man-made accelerators. However several properties of the cosmic ray phenomena and of the detection techniques impose limitations on the available information to study the highest energetic interactions. Experiments are only able to measure the air shower produced by the interaction of the cosmic particles with the nucleus of atoms in the atmosphere. The study of elementary properties of particle physics is done by relating global shower parameters to the properties of the hadronic interactions.

The depth at which the shower reaches its maximum number of particles ($X_{\rm{max}}$) and the muon component are the most used shower features related to properties of particle interactions. At the same time, they are also strongly sensitive to the type of the shower primary particle [1]. Given that the primary particle type cannot be determined on event-by-event basis, different primary compositions must be considered in the study of the interaction properties. Usually, the large possibility of primary particles, from proton to iron nuclei, makes it impossible to disentangle the mass of the cosmic particle from particle interaction properties.

This paper proposes an analysis procedure to constrain hadronic interaction models by the measurements of the muon density by two different detectors. The idea proposed here is to compare the predictions of the models to measurements based on the muon density at ground level. It will be shown in Sec. 4 that the muon ground energy spectrum has valuable information, which allow to discriminate among different hadronic interaction models.

Hadronic interaction models used in Monte Carlo simulations of air showers are limited to phenomenological approaches tested and tuned to collider data  [2]. Three hadronic models for high energy interaction ($>80$ GeV) and three hadronic models for low energy interaction ($<80$ GeV) were investigated. A new and realistic air-shower observable based on two measurements of the muon density by different detectors is proposed in this paper and it is demonstrated to be powerful enough to constrain the hadronic interactions models.

The analysis procedure proposed here does not aim to infer properties of the particle interactions, instead, the new proposed parameter has power to test the predictions of the hadronic interaction models within realistic experimental conditions. By comparing the prediction of the models to future measurements of the new parameter, it will be possible to select the best model and guide the way towards a better understanding of the underlying assumptions taken by that model.

The work is based on Monte Carlo simulations of air showers. The new proposed parameter is studied in the energy range from $10^{17.5}$ to $10^{18.0}$ eV. The energy range is limited in order to minimize the effect of systematic uncertainties in energy reconstruction and also because of the importance of this range, which is the transition energy range from collider data to the extrapolation domain. It is shown that the new parameter can be used to test hadronic interaction models without knowledge on the primary composition.

The paper is organized as follows: Sec. 2 describes the simulations, Sec. 3 shows a study on the energy spectra of muons at the ground level, Sec. 4 proposes the new parameter and studies its behavior with relation to primary mass, energy and detector properties, Sec. 5 quantifies the constraining power of the new parameter and Sec. 6 concludes the work.

2 Simulations

Extensive air showers were simulated using the Monte Carlo code CORSIKA v7.500 [3]. The applied thinning factor was $10^{-2}$ for the electromagnetic component and $10^{-4}$ for the hadronic component, with the maximum weight for any particle set to $100$. It was verified that the choice of this thinning configuration does not result in a bias in the present analysis. The energy of the primary particles was sampled continuously between $10^{17.5}$ and $10^{18.0}$ eV, following a power law energy spectrum with index -3. The arrival directions were sampled following a uniform distribution in solid angle, up to a maximum zenith angle of $60^{\circ}$. Three primaries were simulated: proton, nitrogen nucleus and iron nucleus. The minimum energy of particles simulated in air showers were set to $0.3$ GeV for hadrons and muons, and $0.003$ GeV for electrons and photons. Three hadronic interaction models were used for high energies interactions ($>80$ GeV): QGSJetII-04 [4], EPOS-LHC [5, 6] and Sibyll2.3 [7], and three for low energies ($<80$ GeV): FLUKA [8], GHEISHA [9] and UrQMD [10]. For each selected combination of hadronic interaction model and primaries particle, 1200 showers were generated. The ground altitude (1400 m above sea level) was chosen to be the one at the Pierre Auger Observatory.

The energy spectrum of muons at ground level is the shower feature to be considered in this paper. The energy spectra were built by collecting from the simulated air showers all muons reaching ground in a lateral distance between $425$ and $475$ meters from the shower axis. Although a full detector reconstruction is not applied, in Sec. 4.1 the detector effects are taken into account by artificial smearing the shower observables around its simulated values. Detector thresholds and geometry are considered as well.

To sample showers at the same stage of development, the commonly used $DX$ parameter is used. $DX$ is the slant atmospheric depth between the shower maximum and the ground. It is given by

[TABLE]

where $X_{\rm{gr,vert}}$ is the vertical slant depth of the ground111$X_{\rm{gr,vert}}=870$ g/cm2 for the Pierre Auger Observatory, $\theta$ is the zenith angle of the shower axis and $X_{\rm{max}}$ is the depth in which the shower reaches its maximum. For each simulated event, the $\theta$ corresponds to the true zenith angle as set in the input of the simulation and $X_{\rm{max}}$ was taken directly from the Gaisser-Hillas function fitted to the longitudinal energy deposit profile.

3 Characterization of the muon ground energy spectrum

In this section, the simulations described in Sec. 2 are used to characterize the energy spectrum of muons at ground level and to study its relations with primary mass and hadronic interaction models.

Fig. 1 shows examples of the ground energy spectrum of muons for six simulated events, which differ by the primary particle and the shower geometry. The low and high energy hadronic interaction models are FLUKA and QGSJetII-04, respectively. Examples have been selected to illustrate the general shape of the muon spectrum for different primary particles and extreme values of $DX$. The left-hand column shows deep showers, with relatively small $DX$ values, while the right-hand column shows shallow showers, with relatively high values of $DX$. The normalization and the mean of the distributions are clearly different but the overall shape is very similar.

As shown by the red lines in Fig. 1, the ground energy spectrum of muons is well described by the following asymmetric Gaussian function

[TABLE]

where $x=\log_{10}(E/\rm{GeV})$. $N_{0}$ is the normalization parameter and it is correlated with the total number of muons in the shower. $\eta$ is the mode of the energy distribution, and it is strongly correlated with the average energy of muons reaching ground between $425$ and $475$ meters distance from the shower axis. $\sigma$ and $\alpha$ give the width of the distribution, $\alpha$ being the parameter that measures the degree of asymmetry of the distribution.

Muon energy spectra of all simulated showers were fit by the function presented in Eq. 2, in which $N_{0}$, $\eta$, $\sigma$ and $\alpha$ were taken as free parameters. The fitting was performed using a binned maximum likelihood method with Poissonian probability distribution functions.

Since the ground energy spectrum of muons is well described by the asymmetric Gaussian function, one may study its dependencies on the energy, primary particle and hadronic interaction models through the evolution of the parameters $N_{0}$, $\eta$, $\sigma$ and $\alpha$ with $DX$. Fig. 2 shows the $DX$ evolution of the four parameters for three different energy intervals, Fig. 3 for three primary particles and Fig. 4 for five combinations of high and low energy hadronic interaction models.

It is clear from Fig. 2 that the normalization ($N_{0}$) is the only property of ground muon energy spectrum which shows a significant dependence on the primary energy. As expected, $N_{0}$ also depends strongly on the primary particle and hadronic interaction model (see Figs. 3 and 4), which reflects the very known behavior of number of muons in air showers. Concerning the peak position of the distributions, $\eta$, a strong evolution with $DX$ is observed, revealing the shift of the average energy of muons to higher values with increasing $DX$. Besides normalization and peak positions, changes on the overall shape of the energy spectrum of muons can be evaluated through the parameters $\sigma$ and $\alpha$. These parameters clearly show a very weak evolution with $DX$, and a nearly null dependence on the primary energy and primary particle and a relatively very small dependence on the hadronic interaction models.

The analysis of Fig. 4 also repeats known though less popular lessons: the effect of low energy hadronic interaction models are as important as the high energy one regarding the description of the muonic component in air-shower simulations [11]. The differences of $N_{0}$ between QGSJetII-04/FLUKA and QGSJetII-04/UrQMD or GHEISHA are of the same order or larger than the largest difference between the high energy interaction models. The average energy of muons, correlated with $\eta$, is larger for GHEISHA, followed by UrQMD and then FLUKA. The differences in $\eta$ due to the low energy hadronic interaction models are as large as the differences due to the high energy hadronic interaction models. Regarding $\sigma$ and $\alpha$ parameters, one can see in Fig. 4 they show the same weak $DX$ evolution for all the hadronic model combinations. Furthermore, the effect of low energy hadronic models is again clear as shown by the difference between FLUKA and GHEISHA, or UrQMD.

Finally, the analysis of Fig. 4 also opens new possibilities. It turns out that $\eta$ is strongly sensitive to the hadronic interaction models, and at the same time, it does not show any significant dependence on the primary energy and on the primary particle. The lack of primary energy dependence is important experimentally to eliminate effects from the experimental energy scale, while the lack of primary particle dependence is an essential property to disentangle the hadronic interaction effects from the primary composition determination. The evolution of $\eta$ with $DX$ is simple (linear) and strong, which makes easy to study showers in different evolution stages.

All of these are indications that accessing experimentally the information carried by $\eta$ could be successfully used to constrain the hadronic interaction models. In the next section a new observable which is strongly correlated to $\eta$ is studied and its constraining power in realistic experimental conditions is tested.

4 An observable to test hadronic interaction models

Accessing $\eta$ directly is not possible for any running or planned UHECR experiment. Therefore, instead of proposing an unrealistic parameter, this study starts from a realistic experimental scenario and aims to find an observable which correlates to $\eta$. A generic experimental set-up is considered with two muon detector arrays with different amount of shielding leading to two energy thresholds. With such an experimental set-up, it is possible to measure the integral of the energy spectrum or the energy density of muons above two energy thresholds.

Using the simulations described in Sec. 2 the density of muons in the lateral distance from 425 to 475 m from the shower axis was calculated for a surface ($S_{\mu}^{\rm{sur}}$) and a buried ($S_{\mu}^{\rm{bur}}$) detectors and a new parameter is defined:

[TABLE]

where $\theta$ is the zenith angle of the primary particle. The term $\sec^{\beta}\theta$ compensates the zenith dependence on the energy threshold and the effective area of the buried detectors. The surface detector ($S_{\mu}^{\rm{sur}}$) is considered to have $10$ m2 and the buried detector ($S_{\mu}^{\rm{bur}}$) to have $30/\cos\theta$ m2, which are motivated by water-Cherenkov stations and flat buried scintillators respectively. $\beta$ was determined to minimize the primary mass dependence of $r_{\mu}$. In Fig. 5 the $r_{\mu}$ dependence on $\beta$ is shown for one hadronic interaction models, which justify the choice of $\beta=0.6$ as the value that minimize the difference between primaries. The detector features considered is described in the next section.

Figs. 6 and 7 show the distributions of $r_{\mu}$ for several energy thresholds of the buried detector. Three cases are shown in which the vertical muon energy threshold of the buried detector ($E_{\rm{vert},\mu}^{\rm{th}}$) is changed. The effective energy threshold for a muon with incident zenith angle $\theta_{\mu}$ is $E_{\rm{vert},\mu}^{\rm{th}}/\cos\theta_{\mu}$. For the surface detectors, the energy threshold is kept fixed at $0.3$ GeV. Fig. 6 compares the distributions for different high energy hadronic interaction models and Fig. 7 for different low energy hadronic interaction models. The different degrees of separation of the hadronic interaction models with different energy threshold is clear, pointing to the possibility to constrain the hadronic interaction models using $r_{\mu}$. In the next sections, the constraining power of $r_{\mu}$ is estimated and its correlation with $\eta$ is shown.

4.1 $r_{\mu}$ determination including detector characteristics

In this section, the effect of detector geometry, muon energy thresholds, detector resolutions and systematic uncertainties on $r_{\mu}$ are studied. To include the detector features in the proposed analysis, the general characteristics of the muon detectors of the Pierre Auger Observatory are considered. A surface detector with $10$ m2 of effective area is considered to measure muons with energies above $0.3$ GeV. Its general characteristics are inspired in the AugerPrime [12, 13, 14] design of water-Cherenkov stations with plastic scintillators on top. A second detector is considered based on the general characteristics of the Auger AMIGA detector [15, 16]. Those are $30$ m2 flat scintillator detectors buried $2.5$ m below the ground. The muon energy threshold for the buried detectors depends on the incident zenith angle of the particle and it is given by $E_{\rm{th}}=\beta\rho h/\cos\theta_{\mu}$, where $\beta=1.808$ MeV cm2 g*-1* is the fractional energy loss per depth of standard rock, $\rho=1.8$ g cm*-3* is the soil density, $h=2.5$ m is the vertical depth of the detectors and $\theta_{\mu}$ is the zenith incidence angle of the muon. Because of its flatness, the effective collection area of the detector decrease by a factor $\cos\theta$, where $\theta$ is the zenith angle of the shower. The reconstruction of the muon density by an AMIGA-like detector has been shown to be satisfactorily possible for events with $\theta<45^{\circ}$ [17, 18, 19].

From now on, the results were obtained by considering the detector features shown above to calculate the muon density from the simulated air showers explained in Sec. 2. In this way, the most important properties of the detectors are considered without the need to perform a full detector simulation.

Fig. 8 shows the relation between the average value of $r_{\mu}$ ($\langle r_{\mu}\rangle$) and the average value of $\eta$ ($\langle\eta\rangle$) for showers with $325<DX/(\rm{g/cm}^{2})<375$ according to the detector properties explained above. Five combinations of low and high energy hadronic interaction models are shown in different colors and three primaries are shown in different marker styles. From Fig. 8, one can see a nearly linear relation between $\langle r_{\mu}\rangle$ and $\langle\eta\rangle$. Furthermore, the relatively large separation between the hadronic interaction model combinations by $\langle r_{\mu}\rangle$ is clear, which shows $r_{\mu}$ is a good observable to constrain hadronic interaction properties.

The resolutions on the muon density reconstruction were taken into account by applying a Gaussian smearing on the true signal obtained from the simulations. $S_{\mu}^{\rm{sur}}$ resolution was considered to vary in the range from $10\%$ to $20\%$ [12, 14] and $S_{\mu}^{\rm{bur}}$ resolution in the range from $5\%$ to $10\%$ [18, 19]. Fig. 9 shows the effect of the resolution in the calculation of $\langle r_{\mu}\rangle$. The Gaussian smearing were performed 2000 times and the $\langle r_{\mu}\rangle$ is shown as a function of $DX$. The standard deviation of $r_{\mu}$ distributions is shown by the error bars. Three cases are shown in which the resolution on both $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ vary. The systematic effect of the detector resolution in the determination of $\langle r_{\mu}\rangle$ is smaller than 5%.

Besides the experimental resolutions, systematic uncertainties on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ can be originated from the detection and reconstruction procedures. Typically the most significant systematic uncertainty on muon density measurements are due to systematic uncertainties on the shower energy determination, which affects both $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ in the same magnitude. To evaluate the effect of systematic uncertainties on $r_{\mu}$, the simulated $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ were shifted artificially and the resulting $r_{\mu}$ were calculated. First it was considered the shifts on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ are totally correlated, which means the same magnitude and direction. This case represents the energy reconstruction uncertainty effect. Fig. 10 shows the $r_{\mu}$ as a function of $DX$, for one hadronic interaction models combination, for different cases in which $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ were shifted by a factor $1+\delta^{top}$ and $1+\delta^{bur}$ respectively. In Fig. 10 left panel it is shown the effect of a 10% shift on both muon density at the same direction. Clearly, correlated systematic uncertainties on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ have an insignificant effect on $r_{\mu}$. In Fig. 10 right panel it is shown the effect of systematic shifts of 2.5% on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ in opposite directions. The magnitude of the $\langle r_{\mu}\rangle$ deviation is of order 0.05. The consequences of this deviation on the separation between different hadronic interaction models is discussed in Sec. 5.

4.2 $r_{\mu}$ dependence on primary mass and energy

For the following analysis, the $DX$ range was defined to preserve a good statistics in all $DX$ bins and it goes from $100$ to $600$ g/cm2 divided in $5$ bins of $100$ g/cm2. The upper bound of 600 g/cm2 is highly influenced by the shower zenith angle limitation at $\theta<45^{\circ}$ due to the buried detector features.

Fig. 11 shows the evolution of $\langle r_{\mu}\rangle$ as a function of $DX$ for different energy ranges. All simulated primary particles are include (p, N and Fe). Hadronic interaction model combination here is QGSJetII-04/FLUKA. A better visualization of the differences in $\langle r_{\mu}\rangle$ is seen in the bottom panel of Fig. 11, where $\Delta r_{\mu}$ are the differences with relation to the average value of $\langle r_{\mu}\rangle$ for the three energy ranges considered. The differences of $\langle r_{\mu}\rangle$ in all energy intervals is smaller than 1% for the entire $DX$ range. The energy independence of $\langle r_{\mu}\rangle$ is expected, since $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ evolve similarly with energy in the range from $10^{17.5}$ to $10^{18.0}$ eV. The lack of energy dependence is an advantage because it allows the analysis of events in a large energy interval, increasing significantly the available statistics. Furthermore, it also contribute to diminishing any effect due to the experimental energy scale.

The primary mass dependence of the $\langle r_{\mu}\rangle$ is shown in Fig. 12. The hadronic interaction model combination shown is again QGSJetII-04/FLUKA. Very similar results were obtained for all combinations of models. The bottom panel shows the $\Delta\langle r_{\mu}\rangle$. The observed primary mass dependence is below 2%. The dependence of $\langle r_{\mu}\rangle$ with the primary particle was minimized by choosing $\beta=0.6$. The lack of $\langle r_{\mu}\rangle$ dependence on the primary particle is a great advantage because it disantangles the study of the hadronic interaction properties from the determination of the primary particle type.

5 Results: constraining hadronic interaction models with $r_{\mu}$ parameter

In this section, the capacity to constrain hadronic interaction models by measuring $r_{\mu}$ is demonstrated by studying the $DX$ evolution of $\langle r_{\mu}\rangle$ for different combinations of low and high energy hadronic models. Detector resolution is taken into account as explained above. The effect of a limited number of events is considered here. The total number of simulated air showers used (3500) is approximately the number of hybrid events to be measured by the Pierre Auger Observatory infill array in $2$ years of operation. The infill array consists of $750$ m spaced water-Cherenkov stations spread over an area of $23.5$ km2. In this same area the AMIGA-Grande and AugerPrime muon detectors are going to be installed. Considering the duty cycle of the fluorescence detectors to be $14\%$, the total exposure of this experimental set-up is $4.32$ km2.sr.yr. Taking into account the cosmic-rays flux between $10^{17.5}$ and $10^{18.0}$ eV, the expected number of events to be measured per year is 1806.

Figs. 13 and 14 show the $\langle r_{\mu}\rangle$ as a function of $DX$ for different combinations of hadronic interaction models. In Fig. 13 the three high energy hadronic models are shown in combination with one low energy hadronic interaction model: FLUKA. In Fig. 14 the three low energy hadronic models are shown in combination with one high energy hadronic interaction model: QGSJetII-04. The bottom panels show the $\Delta\langle r_{\mu}\rangle$. The worst case for the detector resolution was considered: $10$% for $S_{\mu}^{\rm{bur}}$ and $20$% for $S_{\mu}^{\rm{sur}}$. Figs. 13 and 14 show that even with a relatively poor detector resolution a clear separation between hadronic interaction models is achieved. In Fig. 13 it is observed that EPOS-LHC can be distinguished from QGSJetII-04 and Sibyll2.3, while in Fig. 14 it is seen that GHEISHA can be distinguished from UrQMD and FLUKA.

To better quantify the discriminating power of $\langle r_{\mu}\rangle$, the commonly used Merit Factor can be used. It is defined as:

[TABLE]

where $a$ and $b$ refer to any two hadronic interaction model combination and the $\sigma$’s are the standard deviations of $\langle r_{\mu}\rangle$.

Fig. 15 shows the Merit Factor as a function of $DX$ for the same experimental resolutions and statistics described above. The left-hand panel refers to hadronic model combinations with different high energy models, and on the right-hand panel with different low energy models. The best Merit Factor is $3.0$ between EPOS-LHC and QGSJetII-04 and $2.3$ between FLUKA and GHEISHA.

Figs. 16 and 17 show the Merit Factor as a function of the number of events and detector resolution for one particular $DX$ bin: $300<DX/(\rm{g/cm}^{2})<400$. The resolutions on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ were re-scaled by a common factor $f$, being that $\sigma_{\rm{bur}}=0.1fS_{\mu}^{\rm{bur}}$ and $\sigma_{\rm{top}}=0.2fS_{\mu}^{\rm{top}}$. The effect of the number of events was calculated by re-scaling the standard deviation of the muon densities by a factor $\sqrt{N_{\rm{sim}}/N}$, where $N_{\rm{sim}}$ is the number of simulated showers and $N$ is the number of showers in each case.

Fig. 16 shows that it is possible to reach large Merit Factor values ($>5$) for the separation between EPOS-LHC and QGSJetII-04/Sibyll2.3 using a reasonably small number of events ($<6000$) considering realistic detector resolutions ($\sigma_{\rm{bur}}/S_{\mu}^{\rm{bur}}<0.06$ and $\sigma_{\rm{top}}/S_{\mu}^{\rm{top}}<0.13$). On the other hand, the separation between Sibyll2.3 and QGSJetII-04 is small for any resolutions and number of events, which is expected because of their similar values of $\eta$. The same conclusions can be drawn about the low energy hadronic interactions models from Fig. 17. $\langle r_{\mu}\rangle$ provides a very good separation between FLUKA and GHEISHA/UrQMD, but the separation power is limited for GHEISHA and UrQMD.

6 Conclusions

This paper studies the ground muon energy spectrum of air showers and proposes an analysis procedure to constrain hadronic interaction models used in Monte Carlo simulations. In Sec. 3, it was shown that the energy distribution of muons at ground level can be well described by an asymmetric Gaussian function with four parameters. The study of the evolution of the four parameters with $DX$ concluded that the overall shape of the energy spectrum of muons does not depend strongly on the combination of low and high energy hadronic interactions models. Also, it was verified that the average muon energy, or the parameter $\eta$, is sensitive to the model combination and presents a strong evolution with $DX$.

However, $\eta$ is not an easy parameter to be measured. Therefore in Sec. 4 a new and experimentally motivated parameter ($r_{\mu}$) is proposed and its correlation with $\eta$ is shown. $r_{\mu}$ dependencies on the primary mass and energy were proven to be insignificant which allows on to disentangle the composition and the hadronic interaction studies as well as minimize the effect of systematic uncertainties in the energy reconstruction.

The general properties of the current and proposed muon detectors of the Pierre Auger Observatory are considered to study $r_{\mu}$ under realistic experimental limitations. Sec. 5 shows that the discrimination power of $r_{\mu}$ is significantly large. EPOS-LHC can be separated from Sibyll2.3 and QGSJetII-04 with large Merit Factor ($>5$) using a reasonably small number of events ($<6000$). Sibyll2.3 and QGSJetII-04 show similar average muon energy at ground and therefore can, in the best case, be discriminated with Merit Factor $\sim 2$ with about $9000$ events. As for the low energy hadronic models, FLUKA can be separated from GHEISHA and UrQMD with large Merit Factor ($>5$) using a reasonably small number of events ($<6000$). GHEISHA and UrQMD can, in the best case, be discriminated with Merit Factor $\sim 0.8$ with about $9000$ events. It was shown that correlated systematic uncertainties on $S_{\mu}^{\rm{sur}}$ and $S_{\mu}^{\rm{bur}}$ have insignificant influence on $\langle r_{\mu}\rangle$. Uncorrelated oppositive systematics of order of 2.5% can generate a systematic uncertainty on $\langle r_{\mu}\rangle$ of the same order of the separation betweem the hadronic interaction models. This means that for a realistic application of $r_{\mu}$ analysis, the systematic uncertainties on the muon densities measured by the two detectors have to be correlated (more realistic case). If the systematic uncertainties of the two detectors are in oppositive directions (less realistic case) they have to be smaller than 2.5% in order to keep the discrimination power of the hadronic interaction models.

It was also shown that $r_{\mu}$ is a very robust parameter which can be used irrespective of ignorance of the primary mass composition to test the hadronic interaction models. Constrains imposed by an analysis based on $r_{\mu}$ can in a short period of time contribute to the solution of the know problems with muon production in extensive air shower [20, 21].

Acknowledgements

The authors would like to thank Roger Clay for the review of the manuscript on behalf of the Pierre Auger Collaboration. RRP thanks the financial support given by FAPESP (2014/10460-1). VdS thanks the support of the Brazilian population via CNPq and FAPESP (2015/15897-1, 2014/19946-4).