A measurement of the muon number in showers using inclined events detected at the Pierre Auger Observatory

The average muon content of measured showers with zenith angles between 62◦ and 80◦ detected at the Pierre Auger Observatory is obtained as a function of shower energy using a reconstruction method specifically designed for inclined showers and the hybrid character of the detector. The reconstruction of inclined showers relies on a comparison between the measured signals at ground and reference patterns at ground level from which an overall normalization factor is obtained. Since inclined showers are dominated by muons this factor gives the relative muon size. It can be calibrated using a subsample of showers simultaneously recorded with the fluorescence detector (FD) and the surface detector (SD) which provides an independent calorimetric measurement of the energy. The muon size obtained for each shower becomes a measurement of the relative number of muons with respect to the reference distributions. The precision of the measurement is assessed using simulated events which are reconstructed using exactly the same procedure. We compare the relative number of muons versus energy as obtained to simulations. Proton simulations with QGSJETII show a factor of 2.13 ± 0.04(stat) ± 0.11(sys) at 1019 eV without significant variations in the energy range explored between 4 × 1018 eV to 7 × 1019 eV. We find that none of the current shower models, neither for proton nor for iron primaries, are able to predict as many muons as are observed.


INTRODUCTION
The Pierre Auger Observatory [1] is a hybrid cosmic ray air shower experiment that uses different techniques to detect extensive air-showers.The longitudinal shower development is observed by the Fluorescence Detector (FD) that provides a nearly calorimetric measurement of the shower energy.The duty cycle of the FD is approximately 13% [2].An array of over 1600 water-Cherenkov detectors that covers an area of 3000 km 2 with a 1.5 km triangular grid, samples the signal, S, as the air shower arrives at the Earth's surface.The Surface Detector (SD) has a 100% duty cycle.The FD detector surrounds the array with four buildings with six telescopes each.It is possible to relate the energy measured with the FD to the shower size which can be obtained with the SD, using events that can be reconstructed with both the SD and FD.This procedure provides the energy calibration [3].With the SD detector of the Pierre Auger Observatory we have recorded 5936 events between 1 January 2004 and 31 December 2010, with energy above 4 × 10 18 eV and zenith angle range 62 • < < 80 • .Those events represent a 25% increase over the exposure obtained with events < 60 • .The analysis of inclined events is important because it increases the sky coverage allowing the study of clustering and anisotropy in a larger region of the sky than it is possible with more vertical events.The dominant particles in these events are muons because most of the electromagnetic component is absorbed in the atmosphere.As the total number of muons depends on primary particle type, inclined showers can be used to study composition.Inclined EPJ Web of Conferences events are also important because they constitute the background in the search for neutrino-induced showers.In this paper we will describe the procedure used to reconstruct the shower size, N 19 , of inclined showers.We will explain how N 19 is calibrated using the energy measured by FD in events detected simultaneously by the SD and FD.Finally we will compare the behavior of the number of muons versus energy as observed with data to that obtained through Monte Carlo (MC) simulations.

