Measurements of the muon shower content at the Pierre Auger Observatory

Several methods for estimating the muonic part of the signal observed in the surface Cherenkov detectors have been developed within the Pierre Auger Collaboration in the recent years. The muon shower content, derived from data with these methods, is found to be significantly larger in comparison with predictions of QGSJET II interaction model.


INTRODUCTION
The Pierre Auger Observatory is a hybrid detector designed for measurements of the characteristics of longitudinal and lateral profiles of extensive air showers (EAS) from ultra-high energy cosmic rays with a fluorescence detector (FD) and a surface array of Cherenkov water detectors (SD) [1]. Though no dedicated equipment is available for the measurement of the muon shower content, several methods for its estimation have been developed within the collaboration. In this work we present the results of the methods based on studies of the temporal structure of the total signal in surface detectors (smoothing method and multivariate muon counter), on EAS universality and on fitting of hybrid events [2][3][4][5]. In addition, we give results for inclined showers (see also [6]) where the muon shower content can be measured almost directly with surface detectors since such showers are practically free from contamination with electromagnetic particles.

A smoothing method
The electromagnetic (EM) contribution to the signal in the surface detectors can be estimated by a smoothing method [4]. In fact, due to the difference in the energy deposition in a water tank between muons (≈ 240 MeV) and electrons or photons ( 10 MeV), muons appear as sudden variations of the signal over a smoother electromagnetic background. Similar increases of the electromagnetic signal can be produced only by high energy photons: above 300 MeV their fraction with respect to the number of true muonic signals is estimated to be less than 10% and is included in the total systematic uncertainty.
A smoothing filter averages the total trace recorded by the FADCs of the 3 PMTs in each SD station over a preset number of consecutive time bins N bin . Any positive difference between the original trace a e-mail: yushkov.alexey@gmail.com b For the full authorlist see Appendix "Collaborations" in this volume. This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. EPJ Web of Conferences and the smoothed signal is assigned to the muon component and subtracted from the signal; then the whole procedure is applied again for a number of iterations N iter .
To avoid severe under-or over-smoothing, a variable smoothing window has been used, according to the functional dependence of the muon signal duration on zenith angle. With this choice, the bias in the evaluation of the EM component for both proton and iron primaries, and for the largest possible angular range, is minimized after a small number of iterations. A residual small dependence on energy, less than 10% in the considered energy range, is removed by using a Monte Carlo based correction factor.
Using showers simulated with QGSJET II [7] for E > 3 EeV (corresponding to full efficiency of the SD) up to 60 • , the EM contribution to the signals has been evaluated for distances between 700 and 1500 m from the core with a resolution of 17% and a systematic uncertainty below 8%, irrespectively of the primary energy and composition. The independence of the results on the hadronic interaction models has been checked on EPOS 1.6 [8].
Having evaluated the EM component, the muon component in each detector can then be derived by subtraction with systematic uncertainties below 8% and a resolution close to 20%.

A multivariate muon counter
In this section, the reconstruction of the muon number at 1000 m from the shower axis is accomplished by first estimating the number of muons in the surface detectors using the characteristic signals created by muons in the PMT FADC traces and then reconstructing the muonic lateral distribution function (LDF) of SD events.
In the first stage, the number of muons in individual surface detectors is estimated. As in the jump method [4], the total signal from discrete jumps was extracted from each FADC signal, where x i is the signal measured in the ith bin in vertical equivalent muon (VEM) units, and the indicator function I {y} is 1 if its argument y is true and 0 otherwise. The estimator J is correlated with the number of muons, but it has an RMS of approximately 40%. To improve the precision, a multivariate model was used to predict the ratio = (N + 1)/(J + 1). 172 observables that are plausibly correlated to muon content, such as the number of jumps and the rise-time, were extracted from each FADC signal. Principal Component Analysis was then applied to determine 19 linear combinations of the observables which best capture the variance of the original FADC signals. Using these 19 linear combinations, an artificial neural network (ANN) [9] was trained to predict and its uncertainty. The output of the ANN was compiled into a probability table P ANN = P (N = N | FADC signal). The RMS of this estimator is about 25%, and biases are also reduced compared to the estimator J .
In the second stage of the reconstruction, a LDF is fit to the estimated number of muons in the detectors for each event, where r is the distance of the detector from the shower axis and ν, , and are fit parameters. The number of muons in each surface detector varies from the LDF according to the estimate P ANN and Poisson fluctuations. The fit parameters, ν, , and , have means which depend on the primary energy and zenith angle as well as variances arising from shower-to-shower fluctuations. Gaussian prior distributions with energy-and zenith-dependent means were defined for the three fit parameters. All the parameters were estimated using an empirical Bayesian approach: three iterations were performed between (i) finding the maximum a posteriori estimate ν i , i , and i for each shower i given the fixed priors, and (ii) re-estimating the priors given the fixed parameter estimates ν i , i , and i .

