1 Introduction

The first clear evidence for the existence of an invisible mass component in the universe was Zwicky’s observation of the dynamics of the Coma galaxy cluster [1]. Subsequently, a broad range of cosmological and astrophysical observations supported the existence of this dark matter (DM) on various scales, from galaxy cluster scales down to galactic scales. Measurements of galactic velocity profiles hint at invisible mass distributed beyond the visible disks [2]. Galaxy cluster dynamics exhibit a similar behavior [3].

Further evidence for the existence of dark matter can be found in galaxy cluster mergers like the Bullet Cluster [4, 5]. Following a collision, the interstellar and intergalactic gas components as seen in X-ray observations are spatially separated from the mass distribution reconstructed by weak lensing. Such a separation strongly disfavors theories of modified gravity.

According to the current understanding of the formation and evolution of large-scale structures, cold (non-relativistic), or warm dark matter is preferred over hot (relativistic) dark matter. Otherwise, the formation of the observed large-structures on time scales of the order of the age of the universe would not have been possible [68].

Though the nature of dark matter is largely unknown, some of its properties may be deduced from the above-mentioned observations. Analyses of temperature fluctuations in the Cosmic Microwave Background (CMB) by the Planck collaboration [8] yield the current best estimate for the total content of DM in the universe: \(\Omega _\mathrm{CDM} h^2 = 0.1199 \pm 0.0027\), with the cold DM density parameter \(\Omega _\mathrm{CDM}\), and \(h=0.673 \pm 0.012\) being the Hubble parameter divided by \(100\,\mathrm {km/s\,Mpc}\).

Besides inference from gravitational interaction, particle DM may also be detected indirectly. A weakly interacting massive particle (WIMP) at roughly GeV-scale masses is a favorable class of DM; it naturally provides the right order of magnitude for the thermal relic abundance of DM in the early universe [9]. Examples of WIMPs are neutralinos in supersymmetric extensions of the Standard Model [10], or the lightest stable excitations in Kaluza-Klein models [11].

If DM decays, or (self-)annihilates, a flux of stable final-state messenger particles, e.g. charged leptons, photons, and neutrinos, may be detected at Earth, making DM experimentally accessible by indirect searches (e.g. [1214]). The neutrino is an attractive messenger particle because it propagates without absorption, and neutrino vacuum oscillations do not alter the energy and direction information. Further, no fore- or background from astrophysical objects has been identified yet.

Regions of increased DM density, like massive celestial objects, dwarf galaxies, galactic halos, and the Galactic Center, provide targets to search for an increased flux of neutrinos. Due to its proximity, the Galactic Center is expected to yield the highest flux of annihilation products [15].

While most of these sources would appear as (nearly) point-like sources in the sky, the Galactic Center is an extended source, and a signal from the Milky Way halo would lead to a large-scale anisotropy in neutrino arrival directions  [16, 17]. With its 4\(\pi \) acceptance, the IceCube neutrino detector [18], is well-suited for DM searches from all of the above-mentioned sources.

In this paper we present the results from a search for a neutrino signal from DM self-annihilation in the Galactic Center, targeting DM masses ranging from \(30~\mathrm {GeV}\) to \(10~\mathrm {TeV}\). Due to the wide range of event topologies associated with neutrinos from this energy range, two event selections are motivated and presented. One event selection focuses on the low-mass region ranging from 30 to 100 GeV, accessible through the low-energy in-fill array DeepCore (DC) [19], with the surrounding parts of IceCube used as veto. The second event selection focuses on the mass range 100 GeV–1 TeV, but extends up to 10 TeV. For this selection a larger part of the IceCube detector is defined as fiducial volume. Throughout this paper we denote the low-mass event selection as LE and the high-mass selection as HE.

2 Dark matter halos

Dark matter halos are considered to be gravitationally self-bound overdensities of DM particles, formed through hierarchical merging of proto-halos from primordial density fluctuations [20]. There is a tension between halo profile fits to DM overdensities in N-body simulations, and fits to observational data (the cusp-core problem) [21]. N-body simulations seem to favor cuspy halos, while observations of e.g. dwarf spheroidal galaxies imply a rather flat central core region. However, the inner part of the halo profile may depend on the host halo mass [22]. Simple models of DM halos describe the density by a smooth spherically-symmetric function of the halo-centric radius, with the maximal density at the center, and a decreasing density with increasing radius. One parametrization of such density profiles is given by (modified from [23])

$$\begin{aligned} \rho _\mathrm{{DM}} (r) = {\rho _0 \over \left( \delta + \frac{r}{r_s} \right) ^\gamma \cdot \left( 1 + \left( \frac{r}{r_s} \right) ^\alpha \right) ^{ (\beta - \gamma )/\alpha } }, \end{aligned}$$
(1)

