Search for Sources of Astrophysical Neutrinos Using Seven Years of IceCube Cascade Events

The procedure for rejecting the atmospheric muon background depends on the event topology of interest. Neutrino-induced muon tracks with energies ≳1 TeV originating in the northern sky can be selected with high efficiency and low atmospheric muon contamination by identifying events reconstructed at declinations δ ≳ 5° with high confidence, as only neutrinos can travel through so much intervening earth and/or ice before producing muons that pass through the detector. Neutrino- and cosmic-ray-induced muon tracks originating in the southern sky and entering the detector from above can only be distinguished probabilistically, and only under the assumption that the neutrino spectrum is harder than the atmospheric muon spectrum. Thus the energy threshold increases to ∼100 TeV in the southern sky, resulting in weaker sensitivity, especially for a soft neutrino spectrum.

In this work, we instead turn our attention to cascade events produced when the neutrino interaction vertex, and hence first observed light, occurs inside the detector. With this approach we accept all neutrino flavors and most interaction types, approximately independent of decl., while efficiently rejecting downgoing atmospheric muons. An added benefit for astrophysical neutrino searches is that for declinations ≲−30° the atmospheric neutrino background is naturally suppressed because many are accompanied by incoming atmospheric muons originally produced in the same cosmic ray shower in the upper atmosphere (Schonert et al. 2009).

Most Cerenkov light from a muon traveling through ice is radiated through stochastic processes, resulting in a dense, linear series of cascade-like signatures that may be observed in our detector. The mean distance between these energy deposits decreases with increasing energy. For energies ≳60 TeV, incoming muons can be rejected with high confidence using a veto region consisting of just the outermost DOMs, reserving the majority of the instrumented region as a fiducial volume for neutrino detection (Aartsen et al. 2014b). To lower the threshold to ∼1 TeV while holding the incoming muon rejection rate constant, the thickness of the veto region must be increased. Below we summarize this method, which is used as described in Wandkowsky & Weaver (2018) and which further optimizes the approach first introduced in Aartsen et al. (2015a).

We begin with all events passing one or more basic filters at the South Pole. A splitting algorithm is applied to each event, identifying ∼75% of unrelated but temporally coincident physical events initially merged in the DAQ output by clustering causally connected sets of pulses. We reject any event in which the first ≥3 pulses appear in the outer layer veto region as described in Aartsen et al. (2013a). An additional veto is applied to reject events in which two or more PEs are observed consistent with a downgoing track passing through the interaction vertex or a major energy deposition. Finally, a cut is applied on the interaction vertex location, scaling with observed charge as described in Aartsen et al. (2015a) such that at 100 PEs the fiducial volume is reduced to just the DeepCore sub-array, while at ≥6000 PEs the fiducial volume consists of all but the outermost layer of DOMs. This final cut enables efficient background rejection down to ∼1 TeV by keeping the probability of observing veto photons approximately independent of energy.

We rely on a traditional maximum likelihood method (Aartsen et al. 2014a) to obtain initial reconstructions used for cascade/track discrimination. The goal of this reconstruction is to unfold the spatial and temporal pattern of energy depositions for each event. Two fits are performed: one that is constrained to find a single dominant cascade-like energy deposition, and one that finds a linear combination of such energy depositions distributed along a possible muon track. Events in which at least 6000 PEs were collected are classified as tracks if the free track fit finds at least two non-negligible depositions more than 500 m apart, or if the free track fit is associated with more charge than the single cascade fit. Events with less total collected light are classified as tracks if at least 1.5 PEs are consistent with an outgoing muon track originating at the reconstructed interaction vertex (Wandkowsky & Weaver 2018). All other events are classified as cascades and are used in the present analysis.

The selection criteria described above were applied to data taken from 2010 to 2017 May as well as to neutrino and atmospheric muon Monte Carlo (MC) simulations used for performance estimates. The first year of data comes from the nearly complete 79-string configuration while the remaining six years make use of the complete 86-string detector. In a total of 2428 days of IceCube livetime, 10,422 events survive until cascade/track discrimination; of these, 1980 are identified as cascades. Note that while the dominant improvement in this data set is the increase from two to seven years of data, the neutrino effective area is also enhanced by applying coincident event splitting and veto criteria to data from every initial South Pole filter. This increases the acceptance by 23% (67%) for a signal following an E⁻² (E⁻³) spectrum.

