Latest Cosmic Ray Results from IceTop and IceCube

The IceCube Neutrino Observatory at the geographic South Pole, with its surface array IceTop, detects three different components of extensive air showers: the total signal at the surface, low energy muons on the periphery of the showers, and high energy muons in the deep InIce array of IceCube. These measurements enable determination of the energy spectrum and composition of cosmic rays from PeV to EeV energies, the anisotropy in the distribution of cosmic ray arrival directions, the muon density of cosmic ray air showers, and the PeV gamma-ray flux. Furthermore, IceTop can be used as a veto for the neutrino measurements. The latest results from these IceTop analyses will be presented along with future plans.


Introduction to Cosmic Ray Physics with the IceCube Neutrino Observatory
The IceCube Neutrino Observatory at the geographic South Pole, completed in December of 2010, is not only a world-class neutrino observatory but is also an excellent instrument to study cosmic rays. IceCube includes multiple detector components, as shown in Figure 1. The IceCube-InIce array consists of 86 strings buried beneath the surface of the Antarctic ice sheet to a depth of 2500 m. Below a depth of 1500 m, these strings are instrumented with 60 digital optical modules (DOMs) apiece [1]. The DOMs are designed to detect the Cherenkov light emitted by charged particles traversing the ice [2]. The strings are arranged in a triangular grid with ∼125 m separation, as shown in Figure 2. Each IceCube-InIce string is topped with two ice-Cherenkov tanks separated by 10 m. These two tanks are referred to as a station, and all the surface stations together comprise the IceTop array [3]. Each tank is viewed by two DOMs apiece, one running at low gain, the other at high gain, to maximize the dynamic range. Both the IceTop and IceCube-InIce DOMs are fully integrated into the data acquisition system of the observatory.

Observation Modes for Cosmic Ray Studies Using IceTop and IceCube-InIce
Since the IceTop and the IceCube-InIce arrays can be operated independently or in coincidence, there are a number of different cosmic ray studies that can be performed utilizing those three possible observation modes. We begin with a discussion of each observation mode: the analyses are discussed in the next sections. * e-mail: karen.andeen@marquette.edu Figure 1: The IceCube Neutrino Observatory, with multiple sub-arrays labeled. Here, analyses using IceTop and the IceCube-InIce arrays are discussed.

Studies using the IceTop Array alone
As discussed in [3], when six tanks in three stations register a signal in coincidence, the IceTop surface array is triggered and the signals from all tanks and the deep-ice detectors are preserved. IceTop data from each air shower event are then reconstructed using a maximum-likelihood procedure to fit the shape and normalization of a lateral distribution function (LDF) of the deposited charges. This reconstruction algorithm takes into account arrival time fluctuations and results in the fitted shower core position (x, y, z), direction (θ, φ), and (S 125 , β) . Here, β is related to the slope of the LDF, while S 125 is the "shower size" parameter, the result of the LDF fit to the signal strength measured in vertical equivalent muons (VEM) at a reference distance of 125 m perpendicular to the shower axis, as shown in Figure 3. At this distance, the shower size  16  17  18  19 20  21  22   23  24  25  26   27 28  29  30   31   32  33  34  35  36  37  38  39  40  41  42  43  44  45  46   47  48  49 50  51  52  53 54  55   56  57  58  59  60  61  62  63  64   65  66  67  68  69  70  71  72  73  74  75  76  77  has the lowest fluctuations and the smallest dependence on the mass of the primary particle, as shown in [4]. S 125 is directly related the energy of the primary cosmic ray, as shown in Figure 4. Snow accumulates at a different rate over each IceTop tank at an average rate of 20 cm per year, which decreases the expected electromagnetic signal. Thus, prior to fitting the lateral distribution function, the measured signals must be corrected. The snow is assumed to attenuate the electromagnetic signals exponentially with the slant depth through the snow to the detector. Therefore, the snow depth is measured biannually at each tank and interpolated across the detector throughout the year, and an absorption factor is applied which is optimized for each year individually, as detailed in [3,5].
The IceTop array is sensitive to air showers with energies ranging from ∼ 300 TeV to ∼ 2 EeV: the low-energy threshold is determined by the spacing of the stations 1 , and the high energy limit is determined by the total size of the array. After a standard selection for event quality, the uncertainty in the direction of events is ∼ 0.2 • at 30 PeV, and the energy resolution for protons at 30 PeV is ∼ 0.05 in log 10 (E/GeV) [3]. This very good energy resolution is partly due to the high altitude of the array (equivalent to a depth of ∼690 g/cm 2 ), which places the detector near the shower maximum in the relevant energy range. Both the all-particle energy spectrum, as discussed in Section 3.1, and the anisotropy of PeV cosmic rays, as discussed in Section 3.4, are studied using data from IceTop alone. The relationship between the signal, measured in vertical equivalent muons (VEM), and the radial distance from the fitted shower axis (in meters) is shown. The colored circles represent IceTop tanks that registered a signal: the color scale ranges from red (early hits) to blue (late hits) and the size of the circles is related to the amount of charge deposited in each tank. The signals are fitted to the LDF function, shown as a blue line, and the open circle is the fit result at 125 m, which is then the quantity S 125 for this event.  The relationship between log 10 (S 125 ) and the primary energy log 10 (E/GeV) for simulated proton events with cos(θ)>0.85. The color-scale on the z-axis denotes number of events. It is clear that S 125 is strongly correlated with the primary energy of the cosmic rays.
In addition to the studies listed above-which focus on the electromagnetic signals near the core of the air showers-IceTop's low trigger threshold of 1/6 VEM allows for the study of the GeV muon signal, which dominates far from the core of the shower. The GeV muon signal is studied as a check of the hadronic interaction models, as discussed in Section 3.3, and is also used to select muon-poor showers, which are candidate events for the PeV gamma-ray studies discussed in Section 3.5.

