Cosmic Ray Spectrum and Composition from PeV to EeV Using 3 Years of Data From IceTop and IceCube

TL;DR

This study uses three years of IceTop and IceCube data to measure the cosmic ray energy spectrum and composition from PeV to EeV energies, employing neural networks and multiple analysis methods.

Contribution

It introduces a combined analysis approach using surface and in-ice detector data with neural networks to simultaneously determine cosmic ray spectrum and composition.

Findings

All-particle spectrum measured across PeV to EeV range.

Elemental group spectra reconstructed with systematic uncertainties analyzed.

Comparison of two analysis methods highlights their respective advantages.

Abstract

We report on measurements of the all-particle cosmic ray energy spectrum and composition in the PeV to EeV energy range using three years of data from the IceCube Neutrino Observatory. The IceTop detector measures cosmic ray induced air showers on the surface of the ice, from which the energy spectrum of cosmic rays is determined by making additional assumptions about the mass composition. A separate measurement is performed when IceTop data are analyzed in coincidence with the high-energy muon energy loss information from the deep in-ice IceCube detector. In this measurement, both the spectrum and the mass composition of the primary cosmic rays are simultaneously reconstructed using a neural network trained on observables from both detectors. The performance and relative advantages of these two distinct analyses are discussed, including the systematic uncertainties and the dependence on the hadronic interaction models, and both all-particle spectra as well as individual spectra for elemental groups are presented.

Tables

Table 1. Table 1: Fit parameters for converting IceTop shower size S 125 subscript 𝑆 125 S_{125} to energy in Eqn.  2 , using the “H4a” composition assumption. Errors are on the order of 0.006 for p 0 subscript 𝑝 0 p_{0} , and 0.0035 for p 1 subscript 𝑝 1 p_{1} .

Zenith range	$p_0$	$p_1$
	6.011	0.933
	6.055	0.924
	6.110	0.915
	6.177	0.907

Table 2. Table 2: Sytematic light yield shift

Effect	Light yield shift
+10% scattering	+3.6%
+10% absorption	$-$11.8%
$-$7.1% scattering and absorption	+7%
30 cm hole ice scattering	+4.5%
100 cm hole ice scattering	$-$2.9%
DOM efficiency	$\pm 3\%$
Total Light Yield Effect	+9.6%,$-$12.5%

Table 3. Table 3: Total flux IceTop-only analysis.

Energy	Bin Width	Flux	$\pm$ Stat.	$+$ Det.Syst. $-$
(log10(Ereco /GeV))		GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1
6.450	0.10	1.730	0.001	0.131	0.151	$\times {10}^{-13}$
6.525	0.05	1.063	0.001	0.085	0.097	$\times {10}^{-13}$
6.575	0.05	7.667	0.006	0.612	0.734	$\times {10}^{-14}$
6.625	0.05	5.508	0.005	0.448	0.512	$\times {10}^{-14}$
6.675	0.05	3.927	0.004	0.330	0.380	$\times {10}^{-14}$
6.725	0.05	2.787	0.003	0.235	0.264	$\times {10}^{-14}$
6.775	0.05	1.967	0.002	0.168	0.190	$\times {10}^{-14}$
6.825	0.05	1.390	0.002	0.120	0.139	$\times {10}^{-14}$
6.875	0.05	9.742	0.015	0.859	0.943	$\times {10}^{-15}$
6.925	0.05	6.797	0.012	0.614	0.694	$\times {10}^{-15}$
6.975	0.05	4.745	0.009	0.433	0.455	$\times {10}^{-15}$
7.025	0.05	3.288	0.007	0.305	0.338	$\times {10}^{-15}$
7.075	0.05	2.289	0.006	0.212	0.225	$\times {10}^{-15}$
7.125	0.05	1.604	0.004	0.142	0.158	$\times {10}^{-15}$
7.175	0.05	1.122	0.004	0.104	0.108	$\times {10}^{-15}$
7.225	0.05	7.893	0.028	0.709	0.709	$\times {10}^{-16}$
7.275	0.05	5.591	0.022	0.488	0.547	$\times {10}^{-16}$
7.325	0.05	3.960	0.018	0.368	0.385	$\times {10}^{-16}$
7.375	0.05	2.828	0.014	0.244	0.248	$\times {10}^{-16}$
7.425	0.05	2.011	0.011	0.188	0.192	$\times {10}^{-16}$
7.475	0.05	1.430	0.009	0.128	0.132	$\times {10}^{-16}$
7.525	0.05	1.031	0.007	0.093	0.091	$\times {10}^{-16}$
7.575	0.05	7.431	0.057	0.599	0.631	$\times {10}^{-17}$
7.625	0.05	5.350	0.045	0.507	0.486	$\times {10}^{-17}$
7.675	0.05	3.754	0.036	0.317	0.423	$\times {10}^{-17}$
7.725	0.05	2.640	0.028	0.258	0.245	$\times {10}^{-17}$
7.775	0.05	1.943	0.023	0.168	0.141	$\times {10}^{-17}$
7.825	0.05	1.369	0.018	0.109	0.140	$\times {10}^{-17}$
7.875	0.05	1.008	0.015	0.121	0.086	$\times {10}^{-17}$
7.925	0.05	6.912	0.115	0.529	0.704	$\times {10}^{-18}$
7.975	0.05	4.948	0.092	0.441	0.393	$\times {10}^{-18}$
8.050	0.10	2.934	0.046	0.291	0.273	$\times {10}^{-18}$
8.150	0.10	1.468	0.029	0.142	0.179	$\times {10}^{-18}$
8.250	0.10	6.991	0.178	0.725	0.551	$\times {10}^{-19}$
8.350	0.10	3.146	0.106	0.282	0.324	$\times {10}^{-19}$
8.450	0.10	1.602	0.068	0.144	0.155	$\times {10}^{-19}$
8.550	0.10	6.468	0.382	0.637	0.798	$\times {10}^{-20}$
8.650	0.10	3.169	0.239	0.390	0.280	$\times {10}^{-20}$
8.750	0.10	1.350	0.139	0.443	0.355	$\times {10}^{-20}$
8.850	0.10	6.091	0.837	0.612	0.842	$\times {10}^{-21}$
8.950	0.10	2.492	0.480	0.514	0.440	$\times {10}^{-21}$
9.050	0.10	8.199	2.472	4.278	2.228	$\times {10}^{-22}$

Table 4. Table 4: Total flux coincidence analysis.

Energy	Bin Width	Flux	$\pm$ Stat.	$-$ Det.Syst. $+$
(log10(Ereco /GeV))		GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1
6.55	0.1	9.005	0.008	0.587	0.668	$\times {10}^{-14}$
6.65	0.1	4.612	0.005	0.317	0.377	$\times {10}^{-14}$
6.75	0.1	2.323	0.003	0.168	0.203	$\times {10}^{-14}$
6.85	0.1	1.158	0.002	0.087	0.103	$\times {10}^{-14}$
6.95	0.1	5.635	0.014	0.446	0.524	$\times {10}^{-15}$
7.05	0.1	2.698	0.008	0.220	0.254	$\times {10}^{-15}$
7.15	0.1	1.308	0.005	0.094	0.117	$\times {10}^{-15}$
7.25	0.1	6.473	0.032	0.473	0.592	$\times {10}^{-16}$
7.35	0.1	3.250	0.020	0.235	0.266	$\times {10}^{-16}$
7.45	0.1	1.657	0.013	0.116	0.133	$\times {10}^{-16}$
7.55	0.1	8.600	0.083	0.596	0.685	$\times {10}^{-17}$
7.65	0.1	4.399	0.053	0.377	0.430	$\times {10}^{-17}$
7.75	0.1	2.194	0.033	0.137	0.219	$\times {10}^{-17}$
7.85	0.1	1.132	0.021	0.098	0.086	$\times {10}^{-17}$
7.95	0.1	5.532	0.133	0.483	0.574	$\times {10}^{-18}$
8.10	0.2	2.062	0.048	0.148	0.205	$\times {10}^{-18}$
8.30	0.2	4.619	0.181	0.429	0.549	$\times {10}^{-19}$
8.50	0.2	9.966	0.666	1.019	1.495	$\times {10}^{-20}$
8.70	0.2	2.237	0.250	0.124	0.206	$\times {10}^{-20}$
8.90	0.2	5.283	0.964	0.720	1.021	$\times {10}^{-21}$