RECONSTRUCTION OF THE SD SIZE PARAMETER N 19
The arrival direction of a cosmic ray is reconstructed from the signal arrival times by fitting a shower front model [4,5].The angular resolution achieved is better than 1 • at energies E > 4 × 10 18 eV.The reconstruction of the shower size requires techniques that are different to those used for more vertical showers.A procedure has been developed in which a set of muon densities at the ground, derived from simulations at different energies and zenith angles, is compared with experimental data.This technique was first used to analyse inclined showers recorded at Haverah Park [6] and has been subsequently adapted for the Auger Observatory [7] which uses water-Cherenkov detectors of exactly the same depth.Two reconstruction procedures have been developed independently.Using the two procedures provides an opportunity to test for systematic uncertainties arising from the different reconstruction processes.
The shower core position and the size parameter are obtained through a fit of the expected number of muons at each detector, , to the measured signal.The expected number of muons can be obtained from: where (x c ,y c ) are the coordinates of the shower core, and are the zenith and the azimuth of the shower direction, 19 is the reference muon distribution (the muon map) corresponding to the inferred arrival direction, A ⊥ ( ) is the detector area projected onto the shower plane.It has been shown that the dependence of the shape of 19 on energy and primary composition can be approximately factorised out and hence the two dimensional distributions of the muon patterns at ground level depend primarily on zenith and azimuthal angles because of the geomagnetic deviation of muons [8].N 19 is thus defined as the ratio of the total number of muons, N , in the shower with respect to the total number of muons at E = 10 EeV given by the reference distribution, N 19 = N (E, )/N map (E = 10 EeV, ).
Two sets of reference distributions have been produced for protons at E = 10 EeV.They are based on two different shower simulation packages and two different hadronic interaction models.One of them has been obtained with a set of histograms in two dimensions based on the model developed in [7], using simulations made with AIRES 2.6.0 [9] with QGSJET01 [10] as the choice for the hadronic interaction.The other set is a continuous parameterizations of the patterns obtained directly from the simulations [11] made with CORSIKA 6.72 [12] and the models QGSJETII [13] and FLUKA [14].Both approaches are valid out to 4 km from the shower axis in the energy range 10 18 eV to 10 20 eV and zenith angle 60 • to 88 • .The Efit (HasOffline) package uses the first (second) set based on QGSJET01 (QGSJETII).In Figure 1  From the measured signal S we obtain the muonic signal S subtracting the average electromagnetic (EM) component.We calculate the EM component in MC simulations in which we obtain r EM = S EM /S , the ratio of the EM to muonic signal.In the left panel of Figure 2 we show some examples of the dependence with the zenith angle a the distance to the core.r EM depends on composition, interaction models and shower energy.We have adopted proton AIRES simulations at 10 19 eV with QGSJET01 as a reference for this work.For inclined showers r EM is largest near the shower axis and its effect practically disappears at zenith angles exceeding about 65 • .Once r EM is parameterised the muonic signal is obtained as S = S/(1 + r EM ).
A maximum likelihood method is used to fit the shower size and the core position.This requires knowledge of the probability density function of the measured signal for all the triggered stations.The probability density p to observe a muonic signal, S , in a given station is obtained from the expected ( Here P oisson (k; ) gives the Poisson probability that k muons go through the detector when are expected, P T r (S ) is the probability of triggering and P st (k, S , ) is the probability density function (p.d.f.) of the muon signal S for k through-going muons.The response of the SD stations to the passage of a single muon is obtained from a detailed simulation using the standard GEANT4 package [16] within the official software framework of the Observatory [17].The p.d.f. for k muons, P st (k, S , ), is obtained by convolution.It is important to include details of the non-triggered stations to avoid biases in the reconstruction arising from upward fluctuations of the trigger.In the right panel of Figure 2 we show some example of p.d.f.distributions as a function of the number of muons.These stations provide upper limits to the signals.This can be achieved within the maximum likelihood method, replacing the 07003-p.3probability density by the probability that the detectors do not trigger:

EPJ Web of Conferences
The likelihood function to be maximised is then the product of these probabilities for all stations.In practice, we only use stations that lie within 5 km of the barycenter.The likelihood function depends on three free parameters that have to be fitted, the shower core position (x c ,y c ) and the shower size N 19 through the expected number of muons ( ) in Eq. ( 1) which determines the Poisson probabilities.The Poisson distribution takes into account fluctuations in the number of muons entering the detector, while the rest of the fluctuations in the station are introduced through the detector response.

ACCURACY OF N 19 RECONSTRUCTION USING MONTE CARLO SIMULATIONS
A sample of 100,000 proton showers were generated using AIRES and the hadronic model QGSJET01, with an E −2.6 energy spectrum in the range log 10 (E/eV) = (18.5, 20).Showers were chosen from a zenith angle distribution that is flat in sin 2 in the range (50 • , 89 • ), and uniform in azimuthal angle.A further 2700 proton showers were generated using CORSIKA, with the hadronic model QGSJETII and FLUKA.The energy spectrum is flat in log 10 (E/eV) in the range (18,20), and the zenith angle has a flat distribution in sin 2 in the range (60 • ,86 • ) and a uniform distribution in azimuthal angle.All SD events were generated within the Offline framework [17] with random impact points on the ground.For simulated showers it is possible to obtain the true value of the size parameter, referred as N MC 19 .In the left panel of Figure 3 we show the difference between the reconstructed N 19 and N MC  19 for both reconstructions.It is plotted as a function of the energy obtained by converting the value of N 19 into energy according to the calibration data (see figure 4 and related discussions).We note that the average bias for both reconstructions is below ∼ 4% level in the range log 10 (E/eV) = (18.5,19.7).We have found agreement between both reconstructions at the 3% level.In the same figure we show the energy resolution going from 25% at log 10 (E/eV) = 18.5 to 12% at high energies where we expect a more accurate reconstruction as there are more triggered stations.
In the right panel of figure 3 we show the average bias in the angular reconstruction as a function of the zenith angle input to the MC calculation.The zenith and azimuth angles are reconstructed with a bias of less than 0.05 • , and the opening angle, the angle between the MC and the reconstructed directions, is always less than 0.5 • .

CALIBRATION OF N 19 : THE TOTAL NUMBER OF MUONS
The energy of each event is obtained by calibrating N 19 with a set of quality hybrid events.In addition the calibration procedure can be used to obtain the number of muons as a function of energy, which is sensitive to cosmic ray composition and to the hadronic interactions in the shower.A fit is done using a power law A(E F D /10 EeV) B for energies above 4 × 10 18 eV where the array is 100% efficient.This considerably reduces the systematic uncertainties as the only use of hadronic models comes through the estimate of the energy carried by muons and neutrinos into the ground, the missing or invisible energy.The results of the fit are shown in the left panel of Figure 4. From the fit we obtain A = (2.13 ± 0.04 ± 0.11 (sys.)) and B = (0.95 ± 0.02 ± 0.03(sys.)).Details of the calibration procedure and the data are given respectively in [18] and [3].A spectrum has been obtained with these data [18].It is possible to relate N 19 and shower energy using MC simulation in an analogous way leading to A MC and B MC .The results depend on the choices made for composition and hadronic model.In the left panel of Figure 4 we also show the extreme cases for protons with QGSJETII and for iron nuclei with EPOS1.99 [19].
Using the formula for the calibration fit [18], it is also possible to derive the muon number in the data compared to the predictions from the simulation N data 19 /N MC 19 as a function of E F D .Simple calculations yield the following expression: We show the observed number of muons in data compared with the number of muons obtained with the two extreme predictions for proton QGSJETII (solid line) and iron EPOS1.99 (dashed line) in the right panel of EPJ Web of Conferences for proton (iron) QGSJETII (EPOS1.99) is 1 (1.7) and that of B MC is 0.934 (0.928).The number of muons deduced from data exceeds that of proton QGSJETII simulations by a factor of 2.1 and that of iron EPOS1.99 by 23%.Withing the errors no significant dependence on the energy is obtained in either case.

SUMMARY AND DISCUSSION
The reconstruction method for inclined showers has been tested with simulated data.It has been shown that the average opening angle between the true and reconstructed directions is below 0.5 • and that the average shower size obtained in the reconstruction reproduces the simulated one within 4%.We have shown that the reconstruction of inclined events can be used to extract the number of muons from the data.This is done through the energy calibration that relates the FD energy to N 19 , the muon size with respect to a reference model.When compared to protons simulated with QGSJETII the ratio of the total number of muons at E F D = 10 EeV is measured to be (2.13 ± 0.04 ± 0.11 (sys.)).The 22% systematic uncertainty in the FD energy measurement [2] has not been included.Several other methods have been developed to obtain the number of muons from the data collected at the Pierre Auger Observatory.All of them use events below 60 • .These methods report similar enhancements of the muon content in the showers [20].
We find that none of the current shower models, neither for proton or for iron primaries, are able to predict as many muons as are observed.
we show as an example two reference muon distribution at 70 • and 88 • .

Figure 3 .
Figure 3. Left panel: the zenith (square) and azimuth (circles) accuracy as a function of the zenith angle input from MC.The triangles shows the opening angle resolution.Right panel: average (points) and RMS (lines) values of the relative difference between reconstructed (N 19 ) and simulated (N MC 19 ) shower sizes.The circles and solid line are the results for Efit.The squares and dashed lines for HasOffline.

Figure 4 .
Figure 4. Left panel: fit of the calibration curve N 19 = A(E/10 EeV) B. The constants A and B are obtained using the maximum-likelihood method.The contours indicate the constant levels of the p.d.f.integrated over zenith angle, corresponding to 10, 50 and 90% of the maximum value[18].Calibration curves for proton QGSJETII (dot line) and iron EPOS1.99 (dashed line) are shown for comparison.Right panel: ratio of the number of muons in data compared to proton QGSJETII (solid line) and iron EPOS1.99 (dashed line) as a function of the energy.The grey bands indicates the systematic uncertanties in N 19 .See text for details.

Figure 4 .
The grey bands corresponds to the N 19 systematic uncertainties obtained by comparison between reconstructed and the true values of N 19 using simulated data.The value of A MC 07003-p.5