Muon Measurements with IceTop

We present the measurement of the density of GeV muons in near-vertical air showers by the IceTop array at the South Pole. The muon density is measured at 600 m and 800 m lateral distance from the shower axis in air showers between 1 PeV and 100 PeV. This result can be used to constrain hadronic interaction models by comparing it with the outcome of Monte Carlo simulations. We show that some models do not produce muon densities in agreement with this result unless an unphysical composition of the primary cosmic ray flux is assumed.


Introduction
One of the challenges in understanding cosmic-ray induced extended air showers is the accurate description of hadronic interactions over several decades in the centerof-mass energy, from a few tens of GeV to more than 10 6 GeV. The relevant interactions are in the forward fragmentation region, with most of the energy beamed into pseudo-rapidity ranges that are difficult to study in colliders because they lie very close to the incident beam. Their cross-sections cannot be computed from quantum chromodynamics, because the strong coupling for these interactions is too large for a perturbative solution. Instead, they are calculated using phenomenological models tuned to a variety of data sets from collider and fixedtarget experiments. Several hadronic interaction models are commonly used. The most recent versions are colloquially called post-LHC, since they take high-energy data from the LHC into account. These are SIBYLL 2.3 [1,2], QGSJet II.04 [3], and EPOS-LHC [4]. SIBYLL 2.1 [5] is a pre-LHC model that is still popular for simulations of air showers initiated by PeV-cosmic rays.
Cosmic-ray air showers can be understood within the Heitler-Matthews model [6]. After each interaction, a fraction of the energy is carried by neutral pions which promptly decay into gamma-rays, producing electromagnetic cascades. The remainder of the energy is carried by hadrons, mostly charged pions, that go on to further interact as part of a hadronic cascade. This process is repeated until the energy per particle is less than a critical energy, at which point the mesons decay, producing muons that can reach the ground. In this simple picture, the number of muons scales sub-linearly with the energy of a primary proton, (1) * e-mail: javierg@udel.edu * * http://icecube.wisc.edu/collaboration/authors/current which is confirmed by detailed air shower simulations. If an air shower initiated by a nucleus with A nucleons is approximated by a superposition of proton showers, with an energy E/A each, we find for the number of muons According to this toy model, an iron primary with 56 nucleons produces about 50 % more muons than a proton shower with the same energy, while in detailed simulations it depends on a complex interplay of many factors, and varies greatly depending on the hadronic interaction model used.
Measurements of muons in air showers, made by the Pierre Auger collaboration, have provided evidence for a discrepancy between data and simulations for primary energies above 1 EeV [7,8]. They show a deficit of muons in simulations. The HiRes-MIA collaboration previously reported a discrepancy between the simulated and measured number of muons in air showers [9]. The HiRes-MIA experiment measured air showers simultaneously with fluorescence telescopes and a ground array of muon detectors. The density of muons at 600 m from the shower axis was larger than the density in simulated air shower produced by iron. At the same time, the composition inferred from the depth of shower maximum showed a composition much lighter than pure iron. Taken together, these measurements point to a deficit of muons in simulations.
Other results are mixed. While the measurement of muons at Yakutsk (E > 2 × 10 10 GeV) appears to support the hypothesis of a deficit of muons in simulations [10], measurements conducted with the EAS-MSU surface detector at energies between 100 and 500 PeV do not yield a discrepancy [11]. Comparisons of measurements of the muon flux at energies between 10 3 and 10 4 GeV with simulations show a deficit or excess, depending on the model [12]. To further complicate the matter, the KASCADE-Grande collaboration has measured the atten- uation of the muon number with zenith angle using a constant-intensity cut analysis and found a significantly larger attenuation compared to several recent hadronic interaction models [13]. Earlier measurements of the muon component with the Akeno extensive air-shower experiment [14] did not show a discrepancy.
The characteristics of the IceCube Neutrino Observatory make it a unique instrument to probe the muon content of air showers. IceTop, its surface component, has been used to measure the energy spectrum of cosmic rays in the energy range between 1.6 PeV to 1.3 EeV [15]. The deep portion of the detector has been used to study the flux of atmospheric muons [16] and the lateral distribution of high energy muons in air showers [17]. Hybrid measurements, which used both detector components, have been used to infer mass composition of the primary cosmic ray flux [18]. The latter measurements were based on the fact that TeV muons, as opposed to the electromagnetic (EM) component of the air showers, are able to penetrate the ice shield above the deep part of IceCube.

