Measurement of the muon content in air showers at the Pierre Auger Observatory

The muon content of extensive air showers produced by ultra-high energy cosmic rays is an observable sensitive to the composition of primary particles and to the properties of hadronic interactions governing the evolution of air-shower cascades. We present different methods for estimation of the number of muons at the ground and the muon production depth. These methods use measurements of the longitudinal, lateral, and temporal distribution of particles in air showers recorded by the detectors of the Pierre Auger Observatory. The results, obtained at about 140 TeV center-of-mass energy for proton primaries, are compared to the predictions of LHC-tuned hadronic-interaction models used in simulations with different primary masses. The models exhibit a deficit in the predicted muon content. The combination of these results with other independent mass composition analyses, such as those involving the depth of shower maximum observable Xmax, provide additional constraints on hadronic-interaction models for energies beyond the reach of the LHC.


Introduction
Within a scheme of the generalized Heitler model of hadronic air showers, muons are produced by decays of charged mesons [1].In the evolution of the cascade, energies of charged pions are reduced to a few dozens of GeV [2], at which point the decay length becomes shorter than the interaction length, and, therefore, the decay of pions into muons starts to dominate.Muon production in air showers depends on intricate details of the hadronic interactions; nevertheless, in the simple Heitler scheme, the number of produced muons N μ approximately follows a power law where E and A are the primary cosmic-ray energy and mass, respectively, c is the critical energy at which charged pions decay, and the value for β, which is approximately 0.9 for a fixed value of A (for fixed composition as well), is obtained from the multiplicity and fraction of charged pions.Measurements of the muon content of extensive air showers produced by ultra-high energy cosmic rays thus give excellent insight into the physics of hadronic interactions and multi-particle productions a e-mail: darko.veberic@kit.eduLeft: Energy-normalized relative muon number R μ /(E/10 19 eV) for various interaction models and two primaries, p and Fe, as obtained from inclined events [5].Note that the Auger data (black points) indicate a heavy composition and exhibit a different slope (black line) than the model predictions (dashed and dotted lines).
Right: ln R μ vs. X max for various interaction models (lines) and primary masses (markers) and the Auger data (black point with systematic-uncertainty brackets) from events with energies around E = 10 19 eV [5].
that take place in such high-energy environments.When the primary energy and/or mass are obtained from independent observables, muon measurements can be used for direct tests of existing hadronicinteraction models and as such, these tests deliver the best constraints on model parameters.Furthermore, the dependence of muon numbers on the primary mass can be used as a sensitive estimator of cosmic-ray composition.
At the Pierre Auger Observatory we are performing extensive muon studies, and the planned upgrade of the detectors [3] is primarily aimed at further increasing the separability of the muonic and electromagnetic components of the showers.The two main parts of the Observatory [4] are the Surface Detector (SD), consisting of a ∼3000 km 2 array of water-Cherenkov stations with an almost 100% duty cycle, and the Fluorescence Detector (FD) which is comprised of 27 telescopes operating on a ∼15% duty cycle.
Due to the limited space available in the following Sections, we only briefly outline the main features of muon studies based on data from the Pierre Auger Observatory.