From MC simulations, we find that 98% of truly cascade-like events that pass all selection criteria are correctly identified as such. The rate at which CC muon neutrino interactions are successfully classified as track events increases with energy as more light is produced by the outgoing muon. For a conventional atmospheric neutrino spectrum, 30% of the cascade channel consists of misclassified CC muon neutrino interactions; for an astrophysical spectrum following E^−2.5 or harder, this contribution reduces to 5% or less. This population of misclassified events results in a tolerable background at lower energies as well as a small signal contribution at higher energies.

Because muon track analyses specifically target events with high-quality track reconstructions and reject events dominated by individual cascade-like energy depositions, we expect the cascade analysis to be largely statistically independent in spite of the small but nonzero misclassification rate. In fact, the final cascade selection shares just a single ∼2 TeV event in common with the latest muon track selection.

In past work, we have used a maximum likelihood method to reconstruct neutrino energy and direction of travel from IceCube cascades (Aartsen et al. 2014a). This approach relies on detailed parameterizations of the position- and direction-dependent light absorption and scattering lengths in the ice, neither of which is large compared to the DOM spacing. This results in a complex multi-dimensional likelihood function with many local optima in the R.A. and decl. coordinates (α, δ), such that it is computationally expensive to find the global optimum for any given event and prohibitive to estimate the per-event statistical uncertainties.

In this work, we introduce a novel cascade reconstruction using a deep NN. An NN is a highly flexible function mapping from an input layer to an output layer via a series of hidden layers, where each successive layer consists of a set of values computed based on the values contained in the previous layer. The functional forms of the layer-to-layer connections (the network architecture) must be designed a priori; the numerical parameters of those connections are optimized through a training procedure to yield good results for a given training data set. NNs are well-suited to problems in high-energy physics for which we are typically able to generate high-statistics MC data sets for use in training.

Our NN-based reconstruction draws from recent advances in image recognition and is implemented using Tensorflow (Abadi et al. 2015). The network architecture used here is largely the same as that introduced previously for muon energy reconstruction (Huennefeld 2018). The method will be described in detail in a separate publication, but here we will outline the main considerations relevant to this analysis.

IceCube data consist of a set of waveforms (represented as a series of pulse arrival times and PE counts) accumulated over time on a number of DOMs distributed throughout the three-dimensional instrumented volume, and thus is in general four-dimensional. Our first step is to compute waveform summary values for use in the input layer. For each DOM, these values consist of the relative time of the first pulse, the time elapsed until 20%, 50%, and 100% of the total charge is collected, the total charge collected, the charge collected within 100 and 500 ns of the first pulse, and the charge-weighted mean and standard deviation of relative pulse arrival times.

The detector is divided into three sub-arrays: IceCube, lower DeepCore, and upper DeepCore. Each sub-array is independently well-approximated by a regular spatial grid suitable for processing by several initial convolutional layers, which are able to exploit symmetries in the structure of the input data to facilitate efficient network optimization and usage⁵⁸ (see Huennefeld 2018 for diagrams of the relevant geometry). The output from the convolutional layers is taken as the input for each of two fully connected NNs (in which each node in a given layer is connected to every node in the preceding layer). One of these networks is optimized to estimate the physical parameters of interest—the R.A., decl., and energy (α, δ, E)—while the other is optimized to estimate the uncertainties on these parameters.

All training was performed using 50% of the relevant signal MC, with the remaining 50% reserved for testing analysis-level performance. Two training passes were performed. The first pass made use of several MC data sets: one with baseline values for key parameters such as DOM quantum efficiency and light absorption and scattering lengths, and several more with modified values within estimated systematic uncertainties. In addition to offering overall increased training statistics, the use of these differing data sets may give the NN some robustness against known systematic uncertainties. The second training pass refined the network to give the smallest errors and, on-average, unbiased reconstructions for the baseline MC. In each pass, a priori per-parameter weighting was applied such that angular resolution is valued over energy resolution by a factor of 5.