with the shape parameters \(\alpha \), \(\beta \), \(\gamma \), the scale radius \(r_s\), and the mass density normalization \(\rho _0\), which is usually determined from the assumed local DM density, \(\rho _{\mathrm {local}}\), in our solar system. We introduced the parameter \(\delta \) to allow for a central core if set to 1, while \(\delta =0\) describes a cuspy halo profile.

The parametrization of Eq. (1) is a combination of power laws, where e.g. the power-law index \(\gamma \) describes the inner slope, while \(\alpha \) and \(\beta \) describe the outer slope. The halo-centric distance of this transition region depends on the scale radius \(r_s\).

In view of the unresolved cusp-core problem, we present results for two halo density profiles. The widely used Navarro-Frenk-White (NFW) profile represents cuspy halos [24], and is chosen for comparability among different experimental results. The Burkert profile is chosen as representative of flat-cored profiles [25]. Based on observation, the latter profile is currently favored for the Milky Way [26]. Table 1 shows the parameter values for the two models used in this work.

Table 1 DM Halo parameters used in this analysis. Taken from [26]

3 Neutrino signal from dark matter annihilation

The flux of final state particles from annihilating dark matter depends on the DM mass density squared, integrated along the line-of-sight (l.o.s.) through the DM halo, and is given by the \(J_{\mathrm {a}}\)-factor. Following e.g. [13, 15], the \(J_{\mathrm {a}}\)-factor is

$$\begin{aligned} J_{\mathrm {a}}(\Psi ) = \int \limits _0^{l_\mathrm{{max}}} \mathrm {d}l ~ \rho _{\text{ DM }}^2 \left( \sqrt{ R^2_{\text{ SC }} - 2l R_{\text{ SC }} \cos {\Psi } +l^2}\right) . \end{aligned}$$
(2)

Here, the density profile along the l.o.s. is parametrized for a given angle between the l.o.s. and the direction of the center of the galaxy, \(\Psi \). The parameters are the radius of the solar circle, \(R_{\text{ SC }} \approx 8.5~\mathrm {kpc}\), and the maximal distance from the observer along the l.o.s., \(l_\mathrm{{max}}\). The latter is

$$\begin{aligned} l_\mathrm{{max}} = \sqrt{R^2_{\text{ MW }} -R^2_{\text{ SC }} \sin ^2{\Psi } } + R_{\text{ SC }} \cos {\Psi }, \end{aligned}$$
(3)

with the assumed radius of the Milky Way, \(R_{\text{ MW }} \approx 50~\mathrm {kpc}\). Typically, radii larger than the scale radius do not contribute significantly to the total value of \(J_{\mathrm {a}}\). Figure 1 (top panel) shows \(J_{\mathrm {a}}\) for the NFW (solid line) and Burkert (dashed line) profiles.

Fig. 1
figure 1

Top: Line-of-sight integral \(J_{\mathrm {a}}(\Psi )\) of the NFW (solid) and Burkert (dashed) DM profile using parameters as given in Table 1. Bottom: Example of DM annihilation spectra generated with PYTHIA8 [27] for a WIMP mass of \(m_\chi = 500\) GeV. Three annihilation channels are shown: \(b\bar{b}\) (solid), \(W^+W^-\) (dashed), and \(\mu ^+\mu ^-\) (dotted). The neutrino line is not shown, but it is modeled as a box at the WIMP mass with a width of 10 % of the wimp mass

The final differential neutrino flux from DM annihilation, \(\mathrm {d}\phi _\nu / \mathrm {d}E\), depends on the neutrino energy spectrum of the actual annihilation channel. The differential neutrino flux is

$$\begin{aligned} \frac{\mathrm {d}\phi _\nu }{\mathrm {d}E} = \frac{\left<\sigma _\mathrm {A} v\right>}{2} ~ \frac{ 1}{4\pi m_{\chi }^2} ~\frac{\mathrm {d}N_\nu }{\mathrm {d}E}~ J_{a}(\Psi ), \end{aligned}$$
(4)

with \(\left<\sigma _\mathrm {A} v\right>\) being the thermally averaged product of self-annihilation cross-section, \(\sigma _\mathrm {A}\) and WIMP velocity v. Further, \(m_\chi \) is the WIMP mass, \(\mathrm {d}N_\nu / \mathrm {d}E\) is the neutrino energy spectrum per annihilating WIMP pair, and \(J_{\mathrm {a}}\) is the DM abundance along the l.o.s..