Table 5. Table 5: Proton and helium group flux.

Energy	Bin Width	Flux	$-$ Stat. $+$		$-$ Det.Syst. $+$
(log10(Ereco /GeV))		GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1		GeV-1 m-2 s-1 sr-1
Proton
6.55	0.1	2.791	0.448	0.471	1.253	1.289	$\times {10}^{-14}$
6.65	0.1	1.350	0.170	0.173	0.637	0.629	$\times {10}^{-14}$
6.75	0.1	7.838	1.055	1.092	3.834	3.790	$\times {10}^{-15}$
6.85	0.1	3.456	0.542	0.563	2.013	2.045	$\times {10}^{-15}$
6.95	0.1	1.393	0.201	0.207	0.785	0.816	$\times {10}^{-15}$
7.05	0.1	6.684	1.462	1.550	4.106	4.335	$\times {10}^{-16}$
7.15	0.1	2.656	0.601	0.635	1.771	1.943	$\times {10}^{-16}$
7.25	0.1	1.040	0.230	0.231	0.643	0.727	$\times {10}^{-16}$
7.35	0.1	5.291	1.556	1.671	3.677	4.417	$\times {10}^{-17}$
7.45	0.1	3.172	0.807	0.851	2.218	2.168	$\times {10}^{-17}$
7.55	0.1	1.294	0.405	0.412	1.082	1.037	$\times {10}^{-17}$
7.65	0.1	8.205	1.609	1.022	4.838	4.461	$\times {10}^{-18}$
7.75	0.1	2.536	0.646	0.691	1.167	1.866	$\times {10}^{-18}$
7.85	0.1	1.194	0.379	0.431	0.796	1.246	$\times {10}^{-18}$
7.95	0.1	8.318	2.957	2.574	6.256	2.819	$\times {10}^{-19}$
8.10	0.2	2.422	1.118	0.703	1.441	1.101	$\times {10}^{-19}$
8.30	0.2	5.257	1.781	1.589	3.663	1.860	$\times {10}^{-20}$
8.50	0.2	1.340	0.437	0.406	0.631	0.845	$\times {10}^{-20}$
8.70	0.2	3.562	2.825	3.316	2.908	5.462	$\times {10}^{-21}$
8.90	0.2	0.016	0.016	2.319	0.000	2.115	$\times {10}^{-21}$
Helium
6.55	0.1	2.232	0.981	0.974	0.761	0.514	$\times {10}^{-14}$
6.65	0.1	1.520	0.401	0.412	0.373	0.367	$\times {10}^{-14}$
6.75	0.1	4.798	2.109	2.168	2.146	1.984	$\times {10}^{-15}$
6.85	0.1	2.801	1.166	1.161	1.360	1.174	$\times {10}^{-15}$
6.95	0.1	1.682	0.406	0.412	0.301	0.219	$\times {10}^{-15}$
7.05	0.1	6.067	3.120	3.049	2.295	1.549	$\times {10}^{-16}$
7.15	0.1	4.146	1.407	1.399	0.804	0.592	$\times {10}^{-16}$
7.25	0.1	2.241	0.649	0.674	0.290	0.108	$\times {10}^{-16}$
7.35	0.1	7.893	3.603	3.526	2.571	1.885	$\times {10}^{-17}$
7.45	0.1	2.194	1.693	1.668	0.738	1.121	$\times {10}^{-17}$
7.55	0.1	1.482	0.708	0.720	0.521	0.718	$\times {10}^{-17}$
7.65	0.1	0.118	0.118	3.143	0.118	2.046	$\times {10}^{-18}$
7.75	0.1	2.316	1.380	1.379	0.948	0.449	$\times {10}^{-18}$
7.85	0.1	1.402	0.922	0.858	1.067	0.110	$\times {10}^{-18}$
7.95	0.1	3.182	3.182	6.307	0.000	8.377	$\times {10}^{-19}$
8.10	0.2	0.081	0.081	2.774	0.081	2.159	$\times {10}^{-19}$
8.30	0.2	0.000	0.000	3.264	0.000	7.476	$\times {10}^{-20}$
8.50	0.2	0.001	0.001	7.633	0.001	0.070	$\times {10}^{-21}$
8.70	0.2	4.096	4.096	4.094	4.096	2.171	$\times {10}^{-21}$
8.90	0.2	2.785	2.785	1.181	2.785	0.165	$\times {10}^{-21}$

Table 6. Table 6: Oxygen and iron group fluxes.

Energy	Bin Width	Flux	$-$ Stat. $+$		$-$ Det.Syst. $+$
(log10(Ereco /GeV))		GeV-1 m-2 s-1 sr-1	GeV-1 m-2 s-1 sr-1		GeV-1 m-2 s-1 sr-1
Oxygen
6.55	0.1	1.975	0.972	0.942	0.857	0.479	$\times {10}^{-14}$
6.65	0.1	6.714	5.305	5.097	5.909	2.597	$\times {10}^{-15}$
6.75	0.1	5.803	2.068	1.978	2.040	1.016	$\times {10}^{-15}$
6.85	0.1	2.766	1.115	1.108	1.137	0.690	$\times {10}^{-15}$
6.95	0.1	1.206	0.395	0.386	0.355	0.112	$\times {10}^{-15}$
7.05	0.1	8.757	3.037	3.055	2.852	1.583	$\times {10}^{-16}$
7.15	0.1	2.831	1.354	1.352	1.531	0.594	$\times {10}^{-16}$
7.25	0.1	9.002	7.871	7.627	9.002	6.196	$\times {10}^{-17}$
7.35	0.1	9.934	3.663	3.712	6.989	3.545	$\times {10}^{-17}$
7.45	0.1	6.723	1.528	1.538	3.490	1.441	$\times {10}^{-17}$
7.55	0.1	2.647	0.625	0.614	2.095	1.069	$\times {10}^{-17}$
7.65	0.1	2.200	0.311	0.203	1.149	0.340	$\times {10}^{-17}$
7.75	0.1	9.188	1.493	1.517	4.586	3.184	$\times {10}^{-18}$
7.85	0.1	3.866	0.849	0.889	2.925	2.045	$\times {10}^{-18}$
7.95	0.1	1.792	0.601	0.522	1.792	0.719	$\times {10}^{-18}$
8.10	0.2	7.187	2.751	1.536	7.187	4.681	$\times {10}^{-19}$
8.30	0.2	2.101	0.452	0.396	1.786	0.731	$\times {10}^{-19}$
8.50	0.2	1.055	1.055	1.058	1.055	3.217	$\times {10}^{-20}$
8.70	0.2	0.001	0.001	3.993	0.001	0.000	$\times {10}^{-21}$
8.90	0.2	0.001	0.001	1.561	0.001	2.915	$\times {10}^{-21}$
Iron
6.55	0.1	2.008	0.436	0.460	0.897	1.733	$\times {10}^{-14}$
6.65	0.1	1.071	0.257	0.274	0.475	0.945	$\times {10}^{-14}$
6.75	0.1	4.795	0.894	0.949	2.073	4.133	$\times {10}^{-15}$
6.85	0.1	2.561	0.477	0.489	1.107	2.061	$\times {10}^{-15}$
6.95	0.1	1.354	0.165	0.176	0.496	0.989	$\times {10}^{-15}$
7.05	0.1	5.468	1.395	1.428	2.907	5.366	$\times {10}^{-16}$
7.15	0.1	3.450	0.548	0.567	1.447	2.791	$\times {10}^{-16}$
7.25	0.1	2.292	0.332	0.352	0.953	1.529	$\times {10}^{-16}$
7.35	0.1	9.378	1.726	1.759	5.213	9.010	$\times {10}^{-17}$
7.45	0.1	4.482	0.679	0.703	2.183	4.418	$\times {10}^{-17}$
7.55	0.1	3.177	0.328	0.353	1.539	2.486	$\times {10}^{-17}$
7.65	0.1	1.367	0.145	0.162	0.630	1.353	$\times {10}^{-17}$
7.75	0.1	7.904	0.842	0.853	4.190	6.535	$\times {10}^{-18}$
7.85	0.1	4.862	0.441	0.445	2.142	3.552	$\times {10}^{-18}$
7.95	0.1	2.590	0.293	0.299	1.135	1.809	$\times {10}^{-18}$
8.10	0.2	1.093	0.121	0.135	0.480	0.683	$\times {10}^{-18}$
8.30	0.2	1.993	0.334	0.354	0.815	1.818	$\times {10}^{-19}$
8.50	0.2	7.570	1.043	1.092	3.116	2.111	$\times {10}^{-20}$
8.70	0.2	1.471	0.251	0.286	0.133	0.419	$\times {10}^{-20}$
8.90	0.2	2.481	1.036	1.266	1.043	1.693	$\times {10}^{-21}$