UHECR 2012
The value of the muonic LDF at 1000 m, exp( ν), is highly correlated with N (1000), the number of muons in surface detectors at 1000 m from the shower axis. The RMS of exp( ν) in showers simulated using QGSJET II is 12% and 5% for proton and iron primaries. To correct several biases that depend on the energy and zenith angle of the showers, a quadratic function f (exp( ν), ) was tuned on a library of showers simulated using QGSJET II with simulated detector response. The final estimator N (1000) = f (exp( ν), ) has a systematic uncertainty below 50 • of 6% from uncertainty in the composition and +10% −0% from uncertainty in the hadronic models, determined by reconstructing showers simulated using EPOS 1.6. The total systematic uncertainty decreases with the zenith angle: at = 55 • it is +9% −3% .

Inclined showers
The shower size N 19 of inclined showers is expressed with respect to a size of a reference shower conventionally chosen to be that from proton at 10 19 eV [10]. Since inclined showers are essentially free from electromagnetic contamination from the hadronic cascade, N 19 is proportional to the total number of muons n . The shower size N 19 was calibrated to the FD energy E FD using the power law relation A (E FD /10 EeV) B for energies above 4 × 10 18 eV, where the SD is 100% efficient. Details of the calibration procedure can be found elsewhere [6,10]. Calibrating both data and simulations, one readily obtains a ratio between number of muons in them: The best fit gives the following parameters: A = 2.13 ± 0.04(stat) ± 0.11(sys) and B = 0.95 ± 0.02(stat) ± 0.03(sys). For simulations of proton showers with QGSJET II we obtain A MC = 1 and B MC = 0.93. A can be interpreted as the measured number of muons relative to the reference one, QGSJET II protons in the given case, used for the reconstruction at the FD energy of 10 EeV. More details on determination of relative number of muons in inclined showers can be found in a dedicated paper in these proceedings [6].

Universality of S /S em behavior on X v max
The ratio of the muonic signal to the electromagnetic (EM) signal, S /S em , at 1000 m from the shower axis exhibits an empirical universal property for all showers at a fixed vertical depth of shower maximum, X v max [11]. S /S em is independent of the primary particle type, primary energy, and zenith angle. The dependence of S /S em on X v max can be described by a simple parameterization which leads to the following expression for the muonic signal where S (1000) is the total signal at 1000 m, is the zenith angle, , A, a and b are parameters [12]. The estimation of the muonic signal in data is complicated by the dependence of the fit parameters in Eq. (4) on the choice of hadronic interaction model. This dependence gives rise to a systematic uncertainty in the measurement of the muonic signal in data which is difficult to quantify due to the uncertainty in properties of hadronic interactions. For zenith angles above 45 • this uncertainty becomes small [13], since S /S em increasingly reflects the equilibrium between muons and EM halo from their decays and interactions, becoming less depending on the properties of particular hadronic interaction code. The fit in Eq. (4) provides an unbiased estimate of both the muonic and EM signals. The RMS of the difference between simulated muonic signal and the parameterization (4) in CORSIKA [14] showers is less than 10% and 5% for proton and iron primaries. The systematic uncertainty of S from the uncertainty in the hadronic models is estimated to be 10% below 45 • and quickly diminishes to 1-2% 07002-p.3  for zenith angles above 55 • . The full systematic uncertainty of S after event reconstruction is 14% at 10 19 eV, determined through complete shower simulation and reconstruction using the Auger Off line software package [15,16]. The systematic uncertainty from event reconstruction is dominated by the systematic uncertainty of S (1000), with only a few percent coming from the uncertainty of X v max and the zenith angle.

Application to data
By applying the multivariate counter, the smoothing method, and the universality property to data collected between 1 January 2004 and 30 September 2010, a significant deficit of muons is found in comparison to simulations with QGSJET II (see Fig. 1). The multivariate method was applied to SD events over the energy range 18.6 < log(E) < 19.4 and 0 • -57 • . The smoothing method was applied to SD events over the energy range 18.7 < log(E) < 19.2 and 0 • -60 • . The universality method was applied to hybrid events over the energy range log(E) > 18.6 and 20 • -65 • . For all methods, the relative number of muons in data N rel is estimated with respect to simulated showers for proton primaries.
The multivariate method is used to determine the number of muons N (1000), while the smoothing finds it by difference, subtracting the electromagnetic component from the total signal. In both cases, the relative number of muons in data is equal to 1.6 − 1.7 when compared to the prediction of QGSJET II for protons and is almost angle-independent until about 40 • , above which it increases to 1.9 − 2.0. A similar conclusion can be reached from applying the universality. In inclined showers the deficit of muons in simulations is somewhat larger, but compatible with the results from the other methods.