We consider several benchmark annihilation channels with 100 % branching ratios (\(\chi \chi \rightarrow b \bar{b}\), \(W^+ W^-\), \(\mu ^+ \mu ^-\), \(\tau ^+ \tau ^-\),\(\nu \bar{\nu }\)) for the calculation of the neutrino energy spectrum \(\mathrm {d}N_\nu / \mathrm {d}E\). The resulting spectra bracket the realistic annihilation neutrino energy spectra with a mixture of different annihilation branching ratios. The annihilation spectrum to neutrino pairs is approximated by a box at the WIMP mass with a width of 10 % of the WIMP mass. For the other channels we generated a neutrino energy spectrum from annihilating DM using the PYTHIA8 (version 8.175) software package [27]. Our PYTHIA8 simulation was set up to simulate a generic resonance with an energy of twice the WIMP mass forming only the final state particle pair (i.e. \(b\overline{b}\), \(W^+ W^-\), etc.) in question. Subsequent processes like hadronization and decay were simulated using the default PYTHIA8 implementations. The generic resonance ensures an isotropic decay of weak bosons, e.g. in the \(W^+ W^-\) channel. Thus, the spin of the annihilating WIMPs is not considered and we don’t assume a specific WIMP model like the lightest neutralino described by supersymmetric models [9]. If the WIMP is indeed the lightest neutralino, the spin (1/2) of the WIMP would affect the generation of the neutrino energy spectra. The spin of such a WIMP would lead to fully transversely polarized W-bosons in the final state of the annihilation process, thus altering the neutrino energy spectrum [28]. The differences in the differential neutrino yield compared to the isotropically decaying W bosons is about \(\pm 40~\%\). Examples of the neutrino energy spectra used here are shown in the bottom panel of Fig. 1 for the \(b\overline{b}\), \(W^+W^-\), and \(\mu ^+\mu ^-\) annihilation channels.

In general, neutrinos are subject to neutrino oscillations on the way from the Galactic Center to the Earth. Due to the very long baseline, we assume a relative neutrino flavor ratio of 1:1:1 at Earth.

Fig. 2
figure 2

A footprint (top) and side view (bottom) of IceCube in detector coordinates. The green circles mark regular IceCube DOMs with 17 m vertical spacing, while the yellow diamonds mark DeepCore DOMs with about half the regular vertical DOM spacing and higher quantum efficiency PMTs. The black squares mark DOMs which are part of the final 86-string configuration, but are not present in the here used 79-string configuration of IceCube. The purple shaded and red shaded areas illustrate the fiducial volumes used by the LE and HE event selections, respectively

These simulations have the detector response folded in and are generated using the ANIS event generator [29] modelling the neutrino-nucleon charged and neutral current interactions via CTEQ5 [30] parton distributions for neutrino and anti-neutrino interactions.

Finally, the neutrino energy distributions are used to weight generic simulated neutrino data to DM annihilation signal.

4 The IceCube neutrino observatory

The IceCube Neutrino Observatory, situated at the geographic South Pole, consists of an in-ice detector array, IceCube, and a surface air shower detector array, IceTop [31], dedicated to neutrino and cosmic ray research, respectively. IceCube [32] is installed in the glacial ice at depths between 1450 and 2450 m below the surface, instrumenting a total volume of one cubic kilometer. IceCube detects neutrinos by optical detection of Cherenkov radiation induced by secondary charged leptons which are produced in neutrino interactions in the surrounding ice or the nearby bedrock.

Construction of the IceCube detector started in the Austral summer of 2004. In January 2010, 79 detector strings with 60 digital optical modules (DOMs) on each string [32, 33] were deployed. Each DOM contains a 25 cm Hamamatsu photomultiplier tube (PMT) and on-board electronics to readout and digitize the signal from the PMT [34]. In December 2010, the IceCube detector construction was completed. The final IceCube detector consists of 86 strings.

A schematic layout of the detector is shown in Fig. 2. The 79 strings used in this analysis are marked by green and yellow markers within the outer shaded grey area. The square markers denote the additional strings constituting the completed IceCube detector. Of the 79 strings used in this analysis, 73 strings (green) have a horizontal spacing of 125 m and a vertical spacing of 17 m between DOMs. The six remaining strings (yellow) are located near the central string of IceCube. The DOMs on these strings are equipped with PMTs with a 30 % increased quantum efficiency. Together with their nearest IceCube strings, these strings constitute the inner detector, DeepCore [19] (for IceCube-79). The vertical distance between DOMs is reduced to 7 m (10 m) for the bottom 50 (upper 10) DOMs. The horizontal distance between strings in DeepCore is less than 75 m.

These two densely-instrumented parts (see Fig. 2) are separated by a region with significantly reduced scattering and absorption lengths for Cherenkov photons due to dust particles. It is located at a depth of about 2050 m.