The expected performance of the NN angular reconstruction (including systematics; see Section 5.2) is shown as a function of energy in Figure 1. Compared to the reconstructions used in our previous analysis of two years of data (Aartsen et al. 2017d), the NN offers significantly improved angular resolution above 10 TeV (a factor of 2 improvement at 1 PeV). While we do not recover the optimal statistics-limited angular resolution described in Aartsen et al. (2014a), we do obtain performance that improves monotonically with increasing energy up to ∼1 PeV. At higher energies, the estimated systematic uncertainty becomes large enough to prevent any further improvement. Note that an additional advantage of the NN angular reconstruction used here is that it naturally provides per-event uncertainty estimates usable in the statistical analysis described in Section 5.1, whereas previous work relied on a parameterization of typical uncertainties derived from signal MC.

Zoom In Zoom Out Reset image size

The performance of the energy reconstruction is comparable to that used in previous work. The estimated energy is within 60% of the true neutrino energy for 68% of events, averaged over all neutrino flavors and interaction types, and approximately independent of spectrum. This performance estimate, like the sensitivities quoted in Section 5.3, assumes a flavor ratio of 1:1:1 with equal contributions from ν and $\bar{\nu }$, detected via a mixture of CC and NC interactions.

The energy and decl. distributions of cascade events in data are compared with neutrino and atmospheric muon MC in Figure 2. The distributions obtained are similar to those observed in the two year sample (Aartsen et al. 2015a, 2017d).

Zoom In Zoom Out Reset image size

One way to search for astrophysical neutrino sources with only a minimal set of a priori assumptions about source position is to search the entire sky for the most significant point-like neutrino clustering in excess of the background expectation on a dense grid of pixels that are small compared to the neutrino angular resolution. This approach has most recently been employed by IceCube using tracks (Aartsen et al. 2017a) and cascades (Aartsen et al. 2017d) as well as by ANTARES using tracks and cascades in combination (Albert et al. 2017a), and we include it in the present analysis as well. However, an all-sky scan is subject to a large trial factor and thus is in general less sensitive compared to analyses that use prior information to restrict the set of hypothesis tests.

An alternative approach is to scan only the positions of a modest number of well-motivated source candidates, which substantially reduces the trial factor. In addition, where multiple analyses report results for the same or overlapping catalogs, direct comparisons can be made. Here we scan the same catalog of 74 source candidates that was studied in the previous IceCube cascade paper (Aartsen et al. 2017d).

We consider one source in more detail: the supermassive black hole at the center of the Galaxy, Sagittarius A*. Based on hints from gamma-ray observations (e.g., Herold & Malyshev 2019), there may be emission up to some unknown high-energy cutoff from a spatially extended region centered approximately on this object. Therefore we evaluate constraints on the flux from this region as a function of possible spatial extension and for several possible spectral cutoffs.

The gamma-ray blazar TXS 0506+056 does not appear in the a priori catalog described above. In light of this, and in anticipation of future identifications of unexpectedly promising source candidates based on neutrino observations, we treat this object as a monitored source to be studied separately from the catalog scan described above.

For source classes for which we can predict approximate relative signal strengths, it may be possible to increase the signal-to-background ratio using a source-stacking method (e.g., Abbasi et al. 2011). Because the present analysis offers good sensitivity in the southern sky, roughly independent of possible spatial extension up to a few degrees, we include stacking analyses for three Galactic supernova remnant (SNR) catalogs derived from SNR Cat (Ferrand & Safi-Harb 2012) and previously studied using IceCube tracks (Aartsen et al. 2017b). These SNRs are categorized based on their environment: those with associated molecular clouds, those with associated pulsar wind nebulae (PWNs), and those with neither. The angular extension of these objects reaches up to 163, and each catalog comprises a preponderance of objects in the southern sky.

Cosmic ray interactions with interstellar gas in the Milky Way are expected to produce neutral and charged pions, where neutral pions would decay to observable gamma-rays and charged pions would yield potentially observable neutrinos. The hadronic gamma-ray emission up to 100 GeV has been identified by Fermi-LAT using a multi-component fit (Ackermann et al. 2012). A corresponding neutrino flux prediction can be obtained by extrapolating this measurement to energies above 1 TeV in the context of Galactic cosmic ray production and propagation models.

