# A new air-shower observable to constrain hadronic interaction models

###### Abstract

The energy spectrum of muons at ground level in air showers is 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 () is defined. Considering realistic values of detector resolutions and number of measured events, it is also shown can be successfully used to constrain low and high energy hadronic interaction models. The study is focused in the energy range between and 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.

###### keywords:

extensive air showers, muons, hadronic interaction models^{†}

^{†}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 () 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 Kampert2012 (). 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 Engel2011 (). Three hadronic models for high energy interaction ( GeV) and three hadronic models for low energy interaction ( 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 to 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 Heck1998a (). The applied thinning factor was for the electromagnetic component and for the hadronic component, with the maximum weight for any particle set to . 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 and 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 . Three primaries were simulated: proton, nitrogen nucleus and iron nucleus. The minimum energy of particles simulated in air showers were set to GeV for hadrons and muons, and GeV for electrons and photons. Three hadronic interaction models were used for high energies interactions ( GeV): QGSJetII-04 Ostapchenko2010 (), EPOS-LHC Werner2007 (); Pierog2013 () and Sibyll2.3 Riehn2015 (), and three for low energies ( GeV): FLUKA Fluka (), GHEISHA GHEISHA1985 () and UrQMD Bleicher1999 (). 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 and 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 parameter is used. is the slant atmospheric depth between the shower maximum and the ground. It is given by

(1) |

where is the vertical slant depth of the ground^{1}^{1}1 g/cm for the Pierre Auger Observatory, is the zenith angle of the shower axis and is the depth in which the shower reaches its maximum. For each simulated event, the corresponds to the true zenith angle as set in the input of the simulation and 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 . The left-hand column shows deep showers, with relatively small values, while the right-hand column shows shallow showers, with relatively high values of . 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

(2) |

where . is the normalization parameter and it is correlated with the total number of muons in the shower. is the mode of the energy distribution, and it is strongly correlated with the average energy of muons reaching ground between and meters distance from the shower axis. and give the width of the distribution, 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 , , and 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 , , and with . Fig. 2 shows the 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 () is the only property of ground muon energy spectrum which shows a significant dependence on the primary energy. As expected, also depends strongly on the primary particle and hadronic interaction model (see Figs. 4 and 3), which reflects the very known behavior of number of muons in air showers. Concerning the peak position of the distributions, , a strong evolution with is observed, revealing the shift of the average energy of muons to higher values with increasing . Besides normalization and peak positions, changes on the overall shape of the energy spectrum of muons can be evaluated through the parameters and . These parameters clearly show a very weak evolution with , 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 Meurer2006 (). The differences of 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 , is larger for GHEISHA, followed by UrQMD and then FLUKA. The differences in due to the low energy hadronic interaction models are as large as the differences due to the high energy hadronic interaction models. Regarding and parameters, one can see in Fig. 4 they show the same weak 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 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 with 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 could be successfully used to constrain the hadronic interaction models. In the next section a new observable which is strongly correlated to is studied and its constraining power in realistic experimental conditions is tested.

## 4 An observable to test hadronic interaction models

Accessing 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 . 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 () and a buried () detectors and a new parameter is defined:

(3) |

where is the zenith angle of the primary particle. The term compensates the zenith dependence on the energy threshold and the effective area of the buried detectors. The surface detector () is considered to have m and the buried detector () to have m, which are motivated by water-Cherenkov stations and flat buried scintillators respectively. was determined to minimize the primary mass dependence of . In Fig. 5 the dependence on is shown for one hadronic interaction models, which justify the choice of as the value that minimize the difference between primaries. The detector features considered is described in the next section.

Figs. 7 and 6 show the distributions of for several energy thresholds of the buried detector. Three cases are shown in which the vertical muon energy threshold of the buried detector () is changed. The effective energy threshold for a muon with incident zenith angle is . For the surface detectors, the energy threshold is kept fixed at 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 . In the next sections, the constraining power of is estimated and its correlation with is shown.

### 4.1 determination including detector characteristics

In this section, the effect of detector geometry, muon energy thresholds, detector resolutions and systematic uncertainties on 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 m of effective area is considered to measure muons with energies above GeV. Its general characteristics are inspired in the AugerPrime Engel2015 (); AugerPDR (); Gonzalez2016 () 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 AMIGA_ICRC (); AMIGA (). Those are m flat scintillator detectors buried 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 , where MeV cm g is the fractional energy loss per depth of standard rock, g cm is the soil density, m is the vertical depth of the detectors and is the zenith incidence angle of the muon. Because of its flatness, the effective collection area of the detector decrease by a factor , where 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 Supanitsky2008 (); Ravignani2015 (); Ravignani2016 ().

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 () and the average value of () for showers with 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 and . Furthermore, the relatively large separation between the hadronic interaction model combinations by is clear, which shows 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. resolution was considered to vary in the range from to Engel2015 (); Gonzalez2016 () and resolution in the range from to Ravignani2015 (); Ravignani2016 (). Fig. 9 shows the effect of the resolution in the calculation of . The Gaussian smearing were performed 2000 times and the is shown as a function of . The standard deviation of distributions is shown by the error bars. Three cases are shown in which the resolution on both and vary. The systematic effect of the detector resolution in the determination of is smaller than 5%.

Besides the experimental resolutions, systematic uncertainties on and 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 and in the same magnitude. To evaluate the effect of systematic uncertainties on , the simulated and were shifted artificially and the resulting were calculated. First it was considered the shifts on and are totally correlated, which means the same magnitude and direction. This case represents the energy reconstruction uncertainty effect. Fig. 10 shows the as a function of , for one hadronic interaction models combination, for different cases in which and were shifted by a factor and 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 and have an insignificant effect on . In Fig. 10 right panel it is shown the effect of systematic shifts of 2.5% on and in opposite directions. The magnitude of the 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 dependence on primary mass and energy

For the following analysis, the range was defined to preserve a good statistics in all bins and it goes from to g/cm divided in bins of g/cm. The upper bound of 600 g/cm is highly influenced by the shower zenith angle limitation at due to the buried detector features.

Fig. 11 shows the evolution of as a function of 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 is seen in the bottom panel of Fig. 11, where are the differences with relation to the average value of for the three energy ranges considered. The differences of in all energy intervals is smaller than 1% for the entire range. The energy independence of is expected, since and evolve similarly with energy in the range from to 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 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 . The observed primary mass dependence is below 2%. The dependence of with the primary particle was minimized by choosing . The lack of 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 parameter

In this section, the capacity to constrain hadronic interaction models by measuring is demonstrated by studying the evolution of 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 years of operation. The infill array consists of m spaced water-Cherenkov stations spread over an area of km. 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 , the total exposure of this experimental set-up is km.sr.yr. Taking into account the cosmic-rays flux between and eV, the expected number of events to be measured per year is 1806.

Figs. 14 and 13 show the as a function of 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 . The worst case for the detector resolution was considered: % for and % for . Figs. 14 and 13 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 , the commonly used Merit Factor can be used. It is defined as:

(4) |

where and refer to any two hadronic interaction model combination and the ’s are the standard deviations of .

Fig. 15 shows the Merit Factor as a function of 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 between EPOS-LHC and QGSJetII-04 and between FLUKA and GHEISHA.

Figs. 17 and 16 show the Merit Factor as a function of the number of events and detector resolution for one particular bin: . The resolutions on and were re-scaled by a common factor , being that and . The effect of the number of events was calculated by re-scaling the standard deviation of the muon densities by a factor , where is the number of simulated showers and is the number of showers in each case.

Fig. 16 shows that it is possible to reach large Merit Factor values () for the separation between EPOS-LHC and QGSJetII-04/Sibyll2.3 using a reasonably small number of events () considering realistic detector resolutions ( and ). 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 . The same conclusions can be drawn about the low energy hadronic interactions models from Fig. 17. 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 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 , is sensitive to the model combination and presents a strong evolution with .

However, is not an easy parameter to be measured. Therefore in Sec. 4 a new and experimentally motivated parameter () is proposed and its correlation with is shown. 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 under realistic experimental limitations. Sec. 5 shows that the discrimination power of is significantly large. EPOS-LHC can be separated from Sibyll2.3 and QGSJetII-04 with large Merit Factor () using a reasonably small number of events (). Sibyll2.3 and QGSJetII-04 show similar average muon energy at ground and therefore can, in the best case, be discriminated with Merit Factor with about events. As for the low energy hadronic models, FLUKA can be separated from GHEISHA and UrQMD with large Merit Factor () using a reasonably small number of events (). GHEISHA and UrQMD can, in the best case, be discriminated with Merit Factor with about events. It was shown that correlated systematic uncertainties on and have insignificant influence on . Uncorrelated oppositive systematics of order of 2.5% can generate a systematic uncertainty on of the same order of the separation betweem the hadronic interaction models. This means that for a realistic application of 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 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 can in a short period of time contribute to the solution of the know problems with muon production in extensive air shower Aab2015a (); MuonsPRL ().

## 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).

## References

## References

- (1) K.-H. Kampert, M. Unger, Measurements of the cosmic ray composition with air shower experiments, Astropart. Phys. 35 (10) (2012) 660–678. arXiv:1201.0018, doi:10.1016/j.astropartphys.2012.02.004.
- (2) R. Engel, D. Heck, T. Pierog, Extensive Air Showers and Hadronic Interactions at High Energy, Annu. Rev. Nucl. Part. Sci. 61 (1) (2011) 467–489. doi:10.1146/annurev.nucl.012809.104544.
- (3) D. Heck et al., CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers, Tech. Rep. FZKA 6019, Forschungszentrum Karlsruhe GmbH, Karlsruhe (1998).
- (4) S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: QGSJET-II model, Phys. Rev. D 83 (1) (2011) 014018. arXiv:1010.1869, doi:10.1103/PhysRevD.83.014018.
- (5) K. Werner, T. Pierog, Extended Air Shower Simulations Using EPOS, in: AIP Conference Proceedings, Vol. 928, AIP, 2007, pp. 111–117. arXiv:0707.3330, doi:10.1063/1.2775903.
- (6) T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C 92 (2015) 034906. arXiv:1306.0121, doi:10.1103/PhysRevC.92.034906.
- (7) F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, T. Stanev, A new version of the event generator Sibyll, in: Proceedings of the 34th International Cosmic Ray Conference (ICRC), The Hague, Netherlands, 2015, p. 558. arXiv:1510.00568.
- (8) A. Ferrari, P. R. Sala, A. Fasso, J. Ranft, FLUKA: a multi-particle transport code, Tech. Rep. CERN-2005-10 (2005), INFN/TC_05/11, SLAC-R-773.
- (9) H. C. Fesefeld, Simulation of hadronic showers, physics and applications, Tech. Rep. PITHA 85-02, RWTH, Aachen, Germany (1985).
- (10) M. Bleicher, E. Zabrodin, C. Spieles, S. A. Bass, C. Ernst, S. Soff, L. Bravina, M. Belkacem, H. Weber, H. Stöcker, W. Greiner, Relativistic hadron-hadron collisions in the ultra-relativistic quantum molecular dynamics model, J. Phys. G 25 (1999) 1859. arXiv:hep-ph/9909407, doi:10.1088/0954-3899/25/9/308.
- (11) C. Meurer, J. Bluemer, R. Engel, A. Haungs, M. Roth, Muon production in extensive air showers and its relation to hadronic interactions, Czech. J. Phys. 56 (1) (2006) A211–A219. arXiv:astro-ph/0512536, doi:10.1007/s10582-006-0156-9.
- (12) R. Engel, for the Pierre Auger Collaboration, Upgrade of the Pierre Auger Observatory (AugerPrime), in: Proceedings of the 34th International Cosmic Ray Conference (ICRC), The Hague, Netherlands, 2015, p. 686. arXiv:1509.03732.
- (13) A. Aab et al. (The Pierre Auger Collaboration), The Pierre Auger Observatory Upgrade - Preliminary Design Report arXiv:1604.03637.
- (14) J. G. Gonzalez, R. Engel, M. Roth, Mass composition sensitivity of combined arrays of water cherenkov and scintillation detectors in the EeV range, Astropart. Phys. 74 (2016) 37–46. doi:10.1016/j.astropartphys.2015.09.003.
- (15) F. Sanchez, for the Pierre Auger Collaboration, in: Proceedings of the 32nd International Cosmic Rays Conference (ICRC), Beijing, China, 2011. arXiv:0906.4007.
- (16) A. Aab et al. (The Pierre Auger Collaboration), Prototype muon detectors for the AMIGA component of the Pierre Auger Observatory, J. Instrum. 11 (02) (2016) P02012. arXiv:1605.01625, doi:10.1088/1748-0221/11/02/P02012.
- (17) A. D. Supanitsky, A. Etchegoyen, I. Allekotte, M. G. Berisso, M. C. Medina, Underground muon counters as a tool for composition analyses, Astropart. Phys. 29 (2008) 461–470. arXiv:0804.1068, doi:10.1016/j.astropartphys.2008.05.003.
- (18) D. Ravignani, A. D. Supanitsky, A new method for reconstructing the muon lateral distribution with an array of segmented counters, Astropart. Phys. 65 (2015) 1–10. arXiv:1411.7649, doi:10.1016/j.astropartphys.2014.11.007.
- (19) D. Ravignani, A. D. Supanitsky, D. Melo, Reconstruction of air shower muon densities using segmented counters with time resolution, Astropart. Phys. 82 (2016) 108–116. arXiv:1606.02531, doi:10.1016/j.astropartphys.2016.06.001.
- (20) A. Aab et al. (The Pierre Auger Collaboration), Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events, Phys. Rev. D 91 (2015) 032003. arXiv:1408.1421, doi:10.1103/PhysRevD.91.032003.
- (21) A. Aab et al. (The Pierre Auger Collaboration), Testing hadronic interactions at ultrahigh energies with air showers measured by the pierre auger observatory, Phys. Rev. Lett. 117 (19) (2016) 192001. arXiv:1610.08509, doi:10.1103/PhysRevLett.117.192001.