For muon-neutrino events, the neutrino arrival directions are inferred from the muon arrival direction. The latter is reconstructed using a likelihood approach, based on the arrival times of photons at DOMs [35]. Thus, a good understanding of the absorption and scattering of photons is necessary for direction reconstruction. The clean glacial ice is a natural medium, built up over tens of thousands of years. Thus, the optical properties exhibit a variation over the 1 km depth of the instrumented volume. A detailed description of the optical ice properties is given in [36].

5 Event selection

The signal for this analysis is muons produced in charged-current neutrino interactions. These muons produce track-like event signatures in the detector, which allow for a reconstruction of the arrival direction.

At the South Pole, the Galactic Center is always \(29^\circ \) above the horizon. Thus, neutrinos from the direction of the Galactic Center will appear as down-going events in IceCube. The backgrounds for this analysis are therefore down-going atmospheric muons and, at a lower rate, muons produced by neutrinos originating from cosmic-ray showers in the atmosphere. The overwhelming majority of the 2500 Hz trigger rate is due to atmospheric muons. The atmospheric neutrino background contributes to the trigger rate at the 1 mHz-level. This event class is an irreducible background, since in the energy range of interest for this analysis the accompanying muon component of the atmospheric shower is absorbed in the ice sheet above the detector.

The approach adopted here to reduce the background is to consider only neutrinos which interact within the detector, and reject the background of penetrating (in-coming) muons. In order to select events which appear to start within the detector we developed several complementary veto techniques, exploiting differences in timing and topology of background and signal events.

This work is based on two independently developed event selections, referred to as low-energy (LE) and high-energy (HE) selections or samples. There are two reasons for such energy-specific optimization. First, the efficiency of these vetoes decreases rapidly with decreasing event energy, because low-energy muons are able to traverse several string-layers without being detected. Second, the likelihood function (described in Sect. 6) does not use the event energy.

Though the individual event selections differ, the general selection techniques and the analysis pipelines are very similar. First, an initial online selection is performed at the South Pole. Second, a set of cuts on the event quality, topology, and arrival directions is applied. Third, a boosted decision tree (BDT) is used to remove remaining background events, using the TMVA software package [37]. The BDT is trained on a representative signal assumption for each sample. Finally, a likelihood analysis is performed, exploiting the different distributions of arrival directions of background events and events originating in dark matter annihilations in the Galactic Center. The differences between the two event selections are highlighted in the following sections.

This analysis uses data collected with IceCube in its 79-string configuration between May 31, 2010 and May 13, 2011 with a total live-time of 319.7 days of stable high-quality data. The LE sample contains 35,538 events, and the HE sample contains 293,043 events. 4706 events appear in both samples; about 13 % of the LE events are in the HE sample, and about 1.6 % of the HE events are in the LE sample. At this level the angular resolution of the event direction ranges from \(\approx 2^\circ \), for high energy neutrinos, to \(\approx 20^\circ \) for the lowest neutrino energies.

5.1 Low-energy event selection

The LE event selection considers events from the DeepCore online-filter [19], and is optimized for low-mass WIMPs below 100 GeV, and thus uses the bottom part of the densely-instrumented DeepCore sub-array as fiducial volume. The remaining instrumented IceCube volume as well as the two bottom DOM layers are used as a veto. The fiducial volume is illustrated in Fig. 2, and corresponds to roughly 27 Mton of ice.

The LE selection cuts are based on experience from the IceCube-79 Solar WIMP analysis [38], which used DeepCore for the first time in low-mass WIMP searches. The signal used for BDT training are events that are fully contained in the fiducial volume, and originate in annihilations of 65 GeV WIMPs to \(b \bar{b}\)-pairs in the NFW halo. The search window for the LE analysis extends to \(\pm 30^\circ \) in right ascension (\(\alpha \)) with respect to the Galactic Center, while the declination (\(\delta \)) width is asymmetric and extends from \(-39^\circ \) to \(-9^\circ \).

The LE sample data rate at the analysis level is 1.4 mHz.

5.2 High-energy event selection

The IceCube array has a trigger energy threshold of about \(\simeq 100\,\mathrm{GeV}\). Therefore, the search for WIMPs in the mass range above a few hundred GeV benefits from the large volume of IceCube in addition to DeepCore at the cost of a decreasing veto efficiency. The HE selection considers events from the dedicated Galactic Center online-filter and the DeepCore online-filter.

The veto for the HE event selection is defined by the upper 12 DOM layers and the two outer string layers, which roughly corresponds to 200 and 125 m of instrumented distance, respectively. The fiducial volume for the HE selection is shown in Fig. 2. The signal assumed for BDT training are events which start in the fiducial volume, and originate in annihilations of 600 GeV WIMPs to \(W^+W^-\)-pairs in the NFW halo. The search window around the Galactic Center is given by \(\pm 15^\circ \) in both declination and right ascension.

The HE sample data rate at the analysis level is 10 mHz.