Studies using the IceCube-InIce Array alone
Although IceCube is designed as a neutrino observatory, signals in the IceCube-InIce Array are dominated by down-going charged cosmic ray events. These signals are deposited by muons that are produced early in the development of the air showers, with sufficient energy at production to reach the deep array (∼500 GeV minimum and typically ∼1 TeV). This energy threshold increases with increased overburden of ice, and therefore with increased zenith angle. Due to the large detector volume, the IceCube-InIce array alone collects ∼10 11 atmospheric muon events per year, which provides a large data-sample for cosmic ray anisotropy studies above 1 TeV, as discussed in Section 3.4.
Furthermore, the IceCube-InIce array provides additional information about the energy losses from high energy muon bundles. For each event, the pattern of hits in the ice can be translated into an energy loss profile as a function of slant depth [6]. The energy loss profile is then fit to provide three composition-sensitive parameters. First, a measurement of the fitted muon energy loss at X = 1500 m slant depth (dE µ /dX 1500 ) is a proxy for the total number of muons in the bundle, which is directly related to the mass of the primary cosmic ray. Next, two parameters are derived from the deviations from the fitted due to stochastic energy losses, which are also composition sensitive since iron-initiated bundles have more stochastics as the bundles contain more muons, while the energy losses from proton bundles can be more extreme. These parameters are used to make a measurement of the cosmic ray composition, as discussed in Section 2.3, while the composition analysis itself is detailed in Section 3.2.

Studies using the IceTop and InIce Arrays in Coincidence
Events with trajectories that pass through IceTop and IceCube-InIce can be reconstructed in both arrays on an event-by-event basis; thus, the combined reconstruction provides a handle on both the primary mass and energy of the cosmic rays simultaneously. The compositionsensitivity of the comparison of dE µ /dX 1500 from the  IceCube-InIce array with S 125 from IceTop for coincident events is illustrated in Figure 5 from simulations of protons and iron. The composition analysis is discussed in Section 3.2.
Additionally, events with a shower axis reconstructed by IceTop as intersecting the IceCube-InIce volume which then have no TeV muon tracks in the deep ice are likely to be initiated by PeV photons rather than hadrons, as shown in Figure 6 and discussed in Section 3.5.
Conversely, IceTop can be used as a veto: events with muons in the deep-ice detector whose track direction passes through the IceTop array, but which have no corresponding air shower in IceTop, are likely to be astrophysical neutrinos interacting in the ice between the two detector components, as shown in Figure 7 and discussed in Section 3.6.

