Search for PeV Gamma-Ray Emission from the Southern Hemisphere with 5 Yr of Data from the IceCube Observatory

An event is recorded in IceTop whenever at least three stations (six tanks) are hit by a shower front within a time interval of 6 μs. Then the data recorded from the collection of tanks for each event are filtered to remove uncorrelated background particle hits. The simultaneous reconstruction of shower size and arrival direction is carried out through the maximization of likelihood functions describing the lateral distribution of the signal and the shower front shape (Abbasi et al. 2013). The signal charge distribution, S, around the shower core in the shower frame of reference is known as the lateral distribution function (LDF). The LDF chosen to describe air shower events in IceTop is an empirically derived double logarithmic parabola. At a lateral distance R from the shower axis, the charge expectation S in an IceTop tank is defined as

$$\begin{eqnarray}&&S(R)={S}_{125}{\left(\displaystyle \frac{R}{125{\rm{m}}}\right)}^{-\beta -\kappa {\mathrm{log}}_{10}(R/125{\rm{m}})},\end{eqnarray} \tag{ 1 }$$

where S₁₂₅, also referred to as shower size, is the fitted signal strength at a reference distance of 125 m; β is the fitted slope parameter correlated with shower age; and κ = 0.303 is a constant determined through simulation studies. The shower size, S₁₂₅, is proportional to the energy of the primary particle and converted to a reconstructed energy, E_reco, using parameterization obtained from cosmic-ray simulations (Rawlins & Feusels 2015).

Accumulation of snow on top of the IceTop tanks suppresses the electromagnetic portion of the air shower, requiring a correction factor to the signal expectation defined as

$$\begin{eqnarray}&&{S}_{i}^{\mathrm{corr}}={S}_{i}\,\exp \left(\displaystyle \frac{{h}_{s}^{i}}{{\lambda }_{s}\cos \theta }\right),\end{eqnarray} \tag{ 2 }$$

where ${S}_{i}^{\mathrm{corr}}$ is the corrected signal expectation, hⁱ_s is the snow height above the tank i, θ is the reconstructed zenith angle of the shower, and λ_s is the effective absorption length in snow. The value for λ_s used for each analysis year is also listed in Table 1. The effective absorption length was optimized by comparing each year's snow-corrected S₁₂₅ spectrum with the spectrum from the first year of operation with the least amount of snow. Snow on top of IceTop tanks is constantly increasing at a rate of ∼20 cm yr^–1, which significantly degrades the detector acceptance. In order to accurately account for this effect, the detector response was simulated for each year of data included in the analysis using the same set of CORSIKA gamma-ray showers. For each data set, the snow level on top of the tanks was set to the snow heights measured during the austral summer season. Figure 3 shows the effective area of IceTop to gamma-rays for each year as a function of the primary energy, illustrating the event rate loss caused by increasing amounts of snow. The effective area is lower for the 2011 data year due to low statistics for events with less than eight stations because only one in three such events were being transmitted north from the South Pole in 2011.

Zoom In Zoom Out Reset image size

To ensure good energy and direction reconstruction, the following quality cuts were applied.

• 1.  
  Events must trigger at least five IceTop stations (i.e., where both tanks in the station have DOMs with deposited charge within 1 μs).
• 2.  
  The energy and directional reconstructions must converge.
• 3.  
  The reconstructed core position must be contained within the surface array geometry. Specifically, the event must satisfy D/d > 1, where D and d are defined in Figure 4.
• 4.  
  The tank with the largest deposited charge must not lie on the edge of the IceTop array.
• 5.  
  At least one IceTop tank must have a signal greater than six VEMs.
• 6.  
  The reconstructed zenith angles satisfy cos(θ) > 0.8.
• 7.  
  The reconstructed shower sizes satisfy ${\mathrm{log}}_{10}({S}_{125})\gt -0.25$.

We limit the event selection to ${E}_{\mathrm{reco}}\leqslant 100$ PeV, as extending to higher energies would have required additional higher-energy gamma-ray simulations for a negligible improvement in sensitivity.

Zoom In Zoom Out Reset image size