6 Analysis method and sensitivity

A maximum likelihood analysis is performed independently on each event selection for a number of different WIMP masses ranging from 30 GeV to 10 TeV, assuming a 100 % branching ratio for each tested annihilation channel.

Considering the large number of events in the two final samples, a binned likelihood method was chosen. To reduce the number of bins, event arrival directions were only considered in a search window around the Galactic Center. The search window shape and size differ slightly between the two event samples as defined in the two previous Sects. 5.1 and 5.2 as well as the background estimations. Due to these differences the likelihood analysis performed on the LE and HE selections is not identical. The LE likelihood has the more complicated form and is defined as a function of signal fraction, \(\xi \), in the following way:

$$\begin{aligned} L(\xi ) = \left( {\begin{array}{c}N\\ n\end{array}}\right) p^n(1-p)^{N-n}\prod ^n_{i=1}f(\varvec{X}_i,|\xi ) \end{aligned}$$
(5)

where the binomial factor in front of the product accounts for the probability of observing n events in the search window given N total events in the event selection, and the shape term, i.e the direction, \(\varvec{X} = (\delta , \alpha )\), of the events is accounted for by the term \(f(\varvec{X}_i,|\xi )\). The binomial probability was chosen in favor of a Poisson probability since the search window covers a non negligible fraction of the declination band. The probability of an event to fall in the search window is defined as \(p = \pi _\mathrm {s}\xi +\pi _\mathrm {bg}(1-\xi )\), where \(\pi _\mathrm {s}\) and \(\pi _\mathrm {bg}\) are the probability for a signal or a background event, respectively, to fall in the search window. Note that the signal originating in a dark matter halo is an extreme case of an extended source as it is present in the whole sky. Any background estimation based on data will be contaminated by signal. \(\pi _\mathrm {bg}\) is determined from the relation between the size of the search window and the size of the background estimation region, while \(\pi _\mathrm {s}\) is determined from simulation.

Fig. 3
figure 3

Example skymaps of the background (top) and signal (bottom) pdfs in equatorial coordinates for the LE event selection, for 100 GeV WIMPs annihilating into \(W^+W^-\). The LE search window is marked with a dashed line. For illustration the search window for the HE event selection is marked with a dash-dotted line

The directional probability density function (pdf), \(f(\varvec{X}|\xi )\), in Eq. (5) is constructed from binned expectations of event directions for background and signal. Figure 3 shows examples of these binned expectations. The software package HEAL-Pix [39] was used to ensure equal-area bins on the sphere for these two-dimensional pdfs.

To determine the signal pdfs, IceCube neutrino simulations were used, weighting events according to Eq. (4) and the corresponding DM annihilation spectrum. Background pdfs were created by scrambling the right ascension of experimental data events in the final event selections.

The same reasoning regarding signal contamination, as stated above, applies to the directional pdf, \(f(\varvec{X}|\xi )\). In effect, the expected arrival directions of background events will depend on the signal strength. This needs to be accounted for in the directional pdf, \(f(\varvec{X}|\xi )\) in the likelihood, as well as in the background simulation during the construction of confidence intervals. The directional pdf is defined as

$$\begin{aligned} f(\varvec{X}|\xi ) = wf_\mathrm {s}(\varvec{X})+(1-w)[(1+w)f_{\mathrm {bg}}(\varvec{X})-wf_{\mathrm {sc}}(\varvec{X})] \end{aligned}$$
(6)

where \(f_\mathrm {s}\) and \(f_\mathrm {bg}\) are the signal and background directional pdfs, respectively, \(f_{\mathrm {sc}}\) is a pdf describing the signal scrambled in right ascension. The signal fraction inside the search window is

$$\begin{aligned} w = \frac{\pi _\mathrm {s}\xi }{\pi _\mathrm {s}\xi +\pi _\mathrm {bg}(1-\xi )}. \end{aligned}$$
(7)

A different approach is used for the HE analysis. The background estimation is performed on off-source data, excluding all events within \(\pm 30^\circ \) of the Galactic Center. Therefore, any signal contamination of the background estimate is ignored, and the likelihood function from Eq. (5) simplifies to:

$$\begin{aligned} L(n_\mathrm {s}) = \frac{ (n_{\text {bg}}+ n_\mathrm {s})^{n} }{ n! }e^{-(n_\text {{bg}}+ n_\mathrm {s})} \prod _i^{n} f(\varvec{X}_i , n_\mathrm {bg}|n_\mathrm {s}), \end{aligned}$$
(8)

where \(n_\mathrm {s}\) is the number of signal events and \(n_\mathrm {bg}\) is the expected number of background events in the search window. \(n_\mathrm {bg}\) is obtained by multiplication of the number of off-source events by the ratio of the on-source and off-source region sizes. The directional pdf consequently becomes:

