Direct Measurement of the Muon Density in Air Showers with the Pierre Auger Observatory

As part of the upgrade of the Pierre Auger Observatory, the Auger Muons and Infill for the Ground Array (AMIGA) underground muon detector extension will allow for direct muon measurements for showers falling into the SD-750 array. We optimized the AMIGA muon reconstruction procedure by introducing a geometrical correction for muons leaving a signal in multiple detector strips due to their inclined angle of incidence and deriving a new unbiased parametrization of the muon lateral distribution function. Furthermore, we defined a zenith-independent estimator ρ35 of the muon density by parametrizing the attenuation of the muonic signal due to the atmosphere and soil layer above the buried detectors and quantified the relevant systematic uncertainties for AMIGA. The analysis of one year of calibrated data recorded with the prototype array of AMIGA confirms the results of previous studies indicating a disagreement between the muon content in simulations and data.


Introduction
Measurements of the Pierre Auger Observatory have led to major advances in our understanding of ultra-high energy cosmic rays. Nevertheless, open questions regarding their origin and properties remain. To make progress in answering these questions, the Pierre Auger Observatory is currently being upgraded. The key part of the Auger-Prime upgrade is the installation of new plastic scintillation detectors on top of the existing water-Cherenkov detectors (WCDs) of the surface detector (SD) array. The combined analysis of the signals of both the WCDs and the new surface scintillator detectors (SSDs) will allow for the disentangling of the electromagnetic and muonic shower components and hence provide additional measurements of mass composition-sensitive observables up to energies of the flux suppression region [1].
In addition, 61 scintillation detectors with an area of 30 m 2 each will be buried at a depth of 2.3 m in the soil next to each of the WCDs of the SD-750 array. The so-called Auger Muons and Infill for the Ground Array (AMIGA) extension will provide direct measurements of the muon content of a sub-sample of extensive air showers falling onto the SD-750 array and serve for the verification and fine-tuning of the methods used to extract muon information from the combined SSD and WCD signals.
A prototype array, consisting of seven muon counters arranged in a hexagonal layout as shown in Fig. 1, * e-mail: sarah.mueller2@kit.edu * * e-mail: auger_spokespersons@fnal.gov * * * Authors list: http://www.auger.org/archive/authors_2018_10.html  has been operational since February 2015. In addition to five 30 m 2 muon counters, two twin counters of double size were installed to assess the muon counting accuracy. Each muon counter of the prototype array is either composed of three 10 m 2 scintillator modules or two 5 m 2 and two 10 m 2 modules that are segmented into 64 scintillator bars each [2]. While the muon counters of the prototype array have been equipped with photomultiplier tubes (PMTs), the final design of AMIGA foresees to re- place PMTs with silicon photomultipliers (SIPMs). Furthermore, each counter will consist of three 10 m 2 modules [3,4].
In this work, we report on recent improvements of the AMIGA muon reconstruction procedure, estimate the main systematic uncertainties and present first physics results for AMIGA that were obtained from the analysis of one year of data recorded with the prototype array.

Muon number reconstruction
The number of muons impinging the buried muon detectors is counted individually for each of the segmented scintillator modules. Here, we describe the counting strategy for the PMT design which was employed in the prototype phase. The counting strategy for the new SIPM design is explained in [4].
The scintillation light which is produced by a muon hitting a scintillator module is guided by wavelengthshifting (WLS) fibers to the central 64-pixel PMT. The amplified analog PMT signals are then digitized to binary time traces by means of adjustable threshold discriminators with subsequent field-programmable gate array (FPGA) sampling at a frequency of 320 MHz corresponding to 3.125 ns. The discriminator thresholds are calibrated individually for each channel by a background calibration method to ≈ 30% of the mean channel-specific single photo-electron (SPE) amplitude V SPE [2].
The muon counting strategy is based on the identification of patterns in the binary time traces for each scintillator bar within time windows of fixed length [5]. Analog muon pulses can have complex time structures leading to null samples in the binary traces after the pulse discrimination. They need to be distinguished from isolated SPE pulses originating from background effects like cross-talk between neighboring PMT pixels or thermal fluctuations. For this purpose, the muon identification pattern has been chosen as 1x1 with x ∈ [0, 1]. Since SPEs produce at most two consecutive positive samples (0110) in the bi-nary time trace for the chosen discriminator thresholds, matches of the 1x1 pattern efficiently reject background and, at the same time, identify muons with possible disconnected time traces.
To prevent that two or more 1x1 patterns are matched within the binary trace produced by one muon, an inhibition time window starting from the first identified positive sample (1) is applied over which the pattern matching process is stopped. It has been shown with simulations that the optimal window size is seven bins corresponding to 22 ns.
Since it is not possible to resolve further muons within the time interval of the applied inhibition window, a statistical pile-up correction for multiple muons hitting the same scintillator bar within the same time window is applied. For this purpose, the digital time traces of all channels are split into N windows with a length of the chosen inhibition window size. The number of channels k i with a muon pattern match starting within the ith window is counted for each window and the true number of muons is estimated byμ where N seg = 64 is the number of detector segments. The total number of muonsμ = N i=1μ i is obtained by summing over all N time windows [6].
In addition to the pile-up correction, a correction for over-counting resulting from corner-clipping muons is applied. These are muons that arrive from an inclined, nonvertical direction w.r.t. the surface of a scintillator module such that a signal is deposited in two neighboring scintillator bars. Besides the muon inclination angle θ w.r.t. the upwards pointing z-axis, the difference in azimuth ∆ϕ m = ϕ − ϕ m between the momentum direction of the muon and the orientation of the module in the ground plane is of crucial importance as illustrated in Fig. 2.
The relative corner-clipping bias b clip = (N Rec −N MC ) /N MC of the reconstructed number of muons N Rec w.r.t. the num- 14.3% ber of simulated muons N MC impinging a detector module has been parametrized with simulations as [5] where the mean muon zenith and azimuth angles are approximated by the angles θ SD and ϕ SD of the shower axis that are obtained from the SD geometry reconstruction. The average dependence of the relative over-counting on the azimuth difference ∆ϕ m is shown in Fig. 2. The found parametrization is used to correct the number of reconstructed muons by for each orientation of the underground detectors on an event-by-event basis.
When an extensive air shower falls onto the detector array, only a part of the shower is sampled by the muon detectors at discrete distances from the shower axis. As a proxy for the muon content of the shower, the muon density at an optimal reference distance r opt = 450 m [7,8] is estimated by fitting a muon lateral distribution function (MLDF) to the muon densities at detector level. For AMIGA, the chosen parametrization follows the KASCADE-Grande experiment [9,10]. The values of r * , α, and γ as well as the parametrization of β as a function of the zenith angle θ have been optimized with simulations.