Angular resolution, defined as the angular radius that contains 68% of reconstructed showers coming from a fixed direction, drives the sensitivity of searches for pointlike and extended sources. Figure 5 shows the angular resolution of simulated gamma-rays as a function of the true primary energy. The first and last years are shown to illustrate the impacts of the snow accumulation on the angular resolution, which amounts to a degradation of around 8%, on average.

Zoom In Zoom Out Reset image size

In order to extract all of the shower information correlated with the type of the primary particle, both surface and in-ice components of the detector are utilized. Section 3.3.1 details the technique used to obtain clean data from the in-ice array that informs on the number of high-energy muons in the shower. Section 3.3.2 describes the implementation of a new likelihood method that optimally retrieves information from the surface array correlated with the number of low-energy muons, shower age, and shower profile. IceTop-based shower discrimination allows us to include showers that do not pass through the in-ice array in the analysis. While the loss of in-ice information for these showers reduces the separation power, there is a large increase in detector acceptance at higher zenith angles, which is displayed in Figure 6. We combine the in-ice and surface components in a single classifier using machine learning as described in Section 3.3.3.

Zoom In Zoom Out Reset image size

3.3.1. High-energy Muons in Ice

Muons that have energy greater than ∼460 GeV at the surface are capable of generating light that can be detected by the photomultipliers of the in-ice array. This is the main parameter used to judge the hadronic content of air showers with trajectories that pass through the in-ice detector component. We use the total charge measured in-ice as a discriminating feature. In order to isolate the signal produced by muons, a cleaning procedure is applied to remove charge deposited by uncorrelated background particles. The cleaning procedures for events with or without an in-ice trigger are described below.

Nearest or next-to-nearest neighboring DOMs on the same string that have deposited charge (a "hit") within a time window of ±1 μs are designated to be in hard local coincidence (HLC). An in-ice trigger is defined as eight or more such HLC hits in an event within a 5 μs time window. An in-ice trigger that falls between 3.5 and 11.5 μs after the start of an IceTop trigger is considered coincident, in which case, hits outside of the coincident time window are removed. Next, HLC hits are used as seeds in a hit selection algorithm that searches for single hits that are within 150 m and 1 μs of each seed DOM. In an iterative manner, single hits that satisfy the criteria are included in the seeds, and the search is performed again. This is repeated three times. The combined charge of the final selection of hits is used as a discriminating feature.

For those events that have only an IceTop trigger, simpler cleaning is applied. For HLC hits, a time window selection is used to reduce the uncorrelated muon background, optimized to be between 3.5 and 9 μs after the trigger in IceTop. Hadronic showers may produce isolated hits in IceCube without an HLC flag. This is illustrated in Figure 7, which displays the number of single hits as a function of vertical DOM number (analogous to depth) and time with respect to the trigger for a collection of events that have no HLC hits. There is a clear excess of hits near the top of the detector that is suitable for classification. A selection window is optimized by a maximization of the ratio of signal and square root noise. We define the front of the time window for charge selection as

$$\begin{eqnarray}&&{t}_{\mathrm{start}}=\displaystyle \frac{4.8\mu {\rm{s}}\ +{d}_{\mathrm{DOM}}/c}{\cos (\theta )},\end{eqnarray} \tag{ 3 }$$

where θ is the reconstructed zenith angle of the event, d_DOM is the depth of the hit in meters, and c is the speed of light. Hits are retained if they fulfill the following criteria, represented by the green box shown in Figure 7.

• 1.  
  The hit is within 130 m of the reconstructed shower axis.
• 2.  
  The hit is within the top 16 layers of in-ice DOMs.
• 3.  
  The hit has a time stamp t, such that ${t}_{\mathrm{start}}\lt t\lt {t}_{\mathrm{start}}+1.8\,\mu {\rm{s}}$.

Zoom In Zoom Out Reset image size

3.3.2. IceTop Shower Footprint