The Detector
IceTop is an air shower array of 81 stations deployed in a triangular grid with a typical separation of 125 m. It was completed in 2008. The detector covers an area of roughly one square kilometer and is located above the IceCube detector at the geographic South Pole. Each station consists of two Cherenkov detectors separated by ten meters. The Cherenkov detectors are cylindrical tanks with an internal radius of 0.91 m. The detectors contain two Digital Optical Modules (DOMs) and are filled with clear ice up to 0.9 m. The 40 cm between the ice surface and the lid is filled with expanded perlite (amorphous volcanic glass, expanded to low density with grain sizes of the order of 1 mm) for thermal insulation and light protection [19]. Each DOM combines a 10 inch photo-multiplier tube (PMT) with electronics for signal processing and readout [20,21].
A local trigger occurs when the voltage in one of the DOMs in a tank exceeds a predefined discriminator threshold. Stations also have two trigger levels. A Hard Local Coincidence (HLC) occurs when both tanks in a station have a local trigger within a time window of 1 µs. If there is a local trigger in only one tank, it is called a Soft Local Coincidence (SLC). The total charge collected at the PMT's anode, after digitization and baseline subtraction, constitutes the tank's signal. In the case of HLCs, this is done offline, after a better estimate of the baseline becomes available, while in the case of SLCs this is done by the DOM firmware, using the best estimate of the baseline at the time. The calibrated charge is expressed in units of Vertical Equivalent Muon (VEM), which is the average charge produced by a vertically through-going muon.
The shower direction, the intersection point of the shower axis with IceTop (the shower core), and the shower size are estimated by fitting the measured charges with a lateral distribution function (LDF) and the times of the signals with a phenomenological model of the shower front [19,22]. The LDF and the shower front model are functions of the lateral distance, the distance of closest approach between the shower axis and the tanks. The LDF includes an attenuation factor due to the snow cover on top of each tank. To estimate the energy of the primary cosmic ray, we use the relationship between the shower size S 125 , defined as the signal at a lateral distance of 125 m, and the true primary energy. This relationship has been presented previously [15], and for cos θ > 0.95 it is the following: log 10 E = 0.938 log 10 S 125 + 6.018. ( It was derived assuming a specific composition model commonly referred to as H4a. This is a 4-component model based on the 5-component model by Gaisser [23], with the Nitrogen and Aluminum components merged into a single component with Oxygen as the representative element.

The Analysis Method
The procedure for this analysis has been described previously in this conference series [24][25][26]. It proceeds in two parts. The first part uses an empirical model of the response of the detector to estimate the muon density at different energies and lateral distances. In the second part, a multiplicative factor, derived using simulated data, is applied to obtain an unbiased estimate of the muon density.
Muons leave a characteristic signal that can be identified at large lateral distances, as shown in Figs. 1 and 2. roughly follows a power law over the entire lateral distance range. The other population, with signals around 1 VEM (log 10 (S /VEM) = 0), is visible only at large lateral distances (log 10 (R/m) 2.5 in the figure). This population is made up mostly of tanks hit by one muon. Measuring this last population allows us to estimate the number of muons. Figure 2 shows histograms of the signals at a fixed lateral distance, where one can see two populations clearly. A fit to these histograms yields the average number of muons at that energy, zenith angle and lateral distance.
Simulated data is used to derive a multiplicative factor that, after applying it, gives an unbiased estimate of the muon density. The resulting muon densities are then interpolated at two lateral distances, depending on the reconstructed primary energy: 600 m in the 1-40 PeV range, and 800 m in the 9-140 PeV.

Dataset
This analysis is based on data covering a 3-year period, totaling around 947 days of data acquisition. During this period, between May 31, 2010 and May 2, 2013 more than 18 million vertical events (sec θ < 1.05 or θ < 18 • ) were collected. The quality and selection cuts used in the analysis are a subset of the cuts used in previous publications [15]: • Events must trigger at least five stations and the reconstruction fit must succeed.
• Reconstructed core must be within the geometric boundary of the array.
• The station with the largest signal must not be at the edge of the detector.
• There must be at least one station with signal greater than 6 VEM.
After cuts, the energy resolution, defined as the standard deviation of the distribution of log 10 (E reco /E true ), is between 0.08 and 0.1 at 1 PeV and around 0.05 at 100 PeV.

Simulated Datasets
To estimate the cosmic ray air shower parameters, we used detailed simulations. The simulated air showers are produced with CORSIKA [27]. The CORSIKA showers are processed through a detailed simulation of the detector response. Each CORSIKA shower is resampled 100 times by uniformly distributing shower cores over an area larger than the detector. The maximum sampling radius is chosen as the largest distance where the shower can trigger the array. The simulation of the detector response includes a simulation of the entire hardware and data acquisition chain. The interactions with the IceTop tanks are simulated using the GEANT4 package [28]. The procedure accounts for interactions in the snow, small drifts in the calibration constants (on the order of 3%), and accidental coincidences. These are signals uncorrelated with the event that are coincident in time. Four primary types (H, He, O, Fe) were simulated. The energies of the simulated showers are distributed according to an E −1 differential spectrum between 0.1 and 300 PeV, and their zenith angles are uniformly distributed in sin 2 θ between 0 and 65 • . This accounts for the projection of the detector area on a plane perpendicular to the air shower direction.
Several datasets were combined in order to reproduce the experimental conditions. There is one collection of datasets for each campaign year covering the pe- The IC79.2010 collection of datasets, 60,000 showers per primary, is the same one used in a previous publication [15]. They were produced with CORSIKA v69900, using atmosphere 12, which is based on the July 1, 1997 South Pole atmosphere with an atmospheric overburden of 692.9 g/cm 2 (680 hPa). The IC86.2011 and IC86.2012 collections, 20,000 showers per primary each, were produced with CORSIKA v73700 using the average April 2011 atmospheric profile. In all these cases, the hadronic model used was Sibyll 2.1 [5] for interactions with energies greater than 80 GeV and FLUKA [29,30] at lower energies.
Smaller datasets were simulated using post-LHC hadronic interaction models QGSJet II.04 [3], EPOS-LHC [4] (CORSIKA v73700) and Sibyll 2.3 [1,2] (CORSIKA v75000). These datasets are used to determine the true muon densities at ground according to different hadronic models, and to verify that the systematic shifts expected from different hadronic models are within the uncertainty of the method.

Comparison to Simulations
The muon densities resulting from the first part of the analysis are compared directly with the densities of muons at  ground given by CORSIKA. The ratio between the simulated and reconstructed values is then used as a multiplicative factor. We have looked at the extreme cases of a cosmic ray flux composed entirely of protons (pure-proton) or iron nuclei (pure-iron), because the primary cosmic ray flux is expected to lie somewhere in between. The multiplicative factor that must be applied in the pure-proton case differs from the one for pure-iron. Therefore, an uncertainty in the composition of the primary cosmic ray flux produces a systematic uncertainty in the multiplicative factor. The multiplicative factor used in the second part of the analysis corresponds to a mix of equal parts proton/iron, with an uncertainty covering the entire range between proton and iron. The result of the analysis applied to air showers simulated with Sibyll 2.1, including the multiplicative factor, is displayed in Fig.3. The inverse of the multiplicative factor is shown in Fig. 4.
The different contributions to the systematic uncertainty are shown in Fig. 5. The statistical uncertainty in the multiplicative factor translates into a systematic uncertainty in the reconstructed muon density. After considering the systematic uncertainty of the factor itself, the multiplicative factor is the largest source of uncertainty in the analysis.
Several cross checks were performed in order to explore systematic effects. The effect of snow attenuation was observed by calculating the multiplicative factor on the three simulated yearly datasets separately. The yearly attenuation at 600 m is on the order of 8%, while the corresponding effect at 800 m was not significant when compared to the statistical fluctuations in the fit. All the other effects considered are of smaller magnitude. Their effect on the mean factor is smaller than the inter-year variation. A possible sampling bias was studied by performing the analysis on simulated events using their true energy and direction, and applying no event selection. The effect of the finite resolution in energy and direction reconstruction was explored by performing the analysis using various combinations of true/reconstructed quantities: true energy and geometry, reconstructed energy and geometry, true energy but reconstructed geometry, and vice versa. The procedure was also repeated with different compositions of the primary cosmic ray flux, resulting in no observable effect.
If one intends to use this technique to study the performance of hadronic models, then the effect of the choice  of hadronic interaction model used in the simulation is the most important cross-check. The result of this is shown in Fig. 6. Here we can see that the multiplicative factor for both EPOS-LHC and QGSJet II.04 are within the systematic errors quoted.

Results
The results of the analysis are the muon densities at 600 m and 800 m. In order to compare to the expectations from hadronic models, Fig. 7 shows the muon density divided by the expected muon density in proton showers. The brackets denote the systematic uncertainty from Fig. 5. The statistical uncertainty is too small to be seen at the lower end of the energy range. Note that the muon density for iron air showers in relation to proton showers is different for each model. Figure 8 shows the muon density transformed according to the following equation: where ρ p and ρ Fe are the muon densities in proton and iron showers respectively. It follows from equation 2 that the quantity in equation 4 is roughly proportional to the logarithm of the mean primary mass. In the same figure, we show the expected muon density according to three cosmic ray flux models. In the case of Sibyll 2.1, the muon density is consistent with the expectation from flux models. Perhaps it starts to deviate above 80 PeV, although this is not statistically significant. When interpreted according to EPOS-LHC, the measurement is well below the expectation for pure proton, which corresponds to an excess of muons in simulated showers. When compared to Sibyll 2.3 and QGSJet II.04 expectations, the measurement lies within the proton-iron range, although around 1 PeV it is lower than expectation from all cosmic ray flux models considered. This also points to an excess of muons in simulations. Around 100 PeV it is consistent with mixed composition models.

Summary
We have shown that IceTop can measure the muon densities in a model-independent way. We have looked at crosschecks of the method to understand the main sources of uncertainty, including the effect of snow, the detector finite resolution, selection effects, cosmic ray flux model, and hadronic interaction model. All of these effects are described by the quoted systematic uncertainties. The measured muon density was interpreted according to several hadronic interaction models using and compared to the expectations from cosmic ray flux models. At 1 PeV, post-LHC models produce less muons than observed, while at 100 PeV they are consistent with observations assuming a mixed composition model.