$$\begin{aligned} f(\varvec{X}, n_\mathrm {bg}|n_\mathrm {s}) = \frac{n_\mathrm {s}}{n_\mathrm {bg}+n_\mathrm {s}}f_\mathrm {s}(\varvec{X})+\frac{n_\mathrm {bg}}{n_\mathrm {bg} +n_\mathrm {s}}f_{\mathrm {bg}}(\varvec{X}) \end{aligned}$$
(9)
Fig. 4
figure 4

Median upper limits, i.e. sensitivities, at 90 % C.L. (w/o systematics) for the two event selections assuming WIMP annihilation to \(\nu \bar{\nu }\) for the NFW DM profile. The black solid line shows the combined best sensitivity for this particular annihilation channel and DM profile, considering both event selections. At 200 GeV the HE selection yields a slightly better sensitivity and is thus used here

All confidence intervals are constructed using the prescription by Feldman and Cousins [40]. The sensitivity is defined as the median upper limit on the number of signal events at 90 % confidence level.

The final limits for each WIMP mass and annihilation channel are obtained from the sample (LE or HE) which gives the best sensitivity (w/o systematics), i.e the lowest median upper limit, for the particular WIMP mass and annihilation channel. Thus, the cross-over point in the WIMP mass between the two event samples depends on the DM halo model and WIMP annihilation channel. This procedure circumvents the necessity to deal with the small overlap of both samples. Figure 4 illustrates how the two event samples contribute to the best sensitivity at different WIMP masses, in this case for WIMP annihilation to \(\nu \bar{\nu }\). Each sample has a WIMP mass range where it outperforms the other.

Figure 5 shows the neutrino effective area for the two event selections. Even though the effective area for the LE event selection is smaller than that of the HE event selection, the sensitivity to the number of signal events for low-mass WIMPs is better, as can be seen in Fig. 4. The reasons are the larger on-source region for the LE event selection, and a lower background event rate due to higher veto efficiency.

In order to avoid confirmation bias throughout the development of the analysis, blindness with respect to the right-ascension information was imposed by scrambling the right-ascension information of the experimental data. Only the declination information of the events was used for cut development. The experimental right ascension distribution was unblinded after fixing all cuts for

Following the optimization of the two individual event selections, meaning all cuts are fixed, the sensitivity was calculated. For each WIMP mass, halo profile, and channel the sample which yielded the best sensitivity was used for the final analysis, and the data was unblinded.

Fig. 5
figure 5

Neutrino effective areas as a function of energy for the two event selections. The effective areas for the LE and HE selection are shown as solid, and dotted lines, respectively. Although the HE effective area is bigger than the LE effective area at low energies, the higher background contamination at low energies in the HE selection makes it less efficient

7 Discussion of uncertainties

The uncertainties relevant for this analysis can be categorized into two classes:

  • Detector systematic uncertainties impacting the signal efficiency

  • Astrophysical uncertainties (choice of halo model, model-specific parameter)

The former are incorporated into the calculated limits, while the latter are studied to estimate model uncertainties. Both classes are discussed in the following sections.

7.1 Detector systematics

The uncertainties in the signal efficiency are mainly governed by the uncertainties in the optical efficiency of the DOMs and the optical properties of the glacial ice, manifested in the absorption and scattering length. To determine the effects on \(\langle \sigma _\mathrm {A} v \rangle \) due to the mentioned uncertainties, the event selections and analysis were applied to sets of simulated data where the optical properties of the DOMs and the ice were changed. The optical efficiency of the DOMs was varied by \(\pm 10\) %. The same was done for the absorption and scattering lengths of the ice. The resulting uncertainties on \(\langle \sigma _\mathrm {A} v \rangle \) generally lie in the range 10–20 % except for the lowest neutrino energies where they reach up to \(\approx 70~\%\). This is due to threshold effects where events with just enough hit DOMs to trigger the detector would fail to do so with increased absorption and scattering or decreased optical efficiency.

The above-described systematic uncertainties are included into the limits by degrading the baseline results by the relative variation of the detector uncertainties with respect to the baseline, as stated above.

7.2 Astrophysical uncertainties

The astrophysical uncertainty is studied by using two different halo profiles, and also by varying the parameters \(\rho _{\mathrm {local}}\) and \(r_\mathrm{s}\) within the uncertainties stated in [26], and summarized in Table 1. Figure 6 compares the sensitivity for WIMPs annihilating to \(\nu \bar{\nu }\)-pairs for both profiles. The bands depict the variation of the sensitivity within each profile that arises from varying the profile parameters within the given uncertainty.

Fig. 6
figure 6