Table 7. Table 7: Correlation coefficients from KDE fits.

Equations

$$S=S_{125}\cdot {\left(\frac{r}{125m}\right)}^{-\beta -\kappa {\log }_{10}\left(\frac{r}{125m}\right)},$$

$$S_{red}=S\cdot e^{\frac{-d}{\lambda \cos \theta }},$$

$$\frac{\mathrm{d}E_{\mu}}{\mathrm{d}X}(X)>a\cdot\bigg{(}\frac{\mathrm{d}E_{\mu}}{\mathrm{d}X}(X)\bigg{)}_{\mathrm{reco}}^{b},$$

$${\log }_{10}(E/\mathrm{GeV})=p_0+p_1{\log }_{10}(S_{125}/\mathrm{VEM}).$$

Full text

IceCube Collaboration

Measurements of the cosmic ray energy spectrum and mass composition in the PeV to EeV energy range provide key information about the origin and propagation of cosmic rays in what is commonly considered to be a transition region from galactic to extragalactic cosmic rays. In the present paper we report on composition and energy spectrum measurements using three years of data from the IceCube Neutrino Observatory.

The IceCube Neutrino Observatory is a versatile particle detector located at the geographic South Pole, with both surface and deeply-buried components. The latter, or InIce, component (described in detail in Aartsen et al. (2017)) consists of 5160 Digital Optical Modules (DOMs) Abbasi et al. (2008a) deployed with 17 m spacing on 86 strings in a 125 m triangular grid formation, at depths from 1450 m to 2450 m below the surface in transparent ice. Each DOM contains a 10 inch Hamamatsu photomultiplier tube (PMT) and electronics for signal processing and readout Abbasi et al. (2009).

The surface component, IceTop (described in detail in Abbasi et al. (2013a)), is an array of pairs of tanks filled with water that has frozen and containing two DOMs each, operating at different PMT gains for increased dynamic range; the pairs of tanks are called stations, and each station is located above a string of the InIce detector (see Figure 1).

As Cherenkov detectors, both components are sensitive to the charged particles in the extensive air showers (EAS) produced by cosmic rays from the Southern Hemisphere sky, and both are used to extract the cosmic ray energy spectrum in different ways. In this paper, we will discuss two methods for extracting this information: the first uses IceTop only and will be referred to as the IceTop-alone analysis, while the other uses both the IceTop and InIce detectors in tandem, and will be referred to as the coincident analysis. The coincident analysis also provides a new measurement of the composition of cosmic rays.

The IceTop detector can be used alone to measure the core position, direction, and size of the air showers at the surface. These observables are utilized to extract a cosmic ray energy spectrum: in particular the shower size, discussed in Section I.4, is strongly correlated with the energy of the incident cosmic ray primary. It is important to note that as the shower size depends slightly on the mass of the cosmic ray primary, the composition must be assumed in order to extract the energy spectrum when using IceTop-alone. Since IceTop collects a large number of showers, rare high-energy events are collected in sufficient numbers to extend the analysis into the EeV range. This analysis was performed most recently in Aartsen et al. (2013a), in which one year of 73-station data was used to measure a spectrum from a few PeV to EeV.

The InIce detector measures the energy loss of the high-energy muons in the deep ice, which is strongly dependent upon the mass of the incident cosmic ray primary, as discussed in Section I.5. Thus, with the IceTop and InIce detectors working in tandem, the high-energy muon component of the air showers are measured in coincidence with the electromagnetic component; thus both the energy spectrum and mass composition are measured without making assumptions about one to determine the other. However, this method requires coincident events between the two components of IceCube: due to the long lever arm between the two arrays, only zenith angles of approximately 0-$30^{\circ}$ on the sky pass the coincident selection criteria, which then yield fewer events from the same data sample as the IceTop-alone analysis. Thus, the coincident analysis cannot reach as high an energy as IceTop-alone. A coincident analysis like this was performed in Abbasi et al. (2013b), in which one month of data from the half-completed 40-station, 40-string detector was used to measure a spectrum and average logarithmic mass ($\langle\ln A\rangle$) from 1 to 30 PeV; that analysis was extended with improved reconstruction techniques in Feusels et al. (2013); Feusels (2014) using one year of data from the nearly complete 73-station, 79-string array, achieving better resolution and reaching to 1 EeV.

In this paper, we will present for the first time the coincident analysis in detail. In addition, the one-year analyses of IceTop-alone Aartsen et al. (2013a) and coincident Feusels et al. (2013); Feusels (2014) data are extended to three years, and we report these improved and updated results.

I Data, Simulation, and Reconstruction

I.1 The 3-year Data Set

The analyses described here use three years of data, from June 1, 2010, through May 2, 2013. The IceTop-alone and the coincident analyses use the same dataset, and thus have both a total livetime of 977.6 days (with a negligible uncertainty of less than half a percent). The IceCube Neutrino Observatory was not running in its complete configuration (81 stations, 86 strings) until May 14, 2011. Therefore, the first year of data included here was taken in the incomplete 73-station, 79-string configuration (IT-73/IC-79). The simulation used in these analyses was also produced using the IT-73/IC-79 detector configuration. Thus, in order to handle the two following years of data in the same analysis, all information from the final 7 deployed strings and 8 stations has been removed from data processing. The configuration used is shown in Figure 1.

I.2 Simulation

Both of the analyses presented here use the same set of Monte Carlo simulations of cosmic ray events to establish relationships between detector observables and cosmic ray energy and mass. (This simulated dataset was also used in Aartsen et al. (2013a).) 30000 air showers of four primary types (protons, helium, oxygen, and iron) were simulated using the CORSIKA air shower generator Heck et al. (1998) with an $E^{-1}$ spectrum between $\log_{10}(E/\mathrm{GeV})$ = 5.0 and 8.0. Additionally, thinned CORSIKA showers Heck and Pierog ; Hillas (1997) were generated at higher energies: 12000 showers of each type, again with an $E^{-1}$ spectrum between $\log_{10}(E/\mathrm{GeV})$ = 7.0 and 9.5. High-energy muons are not thinned in the algorithm. The range of energy overlap between $\log_{10}(E/\mathrm{GeV})$ = 7.0 and 8.0 allows for verification of the un-thinning algorithm Abbasi et al. (2013a).