The original model fits by Fermi somewhat under-predict the measured gamma-ray flux in the Galactic plane, and especially near the Galactic center, above 10 GeV. The ${\mathrm{KRA}}_{\gamma }$ models obtain better agreement with gamma-ray data in this regime by introducing Galactocentric cosmic ray diffusion parameter variability and an advective wind (Gaggero et al. 2015, 2017). Model-dependent neutrino flux predictions are provided assuming cosmic ray injection spectra with exponential cutoffs at 5 PeV or 50 PeV per nucleon; we refer to these as ${\mathrm{KRA}}_{\gamma }^{5}$ and ${\mathrm{KRA}}_{\gamma }^{50}$, respectively.

The latest constraints on diffuse Galactic neutrino emission depend on the ${\mathrm{KRA}}_{\gamma }$ models and were obtained in a joint IceCube and ANTARES analysis (Albert et al. 2018) which made use of complementary features of the IceCube track analysis (Aartsen et al. 2017b) and ANTARES track and cascade combined analysis (Albert et al. 2017b). In this work we search for emission following ${\mathrm{KRA}}_{\gamma }^{5}$ as the primary diffuse Galactic emission result; we also test for emission following ${\mathrm{KRA}}_{\gamma }^{50}$. Finally, we test for emission following the spatial profile of the Fermi-LAT π⁰-decay measurement, assuming an E^−2.5 neutrino energy spectrum.

The Fermi bubbles consist of a pair of gamma-ray emission regions that extend to ∼55° above and below the Galactic center (Su et al. 2010). Most of the Fermi bubble region yields a relatively hard gamma-ray spectrum up to ∼100 GeV, with some evidence for spectral softening above that energy (Ackermann et al. 2014). The gamma-ray emission has been speculated to be of hadronic origin (Crocker & Aharonian 2011), powered by cosmic ray acceleration in the vicinity of the Galactic center; however, the true origin of the Fermi bubbles has not yet been experimentally identified.

We derive constraints on emission from the Fermi bubbles following spectra of the form ${dN}/{dE}\propto {E}^{-2.18}\cdot \exp (-E/{E}_{\mathrm{cut}})$, for ${E}_{\mathrm{cut}}\in \left\{50\,\mathrm{TeV},100\,\mathrm{TeV},500\,\mathrm{TeV}\right\}$—the same spectra tested in recent work by ANTARES (Hallmann & Eberl 2018). If there is neutrino emission from the Fermi bubbles with a significantly softer spectrum or lower cutoff energy, this analysis would not be sensitive to it.

In this work we consider two broad categories of source candidates: point-like and extended template, where the latter include diffuse Galactic emission and emission spanning the Fermi bubbles. Both analysis types are based on the standard likelihood (Braun et al. 2008) given by a product over all events i in the data set

$$\begin{eqnarray}&&{ \mathcal L }({n}_{s},\gamma )=\displaystyle \prod _{i}\left[\displaystyle \frac{{n}_{s}}{N}{{ \mathcal S }}_{i}({{\boldsymbol{x}}}_{i}| \gamma )+\left(1-\displaystyle \frac{{n}_{s}}{N}\right){{ \mathcal B }}_{i}({{\boldsymbol{x}}}_{i})\right],\end{eqnarray} \tag{ 1 }$$

where N is the total number of events, n_s is the expected number of signal events, γ is the signal spectral index, and ${{\boldsymbol{x}}}_{i}$ represents the event R.A., decl., angular uncertainty, and energy $\left\{{\alpha }_{i},{\delta }_{i},{\sigma }_{i},{E}_{i}\right\};$ ${{ \mathcal S }}_{i}({{\boldsymbol{x}}}_{i}| \gamma )$ is the probability density function (PDF) assuming event i is part of the signal population, and ${{ \mathcal B }}_{i}({{\boldsymbol{x}}}_{i})$ is the PDF assuming event i is part of the atmospheric or unrelated astrophysical background populations. For all source types, n_s is free to vary between 0 and N. For point-like sources, the signal spectral index γ is free to vary between 1 and 4, while for extended templates γ is fixed to a source-dependent constant value (γ = 2.5 for diffuse Galactic emission and γ = 2.18 for emission from the Fermi bubbles).