Cosmic Ray Energy Spectrum using IceTop Alone
The all-particle energy spectrum using IceTop alone is derived from the measured S 125 spectrum using an unfolding technique which relies on the relationship between S 125 and primary energy for different groups of nuclei, as derived from Monte Carlo simulations. The Monte Carlo simulations are performed using the CORSIKA [7] air shower generator, with SYBILL 2.1 [8] and FLUKA [9] as high-energy and low-energy hadronic interaction models respectively. A detailed detector simulation is implemented in Geant4 [10], which models the tank response and the effect of snow on top of the tanks. These simulations are then used to determine the efficiency of IceTop as a function of primary energy. The composition assumption (i.e. the relative contribution of different mass groups vs energy) used to derive the unfolding is one of the main systematic uncertainties in extracting this energy spectrum. However, since the true energy spectrum should be independent of the zenith angle, the sensitivity of the spectrum measurement to the composition assumption is checked [5] by making the S 125 to primary energy unfolding independently for four bins in zenith angle under various composition assumptions [5]. The model which produces the most similar spectra at the four different zenith angles is presumed to be the closest to the actual cosmic ray composition. In this analysis, the H4a composition model [11] was the best model of those implemented. The simulated data are therefore weighted using the H4a model and then used to convert S 125 to primary energy. The angular dependence of the spectra reconstructed assuming H4a is used as a measure of the systematic uncertainty from composition.
This technique was first applied to one year of data (2010-11) from the partially completed IceTop detector (IT-73 with 73 of 81 stations in operation) [5]. The same analysis technique has been applied to three years of Ice-Top data (2010-2013). The data from the complete 81 station array is retriggered using only the IT-73 tanks for consistency across the three years (as shown in Figure 2). The three-year analysis includes a number of improvements, including a new treatment of the time-dependent correction for snow above the detector. Figure 8 shows the allparticle spectrum from the IceTop-alone analysis for each of the three years individually (colors) and all three years together (black). The systematic uncertainties are shown as the gray band. The three years of data agree within the statistical and systematic uncertainties. There is a clear hardening of the spectrum around 2 × 10 16 eV and a steepening above 2 × 10 17 eV.

Cosmic Ray Composition Using IceTop and IceCube-InIce in Coincidence
The same three-year data-set used for the IceTop alone energy spectrum described in Section 3.1 is also used for the coincident analysis of the energy spectrum and the By using a neural network technique, a massindependent energy spectrum and individual elemental spectra for primary groups are measured. This technique was first developed using one month of data with the 40string, 40-station partial detector configuration for primary energies up to 30 PeV [12], and was updated for the first year of the present data set in [13]. The analysis technique has now been further refined and applied to three years of data. The latest version of the network, as shown in Figure 9, is generated using the simulated air showers discussed in Section 3.1 for protons, helium, oxygen and iron primaries. Five input parameters are used to train the network: from IceTop, the energy proxy S 125 and the zenith angle; from IceCube-InIce, the muon number proxy dE µ /dX at 1500 m depth and two different selections to quantify the high-energy stochastic energy losses along the muon bundle track, as discussed in Section 2.2 and detailed in [13] and [14]. The two outputs of the neural network are the primary energy and a mass proxy.
The all-particle cosmic ray energy spectrum is derived on an event-by-event basis directly from the neural network energy output, taking into account the effective area and livetime of the detector. Figure 10 compares the total energy spectrum obtained by the IceTop-alone analysis (from Section 3.1) in blue and the coincidence analysis described here in black. The excellent agreement of the two independent analyses confirms that the composition systematic has been dealt with in a reasonable way in the IceTop only analysis, which has higher statistics due to the geometry of the coincidence requirement.
A measurement of the composition as a function of energy was performed on a statistical basis. The mass proxy distributions for each of the four elemental groups simulated (proton, helium, oxygen and iron) are used to generate probability "templates" for small intervals in reconstructed energy using an unbinned kernel density method [15]. These templates are then compared with the experimental data in order to determine the fractional contribution of each of the different mass groups in each interval of reconstructed energy. The resulting elemental energy spectra are shown in Figure 11, together with prediction from theoretical models H3a and H4a [11] and global ex-    : Individual spectra for the four mass groups (protons in red, helium in yellow, oxygen in green, and iron in blue) including total detector systematic compared with various predictions of cosmic ray composition (H3a, H4a, and GST [11] and GSF [16,17]) periment fits GST [11] and GSF [16,17]. The detector systematic uncertainties, which are dominated by the absolute calibration of the light yield in the detector, are shown in gray. The measured composition agrees with all predictions within the statistical and systematical detector uncertainties. The heavy components, represented by oxygen and iron, maintain a hard spectrum up to higher energies than proton and helium.

Low Energy (GeV) Muon Studies Using IceTop Alone
IceTop is sensitive to both the electromagnetic component of the air shower as well as the low-energy muonic part as discussed in Section 2.1 and detailed in [18] and [19]. While the electromagnetic component dominates near the shower axis, at distances far from the shower axis the probability that an individual tank is hit by a particle is much smaller than one. Thus, at large distances, when a single muon deposits its charge (defined as 1 VEM for verti- cal muons), the muon signal is statistically distinguishable from that produced by the electromagnetic particles. Thus the average muon density (ρ µ ) for air showers in a narrow span of primary energy is derived by counting the number of tanks with muon hits in a narrow band a certain distance from the shower core and accounting for the known cross-sectional area of an IceTop tank. Figure 12 shows the average lateral muon density profiles derived for near vertical showers. The muon density is then measured at two lateral distances which are far from the shower axis (in this case 600 m and 800 m) for each band in primary air shower energy. The densities expected from primary iron and protons from a variety of hadronic interaction models are then used to scale the experimentally measured densities between zero (for protons) and one (for iron) for each energy bin, as shown in Figure 13. These results indicate that while the pre-LHC Sibyll-2.1 hadronic interaction model had a deficit of muons at the highest energies, the new post-LHC models predict too many of muons at the lowest energies.