STUDY OF INDIVIDUAL HYBRID EVENTS
At the Auger Observatory, thousands of showers have been recorded for which reconstruction has been possible using both the FD and SD. These hybrid events have been used to construct a library of simulated air-shower events where the longitudinal profile (LP) of each simulated event matches a measured LP. The measured LP constrains the natural shower-to-shower fluctuations of the distribution of particles at ground. This allows the ground signals of simulated events to be compared to the ground signals of measured events on an event-by-event basis.
Hybrid events were selected using the criteria adopted for the energy calibration of the SD [17] in the energy range 18.8 < log(E/eV) < 19. 2008. 227 events passed all cuts. Air showers were simulated using SENECA [18] with QGSJET II and FLUKA [19] as the high-and low-energy event generators. For every hybrid event, three protonand three iron-initiated showers were selected from a set of 200 simulated showers for each primary type. The energy and arrival direction of the showers match the measured event, and the LPs of the selected showers have the lowest 2 compared to the measured LP. The measured LP and two selected LPs of an example event are shown in the left panel of Figure 2. The detector response for the selected showers was simulated using Off line . The lateral distribution function of an observed event and that of two simulated events are shown in the right panel of Figure 2. For each of the 227 events, the ground signal at 1000 m from the shower axis, S (1000), is smaller for the simulated events than that measured. The mean ratio of the measured S (1000) to that predicted in simulations of showers with proton primaries, S(1000) Data S(1000) Sim , is 1.5 for vertical showers and grows to around 2 for inclined events. The ground signal of more-inclined events is muon-dominated. Therefore, the increase of the discrepancy with zenith angle suggests that there is a deficit of muons in the simulated showers compared to the data. The discrepancy exists for simulations of showers with iron primaries as well, which means that the ground signal cannot be explained only through composition.

DISCUSSION
As demonstrated in the analyses, and shown in Figs. 1, 2, simulations of air showers using QGSJET II with proton and iron primaries underestimate both the total detector signal at ground level and the number of muons in events collected at the Pierre Auger Observatory. These discrepancies could be related to an incorrect energy assignment within the 22% systematic uncertainty of the energy scale of the Pierre Auger Observatory and shortcomings in the simulation of the hadronic and muonic shower components.
To explore these potential sources of discrepancy, a simple modification of the ground signals was implemented in the simulated hybrid events of Section 3. The uncertainty in the energy scale motivates the rescaling of the total ground signal by a factor E Resc , and the muon deficit motivates a rescaling of the signal from hadronically produced muons by a factor N ,Resc . The rescaled muonic and EM components of S (1000), S and S EM -both defined for proton primaries -modify the ground signal where the exponent 0.92 is the energy scaling of the muonic signal predicted by simulations. The rescaling factors were applied uniformly to all events. This represents a simplistic modification, and N ,Resc does not reflect any changes in the attenuation (i.e. zenith angle dependence) and LDF of muons.  However, both the attenuation and LDF would change if, for example, the energy spectrum of muons predicted by simulations is not in agreement with the data. E Resc and N ,Resc were determined simultaneously by making a maximum-likelihood fit between the modified, simulated S (1000) and the measured S (1000) for the ensemble of hybrid events. The best fit values of N ,Resc and E Resc are (2.21 ± 0.23 (stat.) +0. 18 −0.23 (syst.)) and (1.09 ± 0.08 +0.08 −0.06 ) respectively; see Figure 3. The systematic uncertainties arise from uncertainty in the composition and event reconstruction.
The signal rescaling in simulated hybrid events is fundamentally different from the other methods. The observational muon enhancement, N rel , which includes all muons, cannot be compared directly to N ,Resc , which represents an increase of only the hadronically produced muons and their decay products. In addition, the potential increase of N rel with zenith angle suggests that a global rescaling of the ground signal from muons is overly simplistic.
In summary, all of the analyses show a significant deficit in the number of muons predicted by simulations using QGSJET II with proton primaries compared to data. Various factors can contribute to this discrepancy and in part it can be related with the current 22% uncertainty on FD energy. The other two evident factors are shortcomings in the hadronic models and the presence of nuclei in primary cosmic rays that are heavier than protons. These factors can be helpful in explaining the dependence of muon deficit on the zenith angle. The average energy of muons at 1000 meters from the core grows from around 1 GeV at 20 • to 5 GeV at 65 • . The observed increase of N rel with zenith angle thus can be interpreted as an indication of harder muon spectra in data and can explain also the higher value of N rel found for inclined showers. Formally, the muon deficit is almost completely negligible for EPOS 1.99 iron showers, in which the signal is 70-80% larger than in proton showers of QGSJET II, but such a heavy composition is in clear contradiction with the measurements of the depth of shower maximum [20] and does not provide an explanation to the growth of N rel with zenith angle.