Systematic Uncertainties
The main sources of systematic uncertainties for AMIGA are the uncertainty of the area-dependent module efficiencies, the uncertainty of the muon density estimate ρ 450 resulting from the unknown true shape of the MLDF for individual events, the uncertainty in the calibration of the PMT discriminator threshold voltages, the uncertainty due to density variations of the soil covering the scintillator modules, and the correction of the zenith angle-dependent attenuation of the muon density. The impact of each source of uncertainty is summarized in Table 1.

Module efficiency correction
Based on laboratory measurements, the efficiency of the 5 m 2 and 10 m 2 scintillator modules to identify individual   muons has been analyzed as a function of the applied inhibition window size in the reconstruction procedure as shown in Fig. 3. As a result of the increased light attenuation in the wavelength-shifting (WLS) fibers of double length, the efficiency of the 10 m 2 modules is reduced compared to the 5 m 2 modules. The estimated muon density for each module is therefore corrected by the areadependent efficiencies 5m 2 win7 = 1.04 and 10m 2 win7 = 0.95 that were obtained for a window size of seven bins according to ρ corr µ = ρ µ/ . The systematic uncertainty of ρ 450 resulting from the efficiency correction has been estimated with data as σ sys,eff/ρ 450 = 9.9% by reconstructing ρ 450 for each event both for the efficiencies derived for a window size of seven bins and for an infinitely large window ( 5m 2 inf = 0.94 and 10m 2 inf = 0.87).

Lateral distribution function
Due to fluctuations of the true function shape for individual events, a further systematic uncertainty arises from us- ing a parametrized slope β(θ) in the MLDF of Eq. (3). We estimated the induced systematic uncertainty of ρ 450 with simulations by reconstructing the same events with slopes β varied by ±σ β (θ) with σ β (θ) = 0.15 · β(θ). For illustration, the mean relative uncertainty σ sys ρ 450/ρ450, obtained with the QGSJetII-04 hadronic interaction model is shown as a function of θ in Fig. 4. On average, a relative systematic uncertainty of σ sys,MLDF/ρ 450 = 8.8% is obtained.

Calibration
The calibration of the AMIGA muon detectors equipped with PMTs consists in setting the discriminator threshold V Thr of each PMT channel to 30% of the mean channelspecific single photo-electron amplitude V SPE . Measurements with the AMIGA prototype array have shown that the spread of the discriminator thresholds after calibration is around 10% of the target values. We estimated the effect of a 2σ variation of V Thr , corresponding to thresholds of (1 ± 0.02)V Thr = (0.3 ± 0.06)V SPE , with simulations for different values of V Thr ranging from 20% to 40% of V SPE . A systematic uncertainty of σ sys,calib/ρ 450 = 3.9% was obtained by a linear fit to the data as shown in Fig. 5.

Soil density
Variations in the soil density are a further source of systematic uncertainty due to the resulting fluctuations in the attenuation of the muonic shower component measured by the buried AMIGA detectors. The soil density in the area of the Pierre Auger Observatory has been measured at three different positions at depths of 1 m, 2 m, and 3 m. We studied the impact of variations of the mean soil density ρ soil = 2.38 g/cm 3 by 3σ ρ soil with σ ρ soil = 0.05 g/cm 3 with simulations. Averaged over all considered angles a systematic uncertainty of σ sys,soil/ρ 450 = 2.8% has been found as shown in Fig. 6.