Anisotropy Studies using IceTop and IceCube-InIce
As discussed in Section 2.2, cosmic-ray induced air showers above TeV energies produce muons that are measured with the IceCube-InIce array. These muons are used to study the anisotropy in arrival directions of the primary cosmic rays. Similarly, as discussed in Section 2.1, IceTop is used to study the anisotropy in the PeV energy range. The latest study of this kind uses 6 years of data from the completed IceCube-InIce and IceTop arrays between May 2009 and May 2015, which resulted in a total of 318 billion cosmic ray induced muon events [20]. In this analysis the IceCube-InIce data are separated into nine independent energy bins with increasing median energy, and the IceTop-alone dataset provides one data  Figures  14 and 15, which show respectively the 13 TeV energy bin from IceCube-InIce, and the IceTop result at 1.6 PeV (which is consistent with the IceCube-InIce result at similar energies, but is derived from significantly more events). Furthermore, the dipole amplitude decreases with energy up to a few hundred TeV, and increases again at higher energies. The origin of this transition is not yet fully understood, but it may indicate a change in dominant sources of cosmic rays at those energies. Therefore combined analyses of the anisotropy with the chemical composition and energy spectrum of cosmic rays is planned for the future.
Additionally, a new analysis combines data from the IceCube and HAWC detectors in order to provide a nearly  [21] and only data that are consistent between the two detectors were selected. The relative intensity as a function of equatorial coordinates is computed by binning the sky into an equal-area grid with a bin size of 0.9 • using the HEALPix library [22]. The expected flux is calculated from the data themselves in order to properly account for rate variations in both time and viewing angle (this is not possible with simulations). The likelihood-based reconstruction developed in [23] is applied, which simultaneously fits cosmic ray anisotropies and detector acceptance. This likelihood method provides an optimal anisotropy reconstruction and the recovery of the dipole anisotropy for ground-based cosmic ray observatories located in the middle latitudes such as HAWC. A smoothing procedure is then applied. (Note the much higher number of events available in 5 PRELIMINARY Figure 16: Relative intensity of cosmic-rays at 10 TeV median energy (the "large scale anisotropy") PRELIMINARY Figure 17: Relative intensity after subtracting the multipole fit from the large-scale map above (the "Small scale anisotropy") years of IceCube data as compared to 1 year of HAWC data.) The combined dataset provides near full coverage of the sky.
The resulting maps show an anisotropy in the arrival direction distribution of 10 TeV cosmic rays that extends across both hemispheres, as shown in Figure 16. In order to eliminate larger structures, fitted multipoles are subtracted. Thus, in addition to a large-scale structure, significant small-scale structure is observed that is largely consistent with previous individual measurements, as shown in Figure 17.

PeV Gamma-Ray Studies using IceTop and IceCube-InIce
The IceCube-InIce array acting in coincidence with IceTop is the only detector in the Southern Hemisphere sensitive to PeV gamma-rays. As mentioned in Section 2.3, since extensive air showers initiated by gamma-rays are characterized by their low muon content and a shower maximum at deep atmospheric depths, combined information from the IceCube-InIce array and IceTop enable discrimination between air showers initiated by PeV gamma-rays,  First, an unbiased scan of the entire observable sky was developed in order to search for point sources of PeV gamma-rays [24]. No evidence of significant gammaray emission above the background expectation was observed; this analysis resulted in a declination-dependent 90% confidence level upper limit of (E/TeV) 2 dΦ/dE ∼ 10 −14 − 10 −13 cm −2 s −1 TeV −1 on source strength, the most stringent for PeV gamma-rays yet reported [24].
Next, the data were searched for correlations with both the 15 H.E.S.S. sources and the 11 neutrino events from the IceCube 4-year high-energy starting event (HESE) sample observable in the field of view of this analysis [24].
Again no significant gamma-ray emission above the background expectation was observed; however, PeV gammaray upper limits have been placed on the H.E.S.S. sources for the first time, thus constraining gamma-ray production scenarios.
Finally, a search was made for diffuse gamma-ray emission from the Galactic plane [25]: this analysis resulted in a 90% confidence level upper limit on the angular-integrated diffuse gamma-ray flux from the Galactic plane at 2 PeV under the assumption of an E −3 spectrum, as shown in Figure 18. This limit is an order of magnitude better than previous IceCube results, and is similar to a model prediction obtained by extrapolating the diffuse flux measured by Fermi-LAT into PeV energies. This limit is placed on a previously unobserved region of the Galactic plane, as illustrated in Figure 19.
IceCube's sensitivity to the gamma-ray flux is expected to improve at a rate lower than the usually expected inverse square root of the livetime due to additional exposure. This is due to the reduced acceptance to gamma-ray air showers with continued snow accumulation on IceTop tanks.