As discussed in Section 2, minimum ionizing muons passing through an IceTop tank deposit charge such that the peak of the muon charge distribution is at ∼1 VEM with a width attributed to the zenith angle distribution of muons (Abbasi et al. 2013). At a sufficiently large distance from the shower core, the electromagnetic component becomes subdominant, and the characteristic ∼1 VEM signals from GeV muons can be discerned. Figure 8 shows a probability distribution function (PDF) describing the distribution of the charges as a function of the lateral distance from the reconstructed shower core, also called the LDF, for simulated gamma-rays and observed cosmic rays. The prominent muon signal can be seen emerging in the cosmic-ray PDF beyond ∼200 m, while it is very diminished in the gamma-ray PDF. The local charge fluctuations, observed as the width of the charge distribution for a given lateral distance in Figure 8, are also a measure of the hadronic content of the shower. The longitudinal development stage of the shower is reflected in the slope of the LDF seen in Figure 8. The curvature of the shower front, i.e., the arrival time distribution of particles as a function of the lateral distance, is also sensitive to the longitudinal stage of the shower as it reaches the IceTop surface.

Zoom In Zoom Out Reset image size

We construct three two-dimensional PDFs that incorporate these shower front properties. For this, we use information from individual IceTop tanks, which are indexed from 1 ≤ i ≤ 162. The PDFs are constructed using tank charges {Q_i}, their lateral distances from the reconstructed shower axis {R_i}, and hit times with respect to the expected planar shower front arrival time {ΔT_i}. Gamma-ray simulations and 10% of cosmic-ray data are used to construct the PDFs for the gamma-ray {H_γ} and cosmic-ray {H_CR} hypotheses, respectively. Unhit and inactive tanks are included in the PDFs by assigning artificial and fixed values to charge and time (Q_i = 0.001 VEM and ΔT_i = 0.01 ns) outside the range of hit tank values. The lateral distance distribution of unhit tanks (as seen in the bottom of the plots in Figure 8) also contributes to the differences between the gamma-ray and cosmic-ray PDFs.

Based on each of the three two-dimensional PDFs, a log-likelihood ratio is calculated for all events. For instance, the log-likelihood ratio using the lateral charge distribution for a given event is defined as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{QR}}={\mathrm{log}}_{10}\left(\displaystyle \frac{{L}_{\mathrm{QR}}(\mathrm{event}| {H}_{\gamma })}{{L}_{\mathrm{QR}}(\mathrm{event}| {H}_{\mathrm{CR}})}\right),\end{eqnarray} \tag{ 4 }$$

where the likelihood L_QR is defined as

$$\begin{eqnarray}&&{L}_{\mathrm{QR}}(\mathrm{event}| H)=\prod _{{\rm{i}}=1}^{162}P\left({Q}_{{\rm{i}}},{R}_{{\rm{i}}}\right|H),\end{eqnarray} \tag{ 5 }$$

with $P\left({Q}_{{\rm{i}}},{R}_{{\rm{i}}}\right|H)$ being the probability of observing a tank with measured charge Q_i and at lateral distance R_i for the hypothesis H. Hit tanks for a sample cosmic-ray event, overlaid on the PDFs in Figure 8 using open squares, illustrate how such an event would collect a greater likelihood from the cosmic-ray PDF as compared to the gamma-ray PDF.

Similarly, one can calculate ${{\rm{\Lambda }}}_{Q{\rm{\Delta }}T}$ and ${{\rm{\Lambda }}}_{{\rm{\Delta }}{TR}}$ from the PDFs that describe the time distribution of charges and the shower front curvature. The sum of all three log-likelihood ratios is then used as an input to a random forest classifier described in Section 3.3.3. The hadronic content and the longitudinal development stage of the shower at the surface depend on the primary energy and zenith angle, in addition to the mass of the primary particle. To reduce this dependence, the construction of PDFs and calculation of the log-likelihood ratio is carried out in ${\mathrm{log}}_{10}({S}_{125})$ bins of 0.1 and $\cos (\theta )$ bins of 0.05.

3.3.3. Random Forest

The final event selection is performed using a random forest classifier implemented using the open-source Python software Scikit-learn (Pedregosa et al. 2011). In total, five features are included in the training process.

• 1.  
  The total charge deposited in the in-ice array after applying the cleaning procedure described in Section 3.3.1.