Two hadronic interaction models were used: FLUKA Battistoni et al. (2007) below 80 GeV, and Sibyll 2.1 Ahn et al. (2009) above 80 GeV. EGS4 Nelson et al. (1985) was used to model the electromagnetic interactions. Other high energy hadronic interaction models are used for systematic studies, as discussed in Section IV.3. The zenith angle from the primary particles is sampled from a $\cos(\theta)\sin(\theta)$ distribution between 0∘ and 40∘, while the azimuth angle is drawn from a uniform distribution over the whole 2$\pi$ azimuth range. To be able to study atmospheric changes and apply corrections where needed, a reference atmosphere was chosen based on the MSIS-90-E parametrization Hedin (1991) of the South Pole Atmosphere on July 1, 1997, which has a ground pressure of 692.2 g/cm2 at the South Pole altitude (CORSIKA atmosphere 12, Heck and Pierog ). To make more efficient use of the CORSIKA showers available, each shower is copied, or resampled, 100 times, and thrown at random locations within a circle of radius $R$ centered on the center of the IceTop array. The resampling radius is the largest possible for the shower to trigger the array Abbasi et al. (2013a).

Next, the particles generated by CORSIKA are propagated into the detectors. The individual responses of the IceTop tanks are simulated using a detailed Geant4 Agostinelli et al. (2003); Allison et al. (2006) model which takes into account the individual IceTop tank properties, including snow and air above the tanks, and the detector electronics. The resulting signal is converted into units of photoelectrons using constants unique to each tank Abbasi et al. (2013a). In this way the same calibration procedure (described in Section I.3) can be performed on both simulated events and experimental data.

For the coincident analysis, the high-energy muons in the CORSIKA air showers must also be propagated through the Antarctic ice and through the deep InIce detector. Muons with energy above 273 GeV111273 GeV is the energy at which 0.1% of the muons are expected to reach the top of the InIce detector and could create a detectable amount of Cherenkov light within the array Feusels (2014) are propagated Chirkin and Rhode (2004) through the Antarctic ice to the bottom of the InIce array. Propagating the Cherenkov photons from the muons through the South Pole ice to the DOMs by directly tracking each one is computationally prohibitive. Therefore, light profiles for GeV emitters (including both pure Cherenkov emission and more diffuse emission from cascade light sources) are tabulated in a software package Lundberg (2007), which includes a model of the full structure of the ice properties Aartsen et al. (2013b). The expected number of photoelectrons and their arrival times at each DOM are retrieved from the tables for each muon and for electromagnetic and hadronic cascades. Simulated noise hits are then added. Finally, the simulated photoelectrons are fed into a simulation of the InIce readout electronics and the detector trigger. The simulated DOM signals then follow the same processing chain as the experimentally measured DOM signals.

I.3 Pulse cleaning and Calibration

Only events which pass the cosmic ray filters Abbasi et al. (2013a) and which contain at least 6 hard local coincidence (HLC222HLC hits occur when both tanks in one station are hit within 1 $\mu$s) DOMs within 6 $\mu$s are processed further. In a first signal cleaning of selected events, signals from all DOMs which were determined to be unreliable at the time of data taking are removed. The remaining signals are calibrated using the procedure described in Aartsen et al. (2017), which returns the number of photoelectrons. At this point the IceTop and InIce data are split.

The IceTop data is cleaned and calibrated as described in detail in Abbasi et al. (2013a). At this stage, only HLC hits are preserved. All signals within one large trigger window are then split into clusters of hits that are likely related to one air shower. Next, an additional noise cleaning procedure is applied to the calibrated HLC signals. This procedure begins with a cluster of 3 hit stations and adds more tanks to the event using a simple distance-per-time requirement Feusels (2014). This extra cleaning procedure improves the removal of noise hits, of tanks with time-fluctuations, and of muon-like signals at large lateral distances from the shower core. These, in turn, improve the IceTop reconstruction procedures described in Section I.4. The specific tank response is then taken into account using a calibration procedure which transforms the signal into that expected from vertical (equivalent) muons (VEMs).

For the coincident analysis, the InIce hits causally connected to the event seen in IceTop are selected based on an allowed time difference window between IceTop and InIce hits. Therefore, before the IceTop and InIce events can be connected, InIce noise hits that affect the trigger time must first be removed in a procedure based on distance-per-time requirements (similar to that in IceTop). Then, InIce triggers matching the IceTop triggers and pulses are selected, and the pulses not connected to the selected trigger are removed. A final noise removal procedure in the InIce detector uses the reconstruction of the track by IceTop whereby only hits connected in time to the track, and in a cylinder around it, are preserved. These steps effectively remove random coincidences between the two detectors, and events that are very close in time.

Finally, after the above event cleaning and calibration procedure, events are required to have hits in at least 5 IceTop (IT-73) stations and, in the coincident analysis, hits in at least 8 InIce (IC-79) DOMs.

I.4 IceTop Reconstruction

Cleaned data from IceTop tanks are processed by a reconstruction software package called Laputop, which has been described in detail in Abbasi et al. (2013a). For each event, Laputop finds the best-fit shower core position ($x_{c}$, $y_{c}$, $z_{c}$) and direction ($\theta$, $\phi$), as well as two parameters describing the shape of the lateral distribution function (LDF) of deposited charge ($S_{125}$, $\beta$). The functional form of the charge LDF is a double logarithmic parabola:

[TABLE]

where $\beta$ is a measure of the steepness of the LDF and $S_{125}$ is the signal expectation at a perpendicular reference distance of 125 meters from the shower axis Klepser (2008). $S_{125}$ will be referred to throughout this paper as the shower size. $\kappa$ scales with the curvature of the parabola which is approximately constant for all hadronic showers; thus $\kappa$ is presently set as a default to 0.303, while the other two parameters are allowed to vary event by event. The best-fit parameters are found using a 3-step maximum-likelihood technique, which compares the timing and charge of the hits to both the expected charge LDF and an expected timing LDF. Both saturation of the tanks, as well as stations which are not hit, are taken into account in the likelihood.

The Laputop reconstruction also takes the actual snow depths on top of the tanks into account. The frozen water tanks of IceTop were deployed flush with the surface of the snow at the site. However, wind-blown snow continuously drifts over the array and covers the tanks with an overburden that increases over time and varies from tank to tank. Figure 2 shows the depths of snow covering the array’s tanks in each of the three years analyzed in this work: 2010, 2011 and 2012. The site accumulates about 20 cm of snow per year on average.

The total signal $S$ observed by the two DOMs in the IceTop tank consists of snow attenuated electromagnetic ($e^{\pm}/\gamma$) particles and unattenuated muons. If the attenuation of the snow is not taken into account, the reduction in total signal amplitude will make an event look less energetic, or “smaller”, than it really is. To take snow attenuation into account, a simple exponential reduction Sred is applied to the expected $S$:

[TABLE]