Using IceTop as a Veto for down-going Neutrino Events
A preliminary study has also been performed to determine the capability of IceTop to act as a veto for the IceCube neutrino analyses [26]. The two main backgrounds to the astrophysical neutrino searches in IceCube are the downgoing muons generated by cosmic rays incident on the atmosphere, and neutrinos generated through interactions in the atmosphere. Therefore, the IceCube extragalactic neutrino searches currently use the outer layer of IceCube-InIce DOMs to veto the muon tracks from the background interactions, resulting in a high purity astrophysical neutrino sample. However, using the outer shell as a veto greatly reduces the effective volume of the detector in these searches. If the IceTop surface array could be used instead to reject the background events arriving from the direction of IceTop (the coincident zenith range is <30 • ), the effective volume for astrophysical neutrino searches would be significantly increased to include any neutrino interactions that occur in the ice in the 1.5 km above the array, thus also extending IceCube's field of view toward the galactic center. In order to study the capability of IceTop to act as a veto for these astrophysical neutrino searches, an analysis was developed using 10% of the data collected in 2012. High quality down-going events were selected from data which have a muon track trajectory extrapolated to the surface well within the footprint of the IceTop array. Astrophysical neutrinos, to act as signal events, are then simulated and passed through the same quality and directional selection criteria. An IceTop trigger is not required, since the goal of the veto analysis is to search for neutrinos, which will not trigger the surface array (as shown in Figure  7); thus, three IceTop observables were developed specifically for this analysis. These observables are based on the air shower traversing the array which is expected if the Figure 20: The three log-likelihood ratios for one energy and zenith bin. The experimental data, which is dominated by the background sample, is shown in red, while the simulated astrophysical neutrino signal is shown in blue. (Note that in the LLHRs, δ represents the distance from the expected shower axis, τ represents the residual time from the expected shower front, and ρ represents the charge at the expected time.) IceCube-InIce signal is not produced by an astrophysical neutrino. These parameters are a lateral distance to the expected shower axis, a residual time in the tanks which is calculated with respect to the expected shower front, and the recorded charge at the expected time. Log-likelihood ratios (LLHR) are then defined from these parameters as shown in Figure 20, which are used to determine cuts that preserve a predefined fraction of the overall astrophysical neutrino signal (in this case 100% and 80%). These cuts are applied to the experimental data. The passing rate following these selections is shown in Figure 21 and the veto efficiency is shown in Figure 22. Where no events remain after the selections, the 68% confidence level upper limit is shown with arrows.
The veto analysis is ongoing with improvements to the log-likelihood discrimination method and a higher statistics sample (which will improve the efficiency at high energy). The method will be used to explore requirements for a future extended veto for IceCube-Gen2 as discussed in Section 4.

Future Plans
A new collaboration, IceCube-Gen2, was founded in 2015 to pursue the design and construction of the secondgeneration IceCube facility at the South Pole [27], which is envisioned to include a ∼10 km 3 high-energy deep-ice array to study the astrophysical neutrino flux. A surface air shower array with an area of ∼10 km 2 corresponding Figure 21: Event rates for cosmic ray and atmospheric neutrino background (data) and astrophysical neutrino signal (calculated with simulation), versus an IceCube-InIce muon energy proxy after applying a LLHR selection which retains 100% (blue) and 80% (green) of the astrophysical neutrinos. Arrows indicate upper limits. to the footprint of the deep detector, as well additional surface detectors beyond the footprint of the deep detector for use as a veto, are also planned. The increased zenith-angle range for coincident events will boost the coincident data rate by a factor of ∼50. The planned surface array will also enable lateral and production depth muon measurements for every point; thus, IceCube-Gen2 will be able to provide the most precise studies of the transition region from galactic to extra-galactic cosmic rays. Studies to optimize the surface component of IceCube-Gen2 for both veto and cosmic-ray physics are ongoing.