• 2.  
  The likelihood sum as described in Section 3.3.2.
• 3.  
  The energy proxy log₁₀(S₁₂₅).
• 4.  
  The cosine of the reconstructed zenith angle.
• 5.  
  A parameter that describes the containment of the shower axis within the in-ice array, defined to be C = R/r, where R and r are defined in Figure 4.

To prevent overtraining, the maximum tree depth is limited to 8. The output of the random forest is a score between zero and 1, with 1 being the most gamma-ray-like. The signal threshold for classifying events as gamma-rays is chosen to be 0.7. These values were optimized through a cross-validation grid search on sensitivity performance. Tuning the additional hyperparameters of the random forest showed negligible impact, and they were kept at their default values.

Due to the differences in flux expectation, the gamma-ray spectrum used for training the classifier is different for the point source and diffuse cases. For point sources, in order to be robust in performance to a range of spectral indices, we use two classifiers: one trained with a relatively hard (E^−2.0) spectrum and the other with a relatively soft (E^−2.7) spectrum. In this case, the event is retained if either classifier returns a score above the signal threshold value. For the diffuse case, a single random forest trained with an ${E}^{-3.0}$ spectrum is used (see Section 5.4 for the motivation behind the choice of spectral index). Figure 9 illustrates the fraction of events that pass the signal threshold cut as function of energy. The total number of events classified as gamma-rays under these criteria for both cases are listed in Table 1 for each analysis year.

Zoom In Zoom Out Reset image size

Sources that are pointlike or extended in TeV gamma-ray astronomy should appear pointlike in this analysis. Hence, we construct a point-source hypothesis for an unbiased source search in our entire FOV, as well as for targeted H.E.S.S. source searches. The likelihood under this assumption takes the form

$$\begin{eqnarray}\begin{array}{rcl}L & = & \displaystyle \prod _{j}\,\displaystyle \prod _{i\in j}\left(\displaystyle \frac{{n}_{s}^{j}}{{N}^{j}}{S}_{i}^{j}\left(\left|{{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{S}\right|,{E}_{i},{\sigma }_{i};\gamma \right)\right.\\ & & \quad \left.+\,\,\left(1-\displaystyle \frac{{n}_{s}^{j}}{{N}^{j}}\right){B}_{i}^{j}\left({\delta }_{i},{E}_{i}\right)\right).\end{array}\end{eqnarray} \tag{ 6 }$$

The likelihood L is a product over i events in each of j data sets, where each data set is comprised of 1 yr of data. For a data set j, n_s^j is the number of signal events originating from the point source, and N^j is the total number of events. Each event has a direction ${{\boldsymbol{x}}}_{i}=({\alpha }_{i},{\delta }_{i})$ consisting of an R.A. α_i and decl. δ_i, an energy E_i, and an angular uncertainty ${\sigma }_{i}$. The events are compared to a point-source hypothesis comprised of a direction ${{\boldsymbol{x}}}_{S}$ and spectral index γ. For the single source case, the signal PDF S^j_i is defined as

$$\begin{eqnarray}&&{S}_{i}^{j}=\displaystyle \frac{1}{2\pi {\sigma }_{i}^{2}}{{\rm{e}}}^{-\tfrac{{\left|{{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{S}\right|}^{2}}{2{\sigma }_{i}^{2}}}{{ \mathcal E }}_{S,i}^{j}\left({E}_{i},{\delta }_{i},\gamma \right),\end{eqnarray} \tag{ 7 }$$

where the angular uncertainty is included using a Gaussian distribution with σ_i. Here ${{ \mathcal E }}_{S,i}^{j}$ is the normalized signal energy distribution. The background PDF is defined as

$$\begin{eqnarray}&&{B}_{i}^{j}=\displaystyle \frac{1}{2\pi }{B}_{\exp }^{j}\left({\delta }_{i}\right){{ \mathcal E }}_{B,i}^{j}\left({E}_{i},{\delta }_{i}\right),\end{eqnarray} \tag{ 8 }$$

where B_exp^j is the decl.-dependent detector acceptance to cosmic rays derived from data, and ${{ \mathcal E }}_{B,i}^{j}$ is the normalized background energy distribution. The background PDF is uniform in R.A. and constructed from cosmic-ray data randomized in R.A.

The likelihood is maximized with respect to n_s and γ, where n_s is the total number of signal events proportionally distributed among the signal events n_s^j of each data set according to the effective area and live time of the samples. This yields the best-fit values ${\hat{n}}_{s}$ and $\hat{\gamma }$ and a test statistic defined as

$$\begin{eqnarray}&&{TS}=-2\mathrm{log}\left(\displaystyle \frac{L({n}_{s}=0)}{L({\hat{n}}_{s},\hat{\gamma })}\right).\end{eqnarray} \tag{ 9 }$$

To evaluate the significance of an observed test statistic, a background test statistic ensemble is constructed from scrambled data using random R.A. values for the event directions. The p-value is the fraction of the ensemble that has a test statistic exceeding the observed value. When relevant, the number of independent trials performed (e.g., the number of H.E.S.S. source locations tested individually) is accounted for, and a posttrial p-value is reported for the search.

To gauge the analysis sensitivity to point sources, a range of simulated fluxes is injected on top of the scrambled data events. The sensitivity is defined as the flux that produces a test statistic above the median of the background-only trial ensemble at the injected direction in 90% of the trials. This is equivalent to the Neyman 90% confidence level construction (Neyman 1941). The discovery potential is defined as the flux that achieves a 5σ detection in 50% of the trials.

When searching for emission from the selected H.E.S.S. point sources, we include a test for signal from all sources combined. This stacking approach requires a modification to the likelihood, which we implement following the method in Aartsen et al. (2017a). For a catalog of M point-source locations, the signal PDF is constructed as

$$\begin{eqnarray}&&{S}_{i}^{j}=\displaystyle \frac{{\sum }_{m}^{M}\tfrac{{R}^{j}({\delta }_{m})}{2\pi {\sigma }_{i}^{2}}{{\rm{e}}}^{-\tfrac{{\left|{{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{S}\right|}^{2}}{2{\sigma }_{i}^{2}}}}{{\sum }_{m}^{M}{R}^{j}({\delta }_{m})}{{ \mathcal E }}_{S,i}^{j}\left({E}_{i},{\delta }_{i},\gamma \right),\end{eqnarray} \tag{ 10 }$$

where R_m^j is the relative detector acceptance to gamma-rays at the location of the source m. In this form, the sources are weighted assuming equal flux at Earth for each source.

The source hypotheses for the Galactic plane and cascade neutrino (see Section 5.3) searches extend spatially over a significant portion of the sky. For these cases, it is no longer valid to treat the signal as having a negligible contribution to the background PDF by averaging the R.A. Instead, a modification to the likelihood is made, following a method first introduced by Aartsen et al. (2015a) in a binned likelihood approach and later applied in an unbinned likelihood by Aartsen et al. (2017c), whose formulation we use here.

The background term B_i^j in Equation (6) is replaced with two terms, ${\tilde{D}}_{i}^{j}$ and ${\tilde{S}}_{i}^{j}$, which are the event densities of the experimental data and gamma-ray simulation, respectively, after averaging over R.A.:

$$\begin{eqnarray}\begin{array}{rcl}L & = & \displaystyle \prod _{j}\,\displaystyle \prod _{i\in j}\left(\displaystyle \frac{{n}_{s}^{j}}{{N}^{j}}{S}_{i}^{j}({{\boldsymbol{x}}}_{i},{\sigma }_{i},{E}_{i};\gamma )+{\tilde{D}}_{i}^{j}(\sin {\delta }_{i},{E}_{i})\right.\\ & & \left.-\,\,\displaystyle \frac{{n}_{s}^{j}}{{N}^{j}}{\tilde{S}}_{i}^{j}(\sin {\delta }_{i},{E}_{i}\right).\end{array}\end{eqnarray} \tag{ 11 }$$

The construction of the signal PDF S begins with a model of the gamma-ray flux distribution over the entire FOV. This raw signal term must be convolved with the detector response to produce the expected observed distribution of emission in direction and energy. This is executed through a bin-by-bin multiplication of the relative detector acceptance to gamma-rays, determined through simulation, and the flux model. The point-spread function of each event, described as a Gaussian distribution of width σ, is accounted for through the convolution of the signal map with a range of σ from 01 to 10 in steps of 005. On an event-by-event basis, the map corresponding to the σ of the event is used as the signal PDF S.

An unbiased search for a point source is accomplished by scanning over the entire FOV. The position of a single point source is assumed to lie in the direction of a given pixel in a HEALPIX map (N_side = 512, pixel diameter of 011; Gorski et al. 2005). A test statistic is calculated using Equation (6) under this hypothesis. This process is repeated for each pixel in the analysis FOV. The resulting p-values of this scan, before accounting for trials, are shown in Figure 13. The region of the sky with a zenith angle >5° is excluded, as within this region, scrambling in R.A. alone is insufficient to build independent background trials.

The hottest spot in the sky, with a pretrial p-value of 4 × 10⁻⁵, is located at −734 in decl. and 1484 in R.A., with n_s = ${67.9}_{-16.6}^{+17.8}$ and a spectral index of ${2.9}_{-0.3}^{+0.3}$. The posttrial p-value, calculated by comparing the observed test statistic to the background ensemble of hottest-spot test statistic values, is 0.18, consistent with the background expectation. Figure 10 shows the sensitivity and discovery potential to point sources of this analysis as a function of decl. The proportion of showers which are contained in IceCube drops sharply at higher shower inclinations, resulting in a decrease in performance at high decl. values.

Zoom In Zoom Out Reset image size

There are a total of 15 TeV gamma-ray sources in the FOV of this analysis that have no evidence of a cutoff in their energy spectrum; all of these sources were reported by the H.E.S.S. collaboration. Table 2 lists the name of each source paired with the citation that provided the H.E.S.S. fit information, along with the best-fit decl. and spectral index observed by H.E.S.S. The projected flux at 2 PeV of these sources, assuming there is no cutoff in their energy spectra, is shown in Figure 10. Attenuation effects were calculated for the Galactic coordinates and distance of each source using model results provided by Lipari & Vernetto (2018). Source distances were estimated from associated X-ray or radio observations when possible (Wakely & Horan 2008). Otherwise, a distance of 8.5 kpc, roughly that of the Galactic center, is assumed. The directions of the sources are shown overlaid on the all-sky scan in Figure 11.

Zoom In Zoom Out Reset image size

In the first test of these sources, each source location is evaluated independently, making no assumption about the spectral index of the source. The pretrial p-value for every considered source resulting from this test is reported in Table 2. The most significant individual source, HESS J1427–608, has a pretrial p-value of 0.07, with a fitted n_s of ${25.9}_{-15.7}^{+16.6}$ and spectral index of ${3.2}_{-0.7}^{+0.8}$. This results in a trial-corrected p-value of 0.65, consistent with the background expectation.

To place constraints on the fluxes of these sources individually, we consider the case where the flux of each source extends with no cutoff in energy and with the spectral index observed by H.E.S.S. Under this flux assumption, we report the 90% confidence level upper limit on the flux at 2 PeV at Earth for each source in Table 2 as Φ_90%. We note that uncertainties in the best-fit H.E.S.S. spectrum have a negligible effect on the upper limits and their relationship to the extrapolated flux.

For a subset of the sources, our limits are strong enough to place constraints on an energy cutoff of the source flux. Here we assume the flux to be the combination of a power law and an exponential energy cutoff at E_cut:

$$\begin{eqnarray}&&{\rm{\Phi }}(E)={{\rm{\Phi }}}_{{\rm{P}}.{\rm{L}}.}(E){e}^{-E/{E}_{\mathrm{cut}}}.\end{eqnarray} \tag{ 12 }$$

The value of E_cut was determined by evaluating, for a range of energy cutoff values, the resulting flux of the source and the analysis sensitivity for the same hypothesis. Here we do incorporate absorption from Lipari & Vernetto (2018) in order to limit the energy cutoff at the source. The minimum energy cutoff value where our sensitivity lies below the flux expectation of the source is reported in Table 2 as a 90% confidence level upper limit on E_cut.

We additionally search for combined PeV emission from all of the H.E.S.S. sources. In this test, all of the sources are stacked using the modified signal PDF from Equation (10), which weights the sources assuming equal flux at Earth for each source. The result of this stacking correlation test is a p-value of 0.08, consistent with the background.

From the 4 yr high-energy starting event (HESE) neutrino sample (Aartsen et al. 2015c), a total of 11 events lie within the FOV of this analysis. The event positions and 1σ uncertainties are overlaid on the pretrial p-value map of the all-sky point-source search in a polar projection in Figure 14. There are two topologically distinct event types, referred to as cascades and tracks. We treat the event types separately.

Cascade neutrino events are produced by neutral-current interactions, as well as charged-current interactions of electron and tau neutrinos. These events have poor angular resolution, typically greater than 10° in the HESE sample. In order to correctly account for the change in detector acceptance over such a large point-spread function, we test for a correlation with these cascade neutrino events using the template likelihood method described in Section 4.2. The signal template is constructed by combining the spatial likelihood contour PDFs for each event. The correlation test resulted in a conservative p-value of >0.5. We place a 90% confidence level upper limit of 1.07 × 10⁻¹⁹ cm⁻² s⁻¹ TeV⁻¹ on the flux at 2 PeV of a source class consistent with the HESE cascade directions and an ${E}^{-2.0}$ spectrum.

Track neutrino events are the product of charged-current muon neutrino interactions. These interactions produce muons that can travel several kilometers through the ice, which allows for reconstructions with an angular resolution of <12 for events in the HESE sample. There is one track event in our FOV, which we treat under the point-source prescription. At a decl. of δ = −8677, scrambling events in R.A. alone is insufficient to build a background test statistic distribution; instead, we scramble events randomly throughout the polar cap region of >5° zenith angle, within which acceptance was found to be approximately even. No evidence of signal was found in the test for a point source at the event location, which resulted in a conservative p-value of >0.5.

To test for a diffuse flux from the Galactic plane, the template likelihood method is employed, with the signal template taken to be the pion decay component of the Fermi-LAT diffuse emission model (Ackermann et al. 2012). The Fermi-LAT template multiplied by the detector acceptance is shown in Figure 15 for the 2012 sample.

The expected observed spectral index of diffuse gamma-rays from the Galactic plane is still an open question due to existing uncertainties in the Galactic cosmic-ray spectrum, interstellar gas distribution, and flux attenuation. The unattenuated flux at PeV energies has been predicted to be as hard as E⁻³ (Ingelman & Thunman 1996; Vernetto & Lipari 2017), although the spectrum could be as soft as E^−3.4 (Ingelman & Thunman 1996) after attenuation. Here we fix the spectral index to 3.

Maximization of the likelihood in Equation (11) returns an observed test statistic corresponding to a p-value of 0.28, which provides no evidence of a diffuse signal from the Galactic plane under the current hypothesis. We place a 90% confidence level upper limit of 2.61 × 10⁻¹⁹ cm⁻² s⁻¹ TeV⁻¹ on the normalization of the spectral energy distribution at 2 PeV with spectral index 3. The normalization energy was chosen to be 2 PeV, as that is the energy for which the analysis was found to be least sensitive to the spectral index assumption.

As previous experimental limits have used a box region around the Galactic plane to describe the diffuse emission, for a limit comparison, we use an angular-integrated scaled flux,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{template}}={\rm{\Phi }}{\rm{\Delta }}{\rm{\Omega }}\,\displaystyle \frac{{\int }_{\mathrm{all}-\mathrm{sky}}{S}_{\mathrm{Fermi}}d{\rm{\Omega }}}{{\int }_{{\rm{\Delta }}{\rm{\Omega }}}{S}_{\mathrm{Fermi}}d{\rm{\Omega }}},\end{eqnarray} \tag{ 13 }$$

where the angular-integrated flux from the observed region is ΦΔΩ, and the second term scales this flux by the fraction of the Galactic plane emission present in the observed region as given by the Fermi-LAT pion decay template S_Fermi. This is, in effect, a limit on the overall normalization of the Fermi template flux. Figure 12 compares the scaled flux limits from this analysis with the existing IC-40 (Aartsen et al. 2013a) and CASA-MIA (Borione et al. 1998) limits, as well as the measured flux by ARGO-YBJ (Bartoli et al. 2015). In addition, the gamma-ray flux from the diffuse Galactic plane in the analysis FOV is shown for two model predictions from Lipari & Vernetto (2018). The first assumes that cosmic-ray spectra have a space-independent spectral shape throughout the Galaxy. The second assumes that the central part of the Galaxy has a harder cosmic-ray spectrum than observed at Earth, an idea supported by some recent analyses (Gaggero et al. 2015b; Yang et al. 2016).

Zoom In Zoom Out Reset image size

There are a number of systematic uncertainties that can affect the estimated sensitivity of the study, including the snow attenuation, the calibration of the charge deposited in IceTop tanks, the hadronic interaction and ice models used in simulations, and the anisotropy of the cosmic-ray flux. For each systematic component, data sets are constructed using the parameter uncertainty bounds. These data sets are propagated through the entire analysis chain to evaluate the corresponding uncertainty in sensitivity. Here we report both point-source and Galactic plane analysis sensitivity uncertainties. For the point-source sensitivity, we assume the median uncertainty based on a scan over the decl. range in 01 bins. These data sets were constructed from events not used in the training of the random forest classifier detailed in Section 3.3.3.

The optimal value for the snow attenuation coefficient λ, used in the charge correction for IceTop tanks as shown in Equation (2), has been found to vary by ±0.2 m over a range of zenith angle and energy values (Aartsen et al. 2013b). A more thorough model of the snow attenuation is a work in progress; for this work, we quote this bound as a conservative estimate. This systematic results in an uncertainty of 11.0% in sensitivity to point sources and 11.2% in sensitivity to a diffuse flux from the Galactic plane.

The calibration of the IceTop tank charge to VEM units requires a fitting of the muon peak in the charge spectrum of the tank (Abbasi et al. 2013). The dependence of this fit to systematic factors was studied in detail by Van Overloop (2011). They found an uncertainty of at most ±3% on the charge calibration, which propagates directly to an uncertainty in the deposited signal. This systematic error results in an uncertainty of 2.1% in sensitivity to point sources and 7.4% in sensitivity to a diffuse flux from the Galactic plane.

The number of muons generated in simulated gamma-ray air showers at energies sufficient to trigger the detectors is governed by the high-energy hadronic interaction model used in CORSIKA. In order to evaluate the magnitude of the model dependence, we perform sensitivity studies with simulations generated using QGSJetII-04 (Ostapchenko 2011) but otherwise identical to the original set. We chose QGSJetII-04 over other post–Large Hadron Collider models because it was the model that produced the most muons in hadronic air showers (Plum et al. 2018). Sensitivities calculated with these systematic data sets resulted in a 23.2% uncertainty for point sources and a 26.2% uncertainty for a diffuse flux from the Galactic plane.

The anisotropy of the cosmic-ray flux is a potential source of signal contamination. While decl.-dependent anisotropy is accounted for due to the use of data to construct the background PDF in the likelihood, any anisotropy in R.A. is not. However, within the analysis FOV, this anisotropy is at a level of at most 0.03% (Aartsen et al. 2016). This is negligible in relation to statistical uncertainties, which are ∼25% in flux at the sensitivity threshold.

In simulation, the uncertainties in the optical properties of the ice can affect the amount of charge measured. While this has a potential for systematic error, an analysis with data sets using ±10% in deposited charge in IceCube showed negligible impact on sensitivity compared to statistical fluctuations. Finally, the method we use naturally corrects for any bias in the energy proxy. Any systematic biases in fitted n_s and γ values that are unaccounted for are found to be negligible when compared to statistical uncertainties.

Under the assumption that the errors discussed are independent and Gaussian-distributed, the overall sensitivity uncertainty resulting from quadrature addition is 25.8% for point sources and 29.4% for the Galactic plane.