Uncertainties on the sensitivity due to the DM halo model parameter uncertainties for the NFW (dashed line, green band) and Burkert (solid line, blue band) DM profile, assuming WIMP annihilation to \(\nu \overline{\nu }\). The reduced width of the bands below 100 GeV is caused by different on-source regions, and thus differences in the integration of the \(J_{\mathrm {a}}\)-factors. The dip below 100 GeV is caused by the under-fluctuation in the LE sample. The natural scale is the self-annihilation cross-section region for WIMPs to be thermal relics from the Big Bang [41]

The relative variation of \(\langle \sigma _\mathrm {A} v \rangle \) due to the halo profile parameter uncertainties was estimated to be 60–100 % for the LE likelihood analysis and 60–200 % for the HE likelihood analysis, and is shown in Fig. 6. However, this is a conservative approach which overestimate on the uncertainty, since the correlation of the two parameters in the fits is neglected.

Different measurements of the local dark matter density yield a wider range of values, e.g. the PDG group states a canonical value for \(\rho _\mathrm{local}\) of 0.3 \(\mathrm {GeV/cm^3}\) within a factor of 2–3 [42]. The limit depends on \(\rho _\mathrm{local}^2\), thus the factor of 2–3 translates to a variation of the limit by almost an order of magnitude.

We refrain from including uncertainties on the halo profile parameters into the limits.

8 Results

Table 2 shows the number of events for the two unblinded event selections. The quantities \(n_\text {obs}\) and \(n_\text {bg}\) are the number of measured on-source events and expected background events in the search windows, respectively. A \(2\sigma \) under-fluctuation of experimental data events was observed for the LE event selection. A systematic origin of this under-fluctuation due to an uneven right-ascension exposure was excluded. A small over-fluctuation was measured for the number of on-source events in the HE selection. However, the likelihood analysis yields underfluctuations in both cases. resulting in upper limits which are lower than the corresponding sensitivities.

Table 2 The number of observed events in the on-source region, the corresponding number of background events expected from data scrambling, and their difference \(\Delta n = n_\text {obs}-n_\text {bg}\). An under-fluctuation is seen in the LE selection while the HE selection shows a slight excess
Fig. 7
figure 7

Sensitivity (dashed) and observed upper limit (solid, w/o systematics) at 90 % C.L. for WIMPs annihilating to neutrinos assuming a NFW (top) and Burkert (bottom) DM halo profile. The statistical uncertainty on the sensitivity is shown at the \(1\sigma \) (green band) and \(2\sigma \) (yellow band) level. The dip below 100 GeV is caused by the under-fluctuation in the LE sample

For the HE selection this implies that despite the over-fluctuation in the number of events, the spatial distribution of these events within the search window is incompatible with the expectation from dark matter annihilation in the halo. As an example, Fig. 7 shows the sensitivity and the observed upper limit after unblinding on \(\langle \sigma _\mathrm {A} v \rangle \) for the NFW and Burkert profiles, assuming WIMP annihilation into neutrinos. In addition the \(\pm 1\sigma \) and \(\pm 2\sigma \) statistical uncertainty bands of the median upper limit are shown as green and yellow shaded areas, respectively. The contribution from the LE and HE event selection can be clearly seen through the upper limit curve (solid black line) being lower than the expected median upper limit (dashed black line) in both cases, but to a different extent, with the switch-over between \(m_\chi = 100\) and 200 GeV for this annihilation channel.

Table 3 summarizes the upper limits for all considered annihilation channels and WIMP masses. The limits are shown separately for the two considered DM halo profiles; the NFW (top) and Burkert (bottom) profile. Figure 8 shows the sensitivity and limit for three different annihilation channels.

To compare the performance of this analysis to previous IceCube analyses and other experiments, we choose the \(\tau ^+\tau ^-\) annihilation channel assuming a NFW DM halo profile. The comparison is shown in Fig. 9. The black solid line shows the limit of this analysis, whereas dashed lines with markers show the limits from previous galactic halo [16, 17] and dwarf spheroidal galaxies [42] analyses with IceCube. The other lines show the limits from gamma-ray experiments, in particular the limit from the dwarf spheroidal galaxy Segue 1 analysis by VERITAS [43] (dash-dotted) and MAGIC [44] (dash-dot-dotted), and the limit from the Fermi analysis of several dwarf spheroidal galaxies [45] (dashed). Also shown is the DM interpretation of the positron-fraction excess reported by the PAMELA collaboration (dark gray shaded region) and the \(3\sigma \) and \(5\sigma \) preferred regions from the \(e^+ + e^-\)-flux excess reported by the Fermi and H.E.S.S. collaborations as dark green and green shaded regions, respectively. All the shaded region data are taken from [46] and rescaled to a local dark matter density of \(\rho _{\mathrm {local}}=0.471\) GeV cm\(^{-3}\) to this DM halo profile parameter with the one considered in the other analyses. For a WIMP mass below 1TeV the present analysis improves significantly in sensitivity on previous IceCube analyses. Furthermore by using the DeepCore detector array, the self-annihilation cross-section for WIMP masses below 100 GeV is probed for the first time by IceCube.