Hybrid observations of inclined events
Since the electromagnetic component of a shower is absorbed faster in the atmosphere than the muonic part, highly inclined showers deliver an almost pure muon "beam", which is suitable for direct studies.Here we use hybrid events, which were recorded by both the SD and FD, and with zenith angles in the range 62 • < θ < 80 • .While the whole procedure is more complicated [5] and involves construction of an unbiased estimator, it can be schematically described as a method where, for each event, the muon density at the ground ρ rec μ , which is reconstructed from measured SD signals, is separated into a footprint shape ρ map μ and a relative muon content R μ as , where the template ρ map μ is derived from simulations of protons with energy E = 10 19 eV and with one selected hadronic-interaction model.The template takes into account most of the asymmetric and elongated shape of the footprint at the ground and depends only weakly on the cosmic-ray energy and mass.
As illustrated in Eq. ( 1), we expect a power-law dependence of the mean muon content on energy, R μ = a(E/10 19 eV) b .A fit to 174 high-quality events gives a = 1.84 ± 0.03 ± 0.32(syst.)and b = 1.03 ± 0.02 ± 0.03(syst.).In Fig. 1-left, the energy-normalized R μ /(E/10 19 eV) is shown for simulations with different primary masses and interaction models, and we can see that the proton and iron showers are well separated.While the large systematic uncertainty (square brackets) limits its power as a mass-sensitive estimator, R μ still delivers valuable insights due to the discrepancy between its absolute scale in the data and the simulations.The diverging slopes between the data and simulated (pure) primaries also suggest a transition from lighter to heavier composition.The discrepancies with hadronic-interaction models become even more obvious when we combine muon results with those of our standard mass-sensitive method [6] involving the depth of shower maximum X max .In Fig. 1-right, simulation results at E = 10 19 eV are shown in a ln R μ vs. X max space together with the data point which lies far away from all model predictions for any primary mass.
The value obtained for ln R μ at E = 10 19 eV and the logarithmic slope d ln R μ /d ln E can be compared to the corresponding results from simulations.In Fig. 2, the results from Auger data (black circles) are compared to simulations performed with various hadronic-interaction models for samples of purely proton (red circles) and purely iron (blue squares) primaries and a mixed composition sample where ln A matches that obtained from conversions of X max (hexagons).In all cases, the models fail to match the measurements and would require substantial increases of muon production.Furthermore, the large value of the measured slope favors a changing composition.

Muon production depth
Under certain conditions (large zenith angles, large distances to the shower core), we can assume that the signal recorded in the SD stations is produced mostly by traversing muons.Through simple max obtained from muon production depth distributions using Auger data (black points) and simulations (lines) for proton and iron primaries [7].Center and right: X μ max converted into ln A values using the two interaction models (red points) shown together with ln A values obtained from conversions of X max measurements (black triangles) [7].geometrical considerations, the muon production distance z can be reconstructed from the time delay t g between the passage of the shower front and the occurrence of the muon signal as where r is the distance of the detector to the shower axis, Δ is the perpendicular distance of the SD station to the position of the shower plane at the time of muon occurrence, and z π is the average distance parent mesons travel before decaying (see [7] for details).For a shower with a zenith angle θ, the production distance z is then converted to the equivalent atmospheric depth with where ρ(h) is the atmospheric density at height h.In this way, signals detected by SD stations are converted into a muon production depth distribution for each event.The maximum, X μ max , of a muon production depth distribution is then obtained with a fit to a Gaisser-Hillas function.
In Fig. 3-left, the mean of the obtained X μ max values from 481 events is shown in several energy bins, together with simulation predictions from two hadronic-interaction models and two primaries [7].Again, the data indicates a change in composition as the energy increases and shows a flatter trend than for pure proton or for pure iron.The X μ max values can be converted to a mean logarithmic mass ln A for each hadronic-interaction model and the rest of Fig. 3 shows this conversion for the two models together with the standard conversions of X max .For both models, results indicate a transition to heavier primaries, but the conversion with the Epos-LHC model predicts implausibly heavy elements as primary cosmic rays.

Top-down reconstruction of hybrid events
In contrast to the inclined-shower case (as described in Sect.2), the estimation of the muon content is more complicated for events with zenith angles θ < 60 • due to the admixture of the electromagnetic component.Nevertheless, for 411 high-quality hybrid events in the energy range 18.8 < lg(E/eV) < 19.2, a so-called "top-down" procedure can be performed [8] where the longitudinal profile, as measured with the FD, is first reproduced by running many simulations and finding  The average ratio R of the SD signal at 1000 m in Auger data and the simulations that match the longitudinal profile measured by the FD [8].Note the increase of R with the increasing zenith angle θ.Right: The 1σ contours of statistical uncertainty for fits of R had and R E rescalings together with gray bands for systematic uncertainties [8].For both plots, shower simulations were performed with several hadronic-interaction models, for pure proton and mixed compositions.the best match.Such simulated events thus have an X max and energy compatible with the measured counterpart.In Fig. 4-left, the ratio of measured and simulated ground signals at 1000 m from the core, R = S (1000)/S MC (1000), is shown as a function of sec θ.The increase of this ratio towards the large zenith angles, where showers have less of the electromagnetic component, may indicate that the measured data at the ground contain more muons than in simulations.Independently of the models and simulated composition, the ratio R is not compatible with 1 within the FD energy-scale uncertainty and actually increases with zenith angle where the muon component of the shower starts to dominate.This indicates that, indeed, it is the muon production that is not reproduced well by the interaction models.
To identify the main contribution to the ground-signal discrepancy, the different components of the shower are, in a second step, self-consistently modified in simulations, allowing rescaling of the hadronic interactions R had and, independently, energy rescaling R E .Note that R had rescales not only the nominal muon component at the ground but also a part of the electromagnetic component originating from decays of muons and the part of the electromagnetic component from hadronic jets produced close to the detectors.In Fig. 4-right, the results of this rescaling analysis show that in order to reproduce the measured data, almost no energy rescaling R E is needed and that most of the discrepancy can be explained by muon (or "hadronic") deficiency in the models.Note that mixed in both plots again refers to a mixed composition which reproduces the X max measurement at this energy.