which only depends on the slant depth of snow overburden for the tank ($d/\cos\theta$), and an effective attenuation length $\lambda$, which takes into account absorption/generation in the shower and particle-type specific absorption behavior. The same attenuation length $\lambda$ is applied throughout the array, in the same way for all showers of any size. However, because of the increasing snow load from year-to-year, the optimal $\lambda$ also changes from year-to-year. Thus each of the three years of data was optimized separately to find the $\lambda$ which best creates agreement in the $S_{125}$ spectrum across different regions in the array (deeply-buried, and sparsely-buried). These best-fit values of $\lambda$ are: 2.1 meters for 2010/11, 2.25 meters for 2011/12, and 2.25 meters for 2012/13. Furthermore, the snow depth, $d$, has two sources of uncertainty De Ridder (2019). First, sastrugi333Sastrugi are irregular waves formed on the surface of the snow by wind erosion. cause variations in the snow depth across a single tank. Sastrugi heights measured in-situ were observed to follow a Gaussian distribution with a standard deviation of 4 cm. Second, the depth at each IceTop tank is measured twice per year, in February and in November, with an occasional third measurement in January; for all other times, a linear interpolation between these measurements is used to estimate the daily snow depth. Finer measurements of accumulation at the South Pole (made monthly by the Antarctic Meteorological Research Center at a site near IceTop) exhibit variations around a smooth interpolation with a $\sigma$ of 27% of the difference between the November and February measurements ($\Delta H_{snow}$, in centimeters). The two sources of uncertainty are combined and the snow depth $d$ in Eqn. 1 (the signal reduction due to the snow) is smeared by a Gaussian with a sigma of $\sqrt{(4\,\text{cm})^{2}+(0.27\cdot\Delta H_{snow})^{2}}$ De Ridder (2019).

I.5 InIce Reconstruction

The IceCube detector observes the Cherenkov light pattern from the high-energy muon bundles that propagate through the Antarctic ice as well as their accompanying energy losses (which create Cherenkov-light-emitting cascades).

These measured observables provide a handle on the primary composition: iron-induced showers are more muon-rich than proton-induced showers of the same energy, due to the superposition model of nucleons (detailed in e.g. Kampert and Unger (2012)). Therefore, iron-induced showers are more muon-rich than proton-induced showers of the same energy, therefore iron-induced showers will result in a greater overall deposit of Cherenkov light in IceCube than proton-induced showers of the same primary energy. Furthermore, proton-induced showers are more likely than iron-induced showers to have extremely high-energy muons in the bundle, which can create larger local energy depositions from Bremsstrahlung. Consequently, proton-induced showers are expected to create fewer but higher-energy stochastic losses in the detector than iron-induced showers (of the same primary energy). On the other hand, since iron-induced showers have more total muons than those induced by protons (of the same primary energy), an iron-induced shower is likely to undergo more lower-energy stochastic losses in the detector than a proton-induced shower (of the same primary energy). Therefore, counting these stochastic fluctuations (henceforth “stochastics”) gives us additional composition-sensitive information.

In order to measure these composition-sensitive parameters for coincident events, the track position and direction as reconstructed by IceTop (with Laputop) is used, and an energy loss reconstruction algorithm (Millipede) uses the timing and charge information of the observed Cherenkov light in the ice to create a profile of energy losses along the track as a function of slant depth, with 20-meter segmentation (as discussed in detail in Aartsen et al. (2014a)).

Figure 3 shows an example energy loss profile from a cosmic ray event. Note that sections of the energy loss profile corresponding to IceCube’s dust layer444The “dust layer” is a distinctive thick layer of dust deposited several millennia ago and currently located $\sim$1950-2050 m beneath the surface of the ice Aartsen et al. (2013b). (where there are consequently few photons) and near the boundaries of IceCube’s volume (where reported energy losses can be irregular) are removed (as shown in gray in the figure). The gaps seen in the energy loss profile are not segments on the muon bundle track where no energy was lost, but rather are pieces of the track that are not well-sampled by the detector.

The energy loss profile is then fit to extract two composition-sensitive parameters: a) the average energy loss behavior, which is indicated as the red line in Figure 3, and b) the size and quantity of deviations from that average behavior (the stochastics). The energy loss observable ($dE_{\mu}/dX$) is defined as the value of the fit to the energy loss profile at a fixed slant depth of X=1500 m, which corresponds roughly to the top of the IceCube detector (marked on the left side of Figure 3).