The details of our signal and background likelihoods, ${{ \mathcal S }}_{i}$ and ${{ \mathcal B }}_{i}$, follow established methods applied previously to IceCube tracks for individual (Aartsen et al. 2017a) and stacked (e.g., Abbasi et al. 2011) point-like sources as well as for extended templates (Aartsen et al. 2017b). We do not require a specialized treatment, in contrast to our previous cascade analysis (Aartsen et al. 2017d), thanks to increased statistics in the experimental data set as well as new per-event angular uncertainty estimates given by the NN reconstruction.

As in previous work, we define the test statistic as the log-likelihood ratio ${ \mathcal T }=-2\mathrm{ln}\{{ \mathcal L }({n}_{s}=0)/{ \mathcal L }({\hat{n}}_{s},\hat{\gamma })\}$, where ${\hat{n}}_{s}$ and $\hat{\gamma }$ are the values that maximize ${ \mathcal L }$, subject to the constraints specified above. This test statistic is used to compute significances, sensitivities, discovery potentials, and upper limits (ULs). For the all-sky (source candidate catalog) scan, we compute a post-trials significance based on the most significant pixel (source candidate) tested, in order to guarantee the reported false positive rates. Sensitivities (90% CL), UL (90% CL), and discovery potentials (5σ) are defined as in our previous analysis (Aartsen et al. 2017d) and are computed using the Neyman construction (Neyman 1937).

The dominant systematic uncertainties in this analysis include the optical properties of the ice, the quantum efficiency of the DOMs, and the neutrino interaction cross section. These uncertainties affect the angular resolution and the signal acceptance. As in our previous cascade analysis (Aartsen et al. 2017d), we treat these effects as approximately separable. However, we have improved our approach to each consideration; we describe our latest method in the following.

The NN reconstruction is trained to yield optimal performance on baseline MC; the angular resolution for real data events will be somewhat worse. To estimate how much worse, we perform dedicated simulations of events similar to those observed, but using depth-dependent ice model variations intended to cover the uncertainties in the model. By comparing the median resolution from these modified simulations with that from the baseline MC, we obtain a function of energy that quantifies how much worse the resolution may be than expected from the baseline. This factor ranges from 10% at 1 TeV to ∼50% at 2 PeV, and is taken as a correction to the angular separation between the reconstructed and true direction for each event in the baseline MC. This factor is similarly applied to the angular uncertainty estimates σ_i for both MC and data events. In this way, we directly account for systematic uncertainties impacting angular resolution in the quantiles shown in Figure 1 as well as in all p-values and sensitivity flux calculations in the analysis.

The above treatment accounts for the analysis-level impact of systematic uncertainties for each observed event. To address the uncertainties in the detection efficiency, and thus in sensitivity, discovery potential, and upper limit fluxes, we compute the energy-integrated signal acceptance, as a function of decl. and for each considered spectrum, based on additional MC data sets produced with varied modeling assumptions (the same modified data sets used in NN training; see Section 3.2). We find that, for plausible ice model and detector variations, the signal acceptance variation ranges from ∼10% for an unbroken E⁻² spectrum to ∼17% for E⁻² with an exponential cutoff at 100 TeV, roughly independent of decl.. As was done in the previous analysis, we estimate an uncorrelated 4% impact from uncertainties in the neutrino interaction cross section. These values are added in quadrature on a per-spectrum basis to obtain a final estimate of uncertainties via signal acceptance effects. In the remainder of this paper, all sensitivity, discovery potential, and upper limit fluxes include this factor.

All sensitivities discussed in the remainder of this paper are per-neutrino flavor (assuming a flavor ratio of 1:1:1 at the detector), but summed over ν and $\bar{\nu }$. The point source sensitivity flux as a function of source decl. is shown in Figure 3 for several spectral scenarios: unbroken power laws following hard (γ = 2) and soft (γ = 3) spectra, and spectral cutoff scenarios ${dN}/{dE}\propto {E}^{-2}\cdot \exp (-E/{E}_{\mathrm{cut}})$ with ${E}_{\mathrm{cut}}\in (100\,\mathrm{TeV},1\,\mathrm{PeV})$. Where published values are available for previous IceCube work with tracks (Aartsen et al. 2017a) or cascades (Aartsen et al. 2017d), or for the most recent ANTARES track and cascade combined analysis (Albert et al. 2017a), these are shown for comparison. We find that the present analysis improves upon the previous IceCube work with cascades at all declinations and across the tested spectra, with the largest improvements reaching a factor larger than 4 in the southern sky. Furthermore, we now obtain the best sensitivity of any analysis for hard sources in the southernmost ∼30% of the sky ($\sin (\delta )\lt -0.4$). This search also achieves sensitivity comparable to that of ANTARES for spectra with cutoffs as low as E_cut = 100 TeV, but with much weaker decl. dependence.