Temporal structure of signals
While all the methods mentioned above were using only total integrated signal and ignored the detailed temporal structure of the SD signal at the ground, application of different filtering techniques on recorded traces can produce a certain degree of separation between the electromagnetic and muonic contributions [9].
Estimating the electromagnetic contribution by a smoothing method and a multi-variate method involving many time-structure related observables shows in either case that relative to the proton simulations, the muon proportion of the signal is up to a third larger in the measured data.Nevertheless, these methods are less powerful due to larger sensitivity to energy-scale systematics and are biased by the aging effects observed in SD response to shower particles from the data that has been gathered in nearly 12 years of the Observatory's operation.

Summary
In these proceedings, we only briefly described efforts within the Pierre Auger Collaboration related to the analysis of muon production in extensive air showers and the influence of muons on energy calibration and composition estimation for primary cosmic-ray energies around E = 10 19 eV [5,[7][8][9].The numbers of muons in measured data relative to the air-shower simulations have been reported for small and large zenith angles, and in all cases indicate a non-negligible dearth in predictions by hadronic-interaction models.These discrepancies cannot easily be reconciled with measurements and even when they can be, then only marginally so.
On the other hand, there is no evidence of increased event-to-event fluctuations in ground signal for showers with fixed X max , which indicates that the observed muon deficit in models cannot be attributed to some exotic scenario (like micro black-holes) taking place in the early stages of shower development that would only increase muon numbers for a certain subset of events [8].
Nevertheless, an important effect can be expected from the hadronic-interaction models when they are tuned [10] to reproduce the rho-meson spectra from collider experiments [11] more accurately, especially in the extreme forward region, where more neutral rho-mesons can be found in the data [12].Due to their decays into charged pion mesons, the neutral rho-mesons can namely keep more energy in the hadronic part of a shower and therefore increase the muon production.

,Figure 1 .
Figure1.Left: Energy-normalized relative muon number R μ /(E/10 19 eV) for various interaction models and two primaries, p and Fe, as obtained from inclined events[5].Note that the Auger data (black points) indicate a heavy composition and exhibit a different slope (black line) than the model predictions (dashed and dotted lines).Right: ln R μ vs. X max for various interaction models (lines) and primary masses (markers) and the Auger data (black point with systematic-uncertainty brackets) from events with energies around E = 10 19 eV[5].

,Figure 3 .
Figure 3. Left: X μ max obtained from muon production depth distributions using Auger data (black points) and simulations (lines) for proton and iron primaries[7].Center and right: X μ max converted into ln A values using the two interaction models (red points) shown together with ln A values obtained from conversions of X max measurements (black triangles)[7].

Figure 4 .
Figure 4. Left: The average ratio R of the SD signal at 1000 m in Auger data and the simulations that match the longitudinal profile measured by the FD[8].Note the increase of R with the increasing zenith angle θ.Right: The 1σ contours of statistical uncertainty for fits of R had and R E rescalings together with gray bands for systematic uncertainties[8].For both plots, shower simulations were performed with several hadronic-interaction models, for pure proton and mixed compositions.