Two methods of selecting a number of high-energy stochastics from a energy loss profile are used in this work: a standard selection (marked as the red dashed line in Figure 3 and a strong selection requiring higher stochastic energy loss (marked as the red dotted line in Figure 3 Feusels (2014); De Ridder (2019). The selection criterion for the stochastics is given by:

[TABLE]

where $a=5$ and $b=0.8$ for the standard selection, and $a=7$ and $b=0.9$ for the strong selection (with appropriate dimensions for $a$ and $b$).

Figure 4 demonstrates the composition-sensitivity of the InIce observables reconstructed by Millipede. The energy loss parameter is directly comparable to the number of high-energy muons in the air shower and is therefore the primary composition-sensitive observable, as shown in Figure 4, left. The stochastics provide additional composition-sensitivity, as shown in Figure 4, center and right: Iron bundles have a larger number of stochastic losses since they have more muons, but the energy losses from proton bundles can be more extreme since the same total energy is transferred to fewer muons. The number of stochastics with the standard selection can separate masses best at low muon multiplicities, below 100 PeV. From about 30 PeV and above, for bundles containing at least 100 muons, the stronger selection performs better. In both cases there are more high energy cascades selected for iron bundles than for proton.

The muon multiplicity in air showers detected by IceCube shows a seasonal variation, which is due to the semi-annual alternation between the polar day and night, and the accompanying temperature changes in the atmosphere. The measured variation of $\log_{10}(\mathrm{d}E_{\mu}/\mathrm{d}X)$ is found to be 10-15% of the proton-iron difference. Since simulations are only performed with one atmosphere, the July,1997 atmosphere, all other months of data need to be corrected with respect to July. This correction correlates the changing temperature profile of the entire atmosphere, weighted with the muon production depth profile, with the measured variation of $\log_{10}(\mathrm{d}E_{\mu}/\mathrm{d}X)$. This energy dependent correlation factor is used as a correction. A small but symmetric variation of $\pm 3\%$ in the proton iron-space remains, which smears the data slightly. For more details, see De Ridder et al. (2013).

I.6 Quality Cuts

Quality cuts are used to ensure a sample of events which will be well reconstructed. The cuts for the IceTop-alone analysis are described in Aartsen et al. (2013a) while the cuts for the coincident analysis are detailed here.

IceTop Selections: Many of the IceTop selections revolve around the success of the Laputop reconstruction algorithm, which does an excellent job of reconstructing contained events (those with a shower core inside the area of the array) but its performance suffers for uncontained events.

• •

The number of stations after cleaning is required to be $\geq 5$.

• •

The largest snow-corrected charge measured in any tank is required to be at least 6 VEM.

• •

The station with the highest deposited charge is not allowed to be at the edge of the detector.

• •

The neighboring tank in the same station as the tank with the largest signal must have at least 4 VEM.

• •

The fraction of hit stations within a circle centered on the center of gravity of the shower with outer radius at the furthest hit station must be greater than 0.2.

• •

The reconstruction algorithm Laputop is required to converge.

• •

The LDF slope parameter $\beta$ is required to be between 1.4 and 9.5.

• •

The core location of the air shower must be reconstructed within a scaled factor of 0.96 of the area of the array.

InIce Selections: The InIce reconstruction begins with a fixed track position and direction from Laputop. Therefore the InIce quality selections focus on ensuring an accurate reconstruction of the energy loss from Millipede.

• •

The track position and direction calculated by Laputop is required to pass within the In-Ice instrumented volume.

• •

A minimum of 8 InIce DOMs are required to be hit.

• •

The Millipede energy loss reconstruction must succeed with log${}_{10}(rlogl)<2.0$ and the total charge predicted (QTOT) must be at least 90% of that measured ( log${}_{10}(\frac{QTOT_{predicted}}{QTOT_{measured}})>-0.03$ ) .

• •

At least 3 reconstructed cascades remain after all previous selections and after removal of cascades in the dust layer and at the edge of the detector (as discussed above).

A note about the zenith angle: As in Aartsen et al. (2013a), the IceTop-alone analysis presented here is divided into four bins of zenith angle, the steepest of which is limited to $cos(\theta)\geq 0.80$. No explicit cut in zenith angle is applied for the coincident analysis, since the solid angle acceptance of the two detectors together limits the zenith angle range to $cos(\theta)\sim 0.85$.

I.7 Performance

The reconstruction procedure and quality cuts described above yield a set of events with a core position resolution of 6-20 meters, and a track direction resolution of 0.3-1.0 degrees. These values depend on energy: Figure 6 shows the position and angular resolutions as a function of energy, for the coincident sample and the IceTop-alone sample (divided into high-zenith and low-zenith angles). Note that the detector reaches 100% efficiency for all particle types at $\sim$3 PeV (Feusels (2014)). As the size of the air showers approach and exceed the size of the IceTop array and the tanks become saturated, the reconstruction becomes less precise.

Figure 6 shows the reconstructed energy loss from IceCube compared to the reconstructed $S_{125}$  parameter from IceTop: in this parameter space there is clearly a strong separation between proton and iron primaries.

II The IceTop-alone analysis

The IceTop-alone analysis is sensitive to the energy spectrum of cosmic rays up to the EeV energy range. The reconstructed shower size parameter ($S_{125}$) in particular is highly correlated to the energy of the primary air shower. Figure 7 shows the relationship between $\log_{10}(S_{125})$ to $\log_{10}(E/\mathrm{GeV})$. In Aartsen et al. (2013a), a function relating $\log_{10}(S_{125})$ to $\log_{10}(E/\mathrm{GeV})$ was derived using Monte Carlo simulations divided into four ranges of zenith angles. For a given zenith range, the distributions of true energy for each slice of 0.05 in $\log_{10}(S_{125})$ were weighted using the H4a model from Gaisser (2012) as a composition assumption, and fitted with a Gaussian. (Since neither silicon nor magnesium were simulated, but both are included in the H4a model, simulated oxygen was weighted by the sum of CNO and MgSi model components.) The fitted mean of the Gaussian was then used as the energy estimate for that slice in $\log_{10}(S_{125})$. The functional form of the conversion is:

[TABLE]

Using updated simulation and reconstruction algorithms, the mapping of $\log_{10}(S_{125})$ to $\log_{10}(E/\mathrm{GeV})$ from Aartsen et al. (2013a) has been re-optimized and applied to the 3-year data set. The updated fit parameters for Eqn. 2 are in the following table:

The energy bias and resolution of this technique are shown in Figure 8. The reduced precision beyond 8.0 (100 PeV) is related to the reduced angular and position resolution shown in Figure 6, which creates an extra smearing effect in $S_{125}$.

Figure 9 shows the result of the IceTop-alone 3-year analysis for each year individually and combined, multiplied by a factor of $E^{3}$ to highlight the details. The gray band represents the total systematic uncertaintity of the IceTop detector, as described in Aartsen et al. (2013a). These results are consistent between the three years and are also consistent with those previously published from one year of data Aartsen et al. (2013a).

III The Coincident analysis

When the surface observables from IceTop are combined with the additional observables from the InIce detector, the high-energy muon component of the shower is measured in coincidence with the electromagnetic component of the EAS. Using this coincident configuration, a mass-independent primary energy spectrum and individual elemental spectra are measured. This technique was developed for the measurement of the cosmic ray composition of one month of IT-40/IC-40 data in the energy range between 1 PeV to 30 PeV Andeen (2011); Abbasi et al. (2013b) using two input variables. Building on this experience, the technique was extended to five input variables over a wider primary energy range, optimized over a larger scan of different network types, and trained on more Monte Carlo simulated events. This updated technique was applied to a single year of data from the nearly complete IT-73/IC-79 detector in Feusels et al. (2013); Feusels (2014). Here, the one year analysis is improved and further expanded to include three years of data.

III.1 Neural Network Mapping Technique

This analysis includes five variables which depend on primary energy and primary mass in a non-linear fashion: the shower size in IceTop ($S_{125}$), the zenith angle ($\cos(\theta)$), the muon energy loss in the ice ($\mathrm{d}E/\mathrm{d}X$), and the number of high-energy stochastics under two selections (standard and strong). There is no theoretical analytical expression that relates our input variables to primary mass and primary energy; thus, an artificial neural network (NN)555In particular, a feed-forward multilayer-perceptron (MLP) neural network is used from the TMVA Hoecker et al. (2007) machine learning package. is trained on simulation to determine the relationships between the five inputs and the two outputs. The network is strongly dependent on the two primary parameters, $S_{125}$ and $\mathrm{d}E/\mathrm{d}X$, but the three other parameters do contribute to the energy and mass reconstruction.

The final high-quality sample of simulated Monte Carlo data is split into three parts. Half of the sample is used to generate the neural network (the network sample). The other half (the verification sample) is used for comparisons of data and simulation in the final analysis steps. The network sample is again split in two: 74357 events are used to train the network (the training sample), the remaining 67399 events (the test sample) serve to test the network and to select the network architecture and optimal activation function based on the network performance. Networks were trained on unweighted events; however, every Monte Carlo sample mentioned above is chosen in such a way that it contains an equal mixture of each of the four primary types (p, He, O and Fe) and covers the full energy range.

During the first 5000 of 10000 minimizer iterations (also called cycles or epochs), only a random selection of 60% of the training data is utilized. After the training converged on this random selection, the training continues on the full training set.

III.2 Optimizing the Neural Network

Many different neural network architectures were evaluated for performance before analyzing any data, as discussed in Feusels (2014). In addition to networks with 5 inputs as described above, alternative networks with the 2 primary inputs ($\log_{10}(S_{125})$ and $\log_{10}(\mathrm{d}E/\mathrm{d}X)$ only), 3 inputs (adding $\cos(\theta)$) and 4 inputs (adding the standard selection of high energy stochastics only) were tested. Three groups of network structures were explored: with one, two, and three hidden layers, and the number of neurons was varied within the hidden layers. Two activation functions (a sigmoid, and a tanh) were explored. In total, 207 networks for each of the two activation functions and for each number of inputs (1656 networks in total) were trained on the simulations.

The performance of each network was assessed according to how well it reconstructed primary energy and primary mass. The assessment process was optimized to find the network with the smallest and most consistent RMS spread and bias over all energies, and which had mass groups that were best-separated and most distinctive (i.e. “peaky”). The final optimized network has 5 inputs, 7 neurons in a first hidden layer, 4 neurons in a second hidden layer, and 2 outputs, with a tanh activation function connecting the neurons and a linear mapping from the last layer to the output neurons. A schematic of this network is shown in Figure 10.

It is important to note that this neural network has two target outputs which are very different in nature: the first output is a continuous energy distribution, the second target output is instead is composed of four discrete numbers corresponding to four elemental masses simulated. Therefore, the neural network energy output (E0,reco) is also a continuous distribution which is expected to reproduce the true primary energy (within some bias and resolution) for each event, as discussed below in Sec. III.3. On the other hand, the neural network mass output results in smeared distributions around the four discrete mass numbers, which require further analysis in order to decompose the primary mass. The mass is therefore not reconstructed on an event-by-event basis but is determined statistically for the entire data set, as discussed below in Sec. III.4.

III.3 Neural Network Primary Energy Reconstruction

The energy dependence of the primary energy bias and resolution as reconstructed by the NN are shown in Figure 11. The energy resolution (Figure 11, lower) ranges from 9% (for iron showers at around 30 PeV) and 18% with the worst resolutions below the energy threshold of $\sim$3 PeV and at the highest energies due to the worsening core position and angular resolution (as discussed in Section I.7). Heavier primaries can be reconstructed more precisely because of their lower intrinsic shower fluctuations. As mentioned in Section II, the overall decrease in precision beyond $\sim$100 PeV is partially caused by the decrease in precision in angular and position resolution shown in Figure 6, which creates an extra smearing effect in $S_{125}$.

In this analysis, events are divided into energy bins of width 0.1 in $\log_{10}(E/\mathrm{GeV})$ , which is larger than both the energy bias and the energy resolution as shown in Figure 11. However, due to the decrease in accuracy, precision, and available statistics at high energies ($\log_{10}(E/\mathrm{GeV})$ > 8.0), bins of width 0.2 are used in this region. Above 1 EeV the energy bias dependence on the primary type becomes too large and limits the energy range over which this analysis is optimal.

Figure 12 shows the all-particle energy spectrum results for the coincident analysis for the three years individually and combined, multiplied by a factor of $E^{3}$ to highlight the details: the results are consistent between the years. The gray band represents the combined systematic uncertainties of the IceTop and InIce detectors for the coincident analysis, as discussed in Section IV.2. These results are included in Table 4 in Appendix .4.

III.4 Composition Reconstruction using kernel density estimation to fit neural network templates

Figure 13 shows histograms for each simulated element (proton, helium, oxygen and iron) in the natural logarithm of the neural network mass output for one slice in reconstructed energy. (The four simulated types (proton, helium, oxygen and iron) are equidistant in $\langle\ln A\rangle$, but not in $A$. Thus, the histograms are expected to be more distinct in logarithmic space.) In every slice in energy, the histogram for each primary element is converted into a template probability density function (p.d.f.) using an adaptive kernel density estimation (KDE) method Cranmer (2001). The template p.d.f.’s are shown as the solid lines in Figure 13. The template p.d.f.’s for all energy slices used in this analysis are given in Appendix .5 in Figure 26. The four primary types exhibit four very distinctive shapes in each slice in energy over the whole energy range. At $\log_{10}(E/\mathrm{GeV})$ = 8.0, the template p.d.f.’s begin to exhibit greater overlap due to the limited statistics in the Monte Carlo sets, reducing the composition sensitivity of the analysis. Beyond $\log_{10}(E/\mathrm{GeV})$ = 9.0 (1 EeV), the analysis becomes unreliable due to the overlap and reduction in data statistics.

Using the Roofit package Verkerke and Kirkby (2003), the set of four template p.d.f.’s were then weighted to find the fractions (which were constrained to add to unity) which best fit the NN mass output for the experimental data in the same slice in reconstructed energy666This serves a similar purpose to the chi-squared minimization approach described in Andeen (2011); Abbasi et al. (2013b); however, the new unbinned extended likelihood technique improves on the previous method by correctly taking into account the Poisson fluctuations of the bin contents in both the data and the templates, which is particularly relevant in bins with few events.. The result of this method applied to the experimental data for the same energy bin as in Figure 13 is shown in Figure 14 (upper). Additionally, the correlation between the fitted weights is shown in the form of uncertainty elipses in Figure 14 (lower). (The fit plots for all energy bins are given in Figure 27 in Appendix .5 and the corresponding correlation coefficients are shown in Table 7 in Appendix .5.) The resulting fractions of neighboring elements (i.e. protons and helium) are anti-correlated, while those from distant elements (i.e. protons and iron) are virtually uncorrelated: this means a proton primary is more likely to be confused for a helium primary than for an iron primary, which is expected.

The KDE template-fitting procedure yields a measurement of the fractions of each of the four nuclear mass groups (represented by H, He, O, and Fe), for each bin in energy. The fractions are shown in Figure 15. Each of these four individual fractions is then translated into an individual spectrum for the corresponding elemental group, as shown in Figure 16 (colors) compared to the all-particle spectrum (black). Recent model predictions are also included in Figures 15 and 16, which will be discussed in Section V. In both figures, the gray band represents the total coincident detector uncertainty from both the IceTop and InIce arrays, which will be discussed in Section IV.2. These results are included in Tables 5-6 in Appendix .4.

Intermediate elements, not part of the four groups listed above, are expected to produce neural network outputs in between the adjacent groups, so will partially contribute to the flux of the groups that bracket it. In order to test this, a small sample of silicon was passed through the NN + KDE chain and treated as “data”. The natural log of the mass of silicon is approximately midway between that of oxygen and that of iron; therefore, as expected (due to the regression-style neural network mass output), the silicon is reconstructed as a nearly 50/50 mixture of oxygen and iron across all energies.

Figure 17 shows the mean log mass, which is derived from the individual fractions shown in Figure 15. Again, the gray band represents the total coincident detector uncertainty from both the IceTop and InIce arrays, which will be discussed in Section IV.2. Each of the three years of data are again shown both separately and combined and agree very well within the statistical and systematic uncertainties.

IV Systematic Uncertainties

The uncertainties in the coincident analysis reported here can be grouped into three categories: analysis method, detector effects, and the hadronic interaction model.

IV.1 Analysis Method

As described above, the shape of the Monte Carlo templates is derived from an adaptive KDE method, which determines the optimal width of the Gaussian kernel function. To check the robustness of the composition fitting results, the optimal kernel width is artificially modified by a factor of $\pm 90\%$, resulting in either very jagged or very smooth templates. These artificially modified templates are then used in place of the optimal templates for the remainder of the analysis in order to measure the effect of the vastly different templates on the results. The variation in the final results due to the modified template shapes is so small that it is not visible in a figure; however, the values for uncertainty are included in the tables in Appendix .5.

IV.2 Detector Uncertainties

Three main detector effects contribute to the uncertainty in the composition results: the snow correction, the absolute energy scale of IceTop, and the light yield in the ice.

IV.2.1 Snow correction

Although the method for snow correction described in Section I.4 works well on average over the entire energy and zenith angle region, it is not perfect. This is predominantly due to its inability to distinguish between the electromagnetic and muonic component of the air shower. A systematic uncertainty of $\pm$0.2 m was assigned to $\lambda$, which covers the variations due to the unknown composition, energy and zenith angle dependence, as discussed in the Appendix of Aartsen et al. (2013a).

IV.2.2 Absolute energy scale

To obtain a reliable energy scale, the simulated response of each IceTop tank is calibrated using the signal from atmospheric muons Abbasi et al. (2013a). This procedure has proven to be very reliable, but due to the unknown composition, atmospheric conditions, etc, a final 3% uncertainty on the charge calibration, and thus on the absolute energy scale, needs to be taken into accountVan Overloop (2011). This translates directly to a 3% shift in S125.

IV.2.3 Light Yield

The uncertainty on the photon detection efficiency by the DOM (referred to as the DOM efficiency) is $\pm$3%, which implies a $\pm$3% possible variation of the total light yield observed. The uncertainty caused by the photon propagation in the ice includes scattering and absorption coefficients in the bulk ice and an effective scattering length of the hole ice, the ice in the drill hole around the DOM Rongen et al. (2016). For the bulk ice scattering and absorption coefficient, 3 points are taken on the 1 $\sigma$ error ellipse around the nominal values, shown in Aartsen et al. (2013c): a +10% scattering coefficient, a +10% absorption coefficient and a -7.1% scattering and absorption coefficient. Alternate effective scattering lengths of 30 cm and 100 cm were used for the hole ice model. For all these individual systematic uncertainties, simulations were produced and their effect was studied. It was found that the main combined effect is to influence the light yield in the DOM, and that this effect is rather independent of the initial light yield and zenith angle of the muon bundle. This means that all those systematic uncertainties can be combined and modeled as a shift on the observed light yield. Furthermore, the systematic uncertainties are (nearly) uncorrelated; thus the various errors and shift on light yield are added in quadrature, which gives a total light yield uncertainty of +9.6% and -12.5%. The individual contribution to the observed light yield shifts, as well as the total light yield uncertainty are given in Table 2.

IV.2.4 Total Detector Uncertainties

The three detector uncertainties discussed above (snow correction, absolute energy scale, and light yield) are added together in quadrature into a total detector uncertainty. Figures 18, 19, and 20 show individual and combined contributions to the uncertainty in the all-particle energy spectra, mean log mass, and the individual elemental spectra. In all figures, the gray bands are the total combined detector uncertainty which match the gray bands shown in Figures 12, 17, and 16, respectively, and are included in the tables in Appendix .4.

IV.3 Hadronic Interaction Model

The influence of the hadronic interaction model on the measurement of the cosmic ray composition arises mainly from variations between the models in the predicted number of high-energy muons. Sibyll 2.1 was used as the hadronic interaction model for the baseline simulation datasets, while the three post-LHC models–QGSJET II-04Ostapchenko (2011), Sibyll 2.3Riehn et al. (2016) and EPOS-LHCPierog et al. (2015)–are used as alternate models. The variation between the models in the high-energy ($>300$ GeV) muon number is smaller than for the surface (GeV) muons, but is still of order 15% at maximum (between QGSJet II-04 and EPOS-LHC). 777Here it is important to emphasize that the muon energy loss, $dE_{\mu}/dX$, as measured in IceCube is a proxy for the number of high energy muons in the muon bundle. These are muons in excess of 300 GeV which are created near the first interaction of the cosmic ray primary with the atmosphere. Recently, a number of studies (summarized here: Dembinski et al. (2019)) show a discrepancy between the experimentally measured number of muons (from many experiments), and the simulated number of muons from latest post-LHC hadronic interaction models. This discrepancy is in the number of lower energy muons, which are produced late in shower development, after many interactions. This discrepancy presently seems to be due to the accumulation of small discrepancies at each interaction; thus, the impact is large for the muons produced after many interactions, and is small/non-existent for those muons produced near the first interaction (Riehn et al. (2019)). Furthermore, this discrepancy has only been measured at cosmic ray primary energies above 1 EeV and it seems to increase with primary cosmic ray energy. ALICE (the only other experiment which has measured the high energy muons in the range between PeV and EeV) reports relatively good agreement between the number of muons in the bundles produced by PeV to EeV cosmic rays and QGSJET II-04 (Adam et al. (2016)), although the number of events is small. Thus, we presently expect little impact on the results here. However, this “muon puzzle”, is still under intense investigation by the community.

It is important to note that the main Sibyll 2.1 dataset included full samples across all energies for each particle type (proton, helium, oxygen, and iron); however, due to the computational time involved in generating a full simulated data set ($\sim$60000 CPU years), the alternate models include 1/10th the number of events, in proton and iron only. As a result of the limited size of the alternate simulated data sets, a full analysis is not repeated using those datasets for the simulated templates. Rather, the difference in the number of high-energy muons and in $S_{125}$ in the alternative sample with respect to Sibyll 2.1 is calculated for each particle type and applied to the experimental data, which is then passed through the full analysis chain. (Neither the stochastics nor the zenith angle are taken into account in this estimate due to their low impact on the results.) The results of this process are then weighted using the reconstructed elemental fractions for each particle type from the baseline Sibyll 2.1 analysis in order to obtain an estimate of what the final result would have been had the alternative interaction model been used as the nominal Monte Carlo simulation instead of Sibyll 2.1. It is probable that the alternate results reported here differ slightly from the actual fractions that would be reconstructed with a full alternative simulated dataset.

The estimated uncertainty due to each alternative model is shown in Figures 21, 22, and 23, which depict the effect of the choice of hadronic interaction model in the all-particle energy spectra, mean log mass, and the individual elemental spectra. Although the uncertainties due to the hadronic interaction models is large, the shapes of the distributions with respect to energy do not change significantly between the models.

V Results, Discussion and Outlook

V.1 Energy Spectra

Figure 24 compares the energy spectra resulting from the IceTop-alone analysis and the coincident analysis as described herein. The analyses are consistent with each other within the statistical and systematic uncertainties (only the smaller IceTop-alone systematic uncertainties are shown, for clarity). This good agreement indicates that the dependence of the IceTop-alone analysis on composition model has been effectively mitigated through the analysis technique, since the IceTop/InIce coincident analysis doesn’t require prior knowledge of the composition.

V.2 Mass composition

The elemental spectra results, shown in Figure 16, agree well with the recent H3a and H4a phenomenological models of the transition region between galactic and extragalactic cosmic rays (Gaisser (2012)), in which heavier elements retain a harder spectral index to higher energies. Within the statistical and systematical uncertainties the elemental spectra are also compatible with the phenomenological data fits, Gaisser-Stanev-Tilav (GST) Gaisser et al. (2013) fit and Global Spline Fit (GSF) Dembinski et al. (2018). These elemental spectra correlate with an increase in the mean-log-mass as a function of energy until about 100–200 PeV, as shown in Figure 17 (which is also derived from the individual fractions). Beyond this energy, given the statistical and systematic uncertainties the data are consistent with a composition that is either unchanging or decreasing. In Figure 25, our results are compared with those from other recent experiments: the results reported here indicate a higher flux in the iron group at high energies. As shown in Section IV.3, the absolute scale of the composition is strongly dependent on which hadronic interaction model is used for the simulations. In fact, the uncertainty due to the choice of hadronic interaction model is the biggest limitation on our analysis.

V.3 Outlook

The all-particle energy spectrum results presented here are consistent with each other and with previously published IceCube results Aartsen et al. (2013a). The limiting systematic effect is the uncertainty of the snow coverage over the tanks. The composition analysis results presented here are significantly improved from previously published results, which included only one month of data taken with a partly completed array Andeen (2011); Abbasi et al. (2013b); however, the present results are still limited by the amount of data on hand, the systematic uncertainty due to detector effects (particularly the light yield in the ice), and the dependence on the choice of hadronic interaction model used for the simulations. For future analyses, we plan to include more years of experimental data, to simulate more intermediate elements, to investigate new composition-sensitive parameters currently under development, and to incorporate results from new internal studies to reduce the detector systematic uncertainties. These updates will improve the precision of both analyses, and enable the extension of the analyses to higher and lower energies. Furthermore, the analyses presented here are well-suited to capitalize on future extensions to the IceCube Neutrino Observatory Aartsen et al. (2014b).

Acknowledgements.

The IceCube collaboration acknowledges the significant contributions to this manuscript from Marquette University, the University of Alaska Anchorage, and the University of Gent. This paper also profited enormously from the contributions of our departed colleague Stefan Westerhoff (1967-2018). USA – U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany – Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia – Australian Research Council; Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark – Villum Fonden, Danish National Research Foundation (DNRF), Carlsberg Foundation; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea – National Research Foundation of Korea (NRF); Switzerland – Swiss National Science Foundation (SNSF); United Kingdom – Department of Physics, University of Oxford.

.4 Table of Results

The values of the all particle energy spectrum from the coincident analysis, and the corresponding statistical and detector-related systematic uncertainties, as shown in Figure 9, are listed in Table 3. The values of the all particle energy spectrum from the coincident analysis, and the corresponding statistical and detector-related systematic uncertainties, as shown in Figure 12, are listed in Table 4. The values for the elemental fluxes, as shown in Figure 16, are listed in Tables 5 and 6.

.5 K.D.E. templates and fit results

As discussed in Section III.4, the individual templates for each energy bin are shown in Figure 26, while the fit to the data for each energy bin is then shown in Figure 27. Correlation matrices for the fit results for all energy bins are shown in Table 7.