Table 3 Final upper limits (including detector systematics) on the self-annihilation cross-section, \(\langle \sigma _{A} v \rangle \), for different annihilation channels and WIMP masses, \(m_\chi \), for the NFW (top) and Burkert (bottom) DM halo profiles
Fig. 8
figure 8

Sensitivity (dashed) and observed upper limits (solid) at 90 % C.L., including detector systematics, for WIMPs annihilating to \(b\overline{b}\) (stars), \(\tau ^+\tau ^-\) (squares), and directly to neutrinos (triangles) assuming a NFW DM halo profile. The shaded areas are guides for the reader’s eyes connecting sensitivity and limit for a particular annihilation channel

Fig. 9
figure 9

Comparison of limits from this work (IC79 GC) to other IceCube (designated by IC + “string number”) searches for dark matter annihilation in self-bound structures [16, 17, 42]. Further, photon search limits from observation of dwarf spheroidals by VERITAS [43], MAGIC [44] and Fermi [45] are shown. The grey-shaded region is a dark-matter interpretation of the positron excess reported by the PAMELA collaboration. The green-shaded regions are the \(3\sigma \) and \(5\sigma \) preferred regions from the \(e^++e^-\)-flux excess reported by Fermi and H.E.S.S. All shaded region data taken from [46]. The region data and the IC22 halo limits are rescaled to the here assumed local dark matter density of \(\rho _\mathrm {local}=0.471\) GeV cm\(^{-3}\). The natural scale is the self-annihilation cross-section region for WIMPs to be thermal relics from the Big Bang [41]. The black dotted line in the upper right part of the figure is the unitarity bound [47]. All shown limits and preferable dark matter regions assume a branching ratio of 100 % into \(\tau \bar{\tau }\). We note that preliminary Galactic Center limits from the ANTARES neutrino telescope have recently been released [48]

9 Conclusion

We have presented limits on the cross-section on dark matter annihilation in the Galactic Center, probing down to \(\left< \sigma _\mathrm {A} v \right> \simeq 4\cdot 10^{-24}\mathrm {cm}^3 \mathrm {s}^{-1}\) at 65 GeV WIMP mass, assuming the NFW halo profile and direct annihilation to neutrino pairs.

This analysis is the first IceCube Galactic Center DM search using the nearly complete detector configuration. Further, it is the first IceCube DM search probing \(\langle \sigma _\mathrm {A} v \rangle \) for WIMP masses below 100 GeV by utilizing the DeepCore infill-array of IceCube.

We have presented methods for a selection of down-going muon neutrinos in IceCube, making the southern hemisphere accessible to low-energy neutrino searches in the energy range 10 GeV–10 TeV. These methods have been applied to create two event selections, that are optimized for neutrino signals from the direction of the Galactic Center. Based on these event selections a likelihood analysis looking for a neutrino flux from annihilating dark matter in the Galactic Center was performed, testing a number of dark matter annihilation channels at different masses.

The results are compatible with the background-only hypothesis, thus upper limits on \(\langle \sigma _\mathrm {A} v \rangle \) were set (c.f. Fig. 9). The limits from the low-energy selection are almost \(2\sigma \) lower than their sensitivity due to an under-fluctuation in the number of background events. The limits presented here for direct annihilation to \(\nu \overline{\nu }\)-pairs are model-independent and conservative upper bounds for dark matter annihilation to Standard Model final states [49]; even small branching ratios to other – more visible – species at the \(\langle \sigma _\mathrm {A} v \rangle \)-level presented here would yield a detectable flux in gamma-ray experiments, or otherwise stronger constraints. Thus, these limits complement gamma-ray detection channels.

Future improvements to this analysis can be expected from improvements in the background rejection in the energy region corresponding to the highest probed WIMP masses. Further, the inclusion of an energy term in the likelihood function is expected to increase the senstitivity, especially for line-like features from direct annihilation to neutrino pairs.

Long-term improvements should also be expected from possible IceCube extensions. The low-energy upgrade PINGU [50] would increase the sensitivity to low-mass WIMPs, and extend the probed mass range below 30 GeV. PINGU is a possible future in-fill array with a denser instrumentation than DeepCore. The high-mass (TeV-PeV) sensitivity would benefit from a future high-energy extension, IceCube-Gen2 [51]. The aim for IceCube-Gen2 is an expanded instrumented volume of the order of 10 km\(^3\) with a larger inter-string spacing, compared to IceCube.