CIC correction
As a result of the longer path in the atmosphere and the increased amount of soil covering the detectors, the muonic component of inclined air showers gets attenuated. We parametrized the attenuation function f att (θ) = 1 + ax + bx 2 (4) with with data from the AMIGA prototype array using the constant intensity cut (CIC) method [11]. The determined parametrization, shown in Fig. 7, is used to define a zenith angle-independent estimator of the muon density at the optimal distance of r opt = 450 m. Treating the parameter uncertainties of f att as systematic uncertainty of the correction, we obtained a mean systematic uncertainty σ sys, f att/fatt = 2.3% for the zenith an-

First Results
We analyzed one year of data recorded by the AMIGA prototype array starting from October 2015. During this period the detectors were calibrated and running stably.
The same SD quality cuts as for previous official reconstructions were applied [12]. Furthermore, to ensure a good sampling of the shower by the AMIGA detectors, the SD station with the largest signal was required to lie within the AMIGA prototype array (see Fig. 1). To avoid large attenuation effects and statistical uncertainties due to the reduced detection areas, the zenith angle range was restricted to θ ≤ 45 • . A total number of 1057 events passed the required quality cuts.
For each event, the attenuation-corrected muon density ρ 35 was reconstructed. First, the reconstructed number of muons measured by the individual 5 m 2 and the 10 m 2 scintillator modules was corrected by the respective efficiencies of ε 5 = 1.04 and ε 10 = 0.95 (see Fig. 3). Subsequently, the muon density ρ 450 at the optimal distance of r opt = 450 m was obtained by a fit of the MLDF of Eq. (3). This density was then corrected for the effect of the zenith-dependent attenuation by calculating ρ 35 according to Eq. (5).
The reconstructed muon densities are shown as a function of the energy, reconstructed by the SD, in Fig. 8. The evolution of ρ 35 has been fitted as a power law with a maximum likelihood approach [13,14] taking into account both the threshold effect that is caused by the application of an energy cut to data, and the uncertainties in the muon density and energy estimates. The best fit solution is a = (1.75 ± 0.05(stat) ± 0.05(sys)) m −2 for the average muon density at 10 18 eV and  for the logarithmic gain. The evolution of the muon content in data is compared to simulations of proton and iron primaries in Fig. 9. To soften the strong energy dependence, the muon densities have been normalized by the energy. The slopes of b = 0.91 for iron and b = 0.92 for proton simulations which are obtained for both hadronic interaction models are slightly steeper compared to data. More strikingly, simulations fail to reproduce the observed muon densities which are between 8% (EPOS-LHC) and 14% (QGSJetII-04) larger than for simulated iron showers at an energy of 10 18 eV.
We quantified the discrepancy of the muon content in simulations and data by combining the AMIGA muon density measurements with independent measurements of the mean depth of shower maximum of the fluorescence detector at fixed energies of 10 17.5 eV and 10 18 eV. Using their linear dependences ln ρ 35 A = ln ρ 35 p + f ρ ln A (10) on the mean logarithmic mass ln A , the mean logarithmic muon densities ln ρ 35 in simulations have been related to the mean depths of shower maximum X max based on proton and iron simulations for both hadronic interaction models and primary energies. The comparison of the muon densities obtained by AMIGA with simulations at the mean depth of shower maximum measured by the FD for the same energy is X max /g cm   shown in Fig. 10. The discrepancy of the muon content in data and simulations is quantified in Table 2. To match the AMIGA measurements, the muon content in simulations would have to be increased between 38% for EPOS-LHC and 50% (10 17.5 eV) to 53% (10 18 eV), for the QGSJetII-04 hadronic interaction model.

Summary
In this work, we described recent improvements of the AMIGA muon reconstruction procedure, estimated the main systematic uncertainties and presented first physics results from the AMIGA prototype array. The largest systematic uncertainties result from the area-dependent efficiency correction of the muon densities recorded by individual detector modules (9.9%) and the lack of knowledge of the shape of the muon lateral distribution function for single events (8.8%). Considering further uncertainties due to the calibration procedure (3.9%), variations in the soil density (2.8%), and the CIC correction (2.3%), an overall uncertainty of 14.3% is obtained.
Analyzing one year of data recorded by the AMIGA prototype array, we found that current hadronic interaction models fail to reproduce the measurements in the considered energy range between 3 × 10 17 eV and 2 × 10 18 eV.
We quantified the disagreement of the muon content measured by AMIGA and simulations with independent FD measurements of the mean depth of shower maximum and found a muon deficit in simulations between 38% (EPOS-LHC) and 50 − 53% (QGSJetII-04) at fixed energies of 10 17.5 eV and 10 18 eV. The reported results will be discussed in detail in a forthcoming paper.