Zoom In Zoom Out Reset image size

The sensitivities of the SNR stacking analyses are listed in Table 1. In this work we obtain a sensitivity below previously set ULs (Aartsen et al. 2017b) only for the SNR-with-PWN catalog, which consists of eight southern SNRs and one northern SNR. It is nevertheless interesting to revisit all three catalogs here because, while they all include southern source candidates, in previous work the results necessarily were dominated by northern candidates due to the strongly decl.-dependent signal acceptance of the IceCube track selection.

The sensitivities of the diffuse Galactic template analyses are listed in Table 2. This analysis obtains ∼30% (40%) better sensitivity to ${\mathrm{KRA}}_{\gamma }^{5}$ (${\mathrm{KRA}}_{\gamma }^{50}$) than the recent joint IceCube+ANTARES analysis (Albert et al. 2018). Compared to the IceCube analysis using seven years of tracks (Aartsen et al. 2017b), this analysis obtains ∼15% better sensitivity to emission following the spatial profile of the Fermi-LAT π⁰-decay measurement. These improvements are possible because the expected emission follows a soft (γ ∼ 2.5) spectrum and is concentrated near the Galactic center at δ ∼ −30°, where IceCube track analyses are subject to a large background of atmospheric muons but the present cascade analysis efficiently rejects this background as well as some of the atmospheric neutrino background; the improvement is larger for the ${\mathrm{KRA}}_{\gamma }$ models than for the Fermi-LAT π⁰ model because the former are specifically tuned to increase the concentration of the expected flux near the Galactic center.

Table 2.  Sensitivity and Results of the Diffuse Galactic Template Analyses, Compared to Latest Previous Work

Template	7 yr Cascades				Previous Work
	p-value	Sensitivity	Fitted Flux	UL	p-value	Sensitivity	Fitted Flux	UL
${\mathrm{KRA}}_{\gamma }^{5}$	0.021	0.58	0.85	1.7	0.29	0.81	0.47	1.19
${\mathrm{KRA}}_{\gamma }^{50}$	0.022	0.35	0.65	0.97	0.26	0.57	0.37	0.90
Fermi-LAT π0	0.040	2.5	3.0	6.5	0.37	2.97	1.28	3.83

Note. Sensitivity, fitted flux, and ULs are given as multiples of the model prediction for ${\mathrm{KRA}}_{\gamma }$ models, and as ${E}^{2}\cdot {(E/100\mathrm{TeV})}^{0.5}\cdot {dN}/{dE}$ in units 10⁻¹¹ TeV cm⁻² s⁻¹ for Fermi-LAT π⁰ decay. The latest previous results are from a joint IceCube-ANTARES analysis (Albert et al. 2018) for ${{\rm{K}}{\rm{R}}{\rm{A}}}_{\gamma }$ and from IceCube tracks (Aartsen et al. 2017b) for Fermi-LAT π⁰ decay.

Download table as:  ASCIITypeset image

The sensitivity flux for the Fermi bubble analyses is ∼30% below the UL shown in Figure 8, approximately independent of spectral cutoff. This analysis obtains sensitivity that is at least one order of magnitude better than the recent ANTARES search (Hallmann & Eberl 2018), with the improvement increasing with spectral cutoff energy, E_cut. Because we assume an even more extended template than ANTARES, covering a total solid angle of about 1.18 sr compared to ∼0.66 sr, this factor is even larger if considered in terms of flux per solid angle. Once again, this improvement is due to efficient rejection of atmospheric backgrounds for the cascade data set used in this work.

Zoom In Zoom Out Reset image size