Search for Neutrinos from Dark Matter Self-Annihilations in the center of the Milky Way with 3 years of IceCube/DeepCore

We present a search for a neutrino signal from dark matter self-annihilations in the Milky Way using the IceCube Neutrino Observatory (IceCube). In 1005 days of data we found no significant excess of neutrinos over the background of neutrinos produced in atmospheric air showers from cosmic ray interactions. We derive upper limits on the velocity averaged product of the dark matter self-annihilation cross section and the relative velocity of the dark matter particles $\langle\sigma_{\text{A}}v\rangle$. Upper limits are set for dark matter particle candidate masses ranging from 10 GeV up to 1 TeV while considering annihilation through multiple channels. This work sets the most stringent limit on a neutrino signal from dark matter with mass between 10 GeV and 100 GeV, with a limit of $1.18\cdot10^{-23}\text{cm}^3\text{s}^{-1}$ for 100 GeV dark matter particles self-annihilating via $\tau^+\tau^-$ to neutrinos (assuming the Navarro-Frenk-White dark matter halo profile).

Abstract We present a search for a neutrino signal from dark matter self-annihilations in the Milky Way using the Ice-Cube Neutrino Observatory (IceCube). In 1005 days of data we found no significant excess of neutrinos over the background of neutrinos produced in atmospheric air showers from cosmic ray interactions. We derive upper limits on the velocity averaged product of the dark matter self-annihilation cross section and the relative velocity of the dark matter particles σ A v . Upper limits are set for dark matter particle candidate masses ranging from 10 GeV up to 1 TeV while considering annihilation through multiple channels. This work sets the most stringent limit on a neutrino signal from dark matter with mass between 10 and 100 GeV, with a limit of 1.18 · 10 −23 cm 3 s −1 for 100 GeV dark matter particles selfannihilating via τ + τ − to neutrinos (assuming the Navarro-Frenk-White dark matter halo profile).

Introduction
With the increasingly strong indications of the existence of extended halos of dark matter surrounding galaxies and galaxy clusters [1], there is much interest within the particle physics community to determine the nature and properties of dark matter [2]. The frequently considered hypothesis is that dark matter consists of stable massive particles interacting feebly with Standard Model particles. The density of dark matter particles today is determined by the 'freeze-out' [3][4][5][6] in the early universe when the thermal equilibrium can no longer be sustained as the universe expands and cools down. This work focuses on a generic candidate particle for dark matter referred to as a weakly interacting massive particle (WIMP) [7][8][9][10], though this search is sensitive to any self-annihilating dark matter particle with a coupling to the Standard Model resulting in a flux of neutrinos. The source considered is the Milky Way galaxy, which is embedded in a spherical halo of dark matter [11][12][13][14][15]. For a given halo density profile, the total amount of dark matter in the line of sight from Earth can be determined [16].
If WIMPs can self-annihilate into Standard Model particles and the dark matter density is sufficiently high, an excess of neutrinos and photons should be observed from parts of the sky with a large amount of dark matter, above the background of muons and neutrinos produced in the Earth's atmosphere. Although photons produced in such annihilations are far easier to detect, it is still of interest to consider scenarios where only neutrinos are produced [17].
The targeted neutrino signal is estimated from a dataset of simulated neutrino events reweighted to the energy and a e-mail: mortenmedici@gmail.com b Earthquake Research Institute, University of Tokyo, Bunkyo, Tokyo 113-0032, Japan.
directional distribution of dark matter in the Milky Way. The background is uniform in right ascension and is estimated from experimental data. A shape likelihood analysis on the reconstructed neutrino direction is used to estimate the fraction of events possibly originating from the targeted signal. From the signal fraction a limit on the signal flux is calculated and the corresponding value of σ A v can be determined for any combination of WIMP mass and WIMP annihilation channel to neutrinos. This search focuses on charged-current muon neutrinos because their directions can be accurately reconstructed. However, other neutrino flavors and events from neutralcurrent neutrino interaction are also present in the final selection (ensuring the most inclusive limits).

IceCube Neutrino Observatory
IceCube detects Cherenkov light from charged particles moving through one cubic kilometer of very transparent ice underneath the South Pole [18,19]. The array consists of 78 vertical strings in a hexagonal grid with 60 digital optical modules (DOMs) [20] spaced evenly on each string every 17 m between 1450 and 2450 m below the surface. The spacing between these nominal strings is approximately 125 m (as shown by the black dots in Fig. 1). In addition there are eight strings in the central area (red dots in Fig. 1) with the DOMs more densely spaced constituting the infill IceCube/ DeepCore [21].
The fiducial volume used in this work is defined by DOMs located 2140-2420 m below the surface situated on the most central strings (indicated with a solid blue region in Fig. 1). The rest of IceCube is used as a veto volume to reject incoming and through-going atmospheric muons.
The strings outside the DeepCore sub-detector volume (indicated with a blue line in Fig. 1) are only used in the initial filtering of triggered data, and are chosen to be shielded by three rows of DOMs from the edge of the array.

Signal expectation
For WIMPs self-annihilating to various Standard Model particles (leptons, quarks, or bosons), the decay chain of the particles will ultimately produce leptons and photons. Depending on the WIMP mass (m DM ) and annihilation channel, a number of neutrinos will be produced in the decay chain, propagate to Earth, and can be detected in neutrino observatories.
Using PYTHIA [22,23], a generic resonance with twice the WIMP mass is forced to decay through one of the particle pairs (annihilation channels) considered and the energy spectra of the resulting neutrinos are recorded for all three neu- Fig. 1 The horizontal position of the deployed strings in the IceCube coordinate system. The blue line shows the strings constituting the DeepCore subdetector, strings outside of this region are used in the initial event rejection. The fiducial volume used in the final analysis is indicated with the solid blue region consisting of both nominal and dense strings trino flavors. This work considers WIMPs with masses from 10 to 1000 GeV self-annihilating through either b-quarks (bb), W -bosons (W + W − ), muons (μ + μ − ), or taus (τ + τ − ) to neutrinos. Annihilation directly to neutrinos (νν) is also considered. In Fig. 2 the energy spectrum, d N/d E, of muon neutrinos from a pair of 100 GeV WIMPs is presented for the annihilation channels considered in this analysis. The energy spectrum is shown after applying long baseline oscillations (determined from parameters in [24]).
For the W + W − -channel only WIMP masses above the mass of the W boson are probed. The energy spectrum of the νν-channel is dominated by the line at m DM , which is modeled with a Gaussian distribution with a width of 5% of m DM . This width provides the possibility to use the same simulated dataset, while still being consistent with a line spectrum after smearing by the event reconstruction. For the signal from the νν-channel a flavor ratio produced at the source of (ν e : ν μ : ν τ ) = (1 : 1 : 1) is used (though the most conservative limits are found for a flavor ratio of (1 : 0 : 0) at source resulting in 10-15% weaker limits). The results will be presented with a 100% branching ratio for each annihilation channel considered.
The rate of WIMP self-annihilation seen in a given solid angle is determined from the integrated dark matter den-  [25], a spherical profile is assumed with one of two standard radial distributions: Navarro-Frenk-White (NFW) [13] and Burkert [14] with parameter values from [26]. The resulting rate of dark matter self-annihilations along the line of sight is strongly dependent on the assumed halo density, with the largest discrepancies near the center of the Milky way where the density is largest. Because of the large uncertainty on the model parameters the dark matter halo model constitutes the largest systematic uncertainty.
The resulting differential flux of signal neutrinos produced by WIMP self-annihilation in the dark matter halo of the Milky Way from a solid angle of the sky, , is given as where the 4π arises from a spherically symmetric annihilation, l is the line of sight through the dark matter halo with density profile ρ(r ) as a function of radius r , and the factor of 1/2 and the squared WIMP mass and halo density profile arise from the fact that two WIMPs are needed in order to annihilate. A sample of neutrino events of each flavor is generated with energies between 1 and 1000 GeV using GENIE [27] and weighted to the targeted flux of Eq. (1) according to their flavor, energy, and arrival direction for each combination of m DM , annihilation channel and dark matter halo density profile. This neutrino sample provides the distribution of the targeted signal that is used in the shape likelihood analysis to determine the fraction of possible signal events in the experimental data.

Background estimation
The background consists of neutrinos with other astrophysical origin, atmospheric neutrinos, and atmospheric muons. At the energies considered, the event sample is dominated by atmospheric neutrinos and muons produced in cosmic ray induced air showers. The cosmic ray flux is isotropic in right ascension, so the atmospheric background can be estimated from experimental data by randomizing the arrival times of each event. Since IceCube has a uniform exposure this corresponds to randomizing the right ascension values, which has shown in a previous analysis to be an unbiased approach to estimate the background [28].
The largest expected background contribution is from down-going atmospheric muons. This is because IceCube is located at the South Pole, so the center of the Milky Way (corresponding to the direction with the strongest signal) will be above the horizon, where there will also be the highest rate from atmospheric muons. Therefore the goal of the initial event selection is to reduce the rate of atmospheric muons. The overall analysis is verified using a simulation of atmospheric muons generated with CORSIKA [29] compared to the experimental data. The rate of simulated background is within 5% of the experimental data (see Table 1).
The other significant background contribution is atmospheric neutrinos. They arrive at IceCube from all directions and cannot be distinguished from extraterrestrial neutrinos event-by-event. However, from the full statistical ensemble the distributions can be distinguished by their energy and arrival direction. Simulated GENIE neutrino datasets are used for estimating the fraction of atmospheric neutrinos in the final selection of the experimental data, using the atmospheric neutrino flux model described in [30]. The simulated atmospheric neutrinos do not impact the result, as the combined background is estimated from experimental data.
The extra-galactic neutrino background can be distinguished from the WIMP neutrino signal by the arrival distri-bution, which is not necessarily the case for galactic neutrinos. But at the energies considered, both are expected to be more than three orders of magnitude below the background of atmospheric neutrinos.

Event selection
The event selection was optimized for the signal of muon neutrinos from 100 GeV WIMPs self-annihilating through the W + W − -channel (benchmark channel) and is applied event wise on the experimental data and the simulated event samples. The aim is to select high quality neutrino induced muons, signified by elongated event topologies (referred to as tracks) starting inside IceCube/DeepCore.
The neutrino induced muons need to be distinguished from the muons produced in the atmosphere. All atmospheric muons detected in IceCube penetrate through the veto volume. The corresponding hits (reconstructed pulses from one or more detected photons) can therefore be used to identify and remove these through-going tracks.
The event selection is a multi-step background rejection procedure that reduces the atmospheric muons by seven orders of magnitude.
The first step is to clean the DOM hits to remove noise so that the precision of the reconstruction is not degraded. Next, events with more than one hit in the volume outside the Deep-Core sub-detector volume causally connected to a charge weighted center of gravity in the fiducial volume within a predefined time window and distance are removed. This filters out atmospheric muons with very basic event information.
By requiring more than ten hits distributed on at least four strings nearly all noise-only events are removed. In addition, this requirement ensures that the events can be well reconstructed. The three first hits in the event are required to be in the fiducial volume, as that is more likely to indicate a starting event and thus reduce the rate of penetrating atmo- Table 1 Event rates for the various components expected in the experimental data given in mHz, and the signal neutrinos are presented as percentage of the events at filtered level for the benchmark signal (annihilation of a 100 GeV WIMP to W + W − ). Everything but the experi-mental data is based on simulation. The atmospheric muons rates are based on the GaisserH3a energy spectrum [34]. The atmospheric neutrinos rates are based on neutrino oscillation parameters in [35]. Due to vanishing rates at higher levels the rate of atmospheric ν τ are not listed spheric muons. The events are reconstructed to preliminarily estimate the direction and interaction point of the candidate neutrino-induced muon. The events with a preliminary zenith angle for the arrival direction of zen > zen GC + 20 • or zen < zen GC − 10 • are rejected, where zen GC denotes the zenith of the Galactic center. The cut is asymmetric because the atmospheric muon background is increasingly larger towards a zenith of zero (i.e. the southern celestial pole). A containment cut is used to keep only events that have a reconstructed interaction vertex within a cylinder with a radius corresponding to the analysis volume depicted on Fig. 1. In addition cuts are applied on track quality [31]. By considering the hits in the veto volume that are cleaned away (as possible noise), clusters are determined for hits that are within 250 m and 1000 ns from each other and are registered earlier than the first quantile of cleaned hits. These clusters are required to have fewer than three hits, as larger clusters are generally observed more often for penetrating atmospheric muons.
A cone with a 20 • opening angle aimed towards the arrival direction is used to check for hits in the uncleaned hit series within 1 µs of the interaction. At most one hit is allowed, since events starting within the fiducial volume should have zero hits within the cone, but one accidental noise hit is allowed. Due to the high rate of atmospheric muons versus possible signal neutrinos, there is a class of background muon events where sparse hits in the veto volume are removed during the hit cleaning. The uncleaned hits in a cylinder with a radius of 250 m pointed towards the arrival direction starting behind the interaction vertex, are used to calculate the likelihood value for the reconstructed track. A high likelihood value indicates that the track probably originated from a penetrating muon, for which the hits deposited in the veto volume are erroneously cleaned away.
At the energies considered in this analysis, the reconstruction must take into account both the hadronic cascade and the muon produced in a typical muon neutrino charged current interaction. With the experimental data event rate reduced by six orders of magnitude from 2 kHz to 3.7 mHz by the cuts described, a more specific event reconstruction can be run. This low energy specialized event reconstruction fits all relevant parameters (direction, interaction vertex, muon track length, and hadronic cascade energy) simultaneously and takes into account both DOMs that did and did not detect any light. In order to thoroughly sample the complex likelihood space of the full 8-dimensional parameter space the Bayesian sampling inference tool MultiNest [32] is used.
The final step of the event selection is a multivariate analysis using a Boosted Decision Tree (BDT) [33]. First of all, the direction and vertex information from the specialised event reconstruction are used along with the number of hits in a 10 degree opening angle veto cone, updated with the specialised event reconstruction. Further, the difference in likelihood in reconstructing the event with a finite track (expected from a neutrino induced starting muon) compared to an infinite track (expected for a through-going atmospheric muon) is used. An additional veto technique traces back in the direction of arrival from the interaction vertex to look for charge on DOMs that would identify the event as a through-going muon misidentified as a starting event. Both the number of hits and the total charge identified by the veto are used in the BDT.
The events are selected based on the BDT score, optimized for the best sensitivity to the benchmark signal of a 100 GeV WIMP annihilating through W + W − . The same cut value is used across multiple WIMP masses and annihilation channels.
The median resolution in azimuthal angle is presented in Fig. 3 as a function of true neutrino energy. Because the azimuthal angle maps directly to right ascension, it provides the dominating separation between signal and background. A comparison of three combinations of WIMP mass and annihilation channel is presented in Fig. 4, illustrating a better resolution for cases where the neutrino spectrum continues to higher energies.
The final event selection results in a data rate of 0.27 mHz, corresponding to a reduction by 7 orders of magnitude from the initial triggering of the data, while retaining 6% of the benchmark signal of muon neutrinos. No cuts have been incorporated to explicitly remove non-muon neutrino flavors. In the final event sample the non-muon neutrinos of the targeted signal are present with a combined rate comparable to that of muon neutrinos. Using the GENIE neutrino simulation weighted to the atmospheric flux model, it is estimated that atmospheric neutrinos constitute one quarter of the final experimental data. A summary of the event selection rates and signal efficiency is given in Table 1.
In Fig. 5 the effective area at the final level is presented for the individual neutrino flavors combining both neutraland charged-current neutrino interactions.

Analysis method
The final event sample is filled into 2D histograms with bins covering the range [0, 2π ] rad in right ascension (RA) and [−1, 1] rad in declination (Dec) using the reconstructed values from the specialised event reconstruction. The bin width is chosen to be 0.4 and 0.63 radians for RA and declination, respectively, based on the resolution of the event reconstruction. In order to ensure a consistent analysis the same bin width is chosen for the combination of WIMP mass and annihilation channel that exhibits the worst resolution. The 2D distributions constitute the probability density functions (PDFs) used in the shape likelihood analysis described below. The shape of the 2D distribution of experimental data produces the data PDF which is compared to the expectation from the weighted signal distributions (or signal PDF) and The experimental data scrambled in RA (assigned a random RA value for each event) consist of a component of scrambled background and potential signal (also scrambled): where μ ∈ [0, 1] parametrizes the fraction of signal in the total sample. From Eq. 2 the background PDF can be estimated from the experimental data (by subtracting the scrambled signal) under the hypothesis that the background is uniform in RA and hence invariant under scrambling.
The total fraction of events within a specific bin i ∈ [bin min , bin max ] is calculated as a function of the signal fraction as In Fig. 6 an example of the relevant PDFs is presented over the full range in right ascension for a single bin in declination (dec ∈ [−1/3, −2/3]) where the largest difference between signal and background is expected. Since the background is uniform in right ascension and the signal is peaked around the position of the center of the Milky Way, it is in right ascension that the difference between signal and background can be found. Figure 6 also illustrates the difference in the targeted signal between the NFW and Burkert models of the dark matter halo density profile.
With a 2D binned shape likelihood analysis, the data PDF is compared to the expectation from the background PDF and the signal PDF, for multiple combinations of WIMP mass, annihilation channel, and halo profile. This way the most probable signal fraction is determined from the experimental data. The likelihood is calculated by comparing the number of observed events in the individual bins n obs (i), assuming a Poisson uncertainty on the number of events expected, determined from the total number of events filled in the histogram n total obs and f (i|μ) calculated in Eq. 3. This results in the following formulation of the likelihood function L(μ): Using the likelihood analysis, the best estimate of the signal fraction can be found by minimizing − log L, and if it is consistent with zero the 90% confidence interval is determined applying the Feldman-Cousins approach [36] to estimate the upper limit on the signal fraction μ 90% . Using the simulated signal neutrinos the signal fraction can be related to σ A v . The expected limit on σ A v in the absence of signal is calculated from 10,000 pseudo experiments sampled from the background-only PDF, from which the median value of the resulting 90% upper limits is quoted as the sensitivity.

Systematic uncertainties
The statistical uncertainty due to the limited number of events in the simulated datasets is insignificant compared to the systematic uncertainties, as the simulation holds 20 times more events than in the experimental data, after cuts. However, all systematic uncertainties are effectively negligible compared to the astrophysical uncertainties associated with the parameters of the dark matter halo models.
The biggest systematic uncertainty arises from the modelling of the ice properties and the uncertainty on the optical efficiency of the DOMs, which increase with lower neutrino energies, and therefore for lower WIMP masses. The precision of the detector geometry and timing are so high that the associated systematic uncertainty is negligible and therefore not included in this study.
The effect of experimental systematic uncertainties on the final sensitivity is estimated using Monte Carlo simulations of neutrinos with uncertainty values varied by ±1σ from the values used in the baseline sets. Each of the datasets with variations is run through the event selection and analysis, providing a different value for the sensitivity on σ A v . The difference between the baseline and the variation will be quoted as the systematic uncertainty on σ A v , for each of the variations. The systematic uncertainties are dependent on the neutrino energy, and hence on the targeted WIMP mass. Since the background is estimated from experimental data, the variations are applied to the signal simulation only.
The optical properties of the ice in IceCube have been modelled and show an absorption and scattering length that vary with depth, generally becoming more clear in the deeper regions of IceCube. For the experimental data there will always be a discrepancy between the ice the photons are propagating through, and the ice [37] assumed in the reconstruction (as the complicated structure of the real ice can not be perfectly modeled). This is also the case in simulation, where the latest iteration of the ice model is used in the Monte Carlo event simulation, but because of its complexity, cannot currently be used for reconstruction. While estimating the impact of using a different ice model for event reconstruction than used in the photon propagation simulation, it additionally accounts for the fact that the ice model in simulation is different from that used in simulation. The effect is calculated using a variant Monte Carlo simulation with a different ice model used for the photon propagation (the same as used in the event reconstruction). This results in a 5-15% (depending on WIMP mass, 10% for the benchmark channel) improvement in sensitivity on σ A v , compared to the baseline simulation.
The ice in the drill hole columns has different optical properties from the bulk ice. The scattering length is greatly reduced due to the presence of impurities. One effect of this column is to increase the detection probability for downgoing photons. Since the DOMs are facing downwards, no down-going photons would be observed without scattering.
The column ice is treated as having a much shorter geometrical scattering length: 50 cm as a baseline [37], implemented in simulation as photons approach the DOMs. The uncertainty on the scattering length is covered by including variations of 30 and 100 cm. This variation results in a 25-30% reduction or 5-10% improvement of the sensitivity on σ A v respectively (depending on WIMP mass, 25 and 8% for the benchmark channel).
The photon detection efficiency of the DOMs (combining the effect of the quantum efficiency of the PMT, photon absorption by the cables in the ice, and other subdominant hardware elements) is determined to 10% accuracy. Increasing or decreasing the DOM efficiency in the simulation corresponds to a 5-40% (depending on WIMP mass, 15% for the benchmark channel) effect that symmetrically improves or reduces the sensitivity on σ A v .
The systematic uncertainties are considered to be independent and the ±variation that results in the largest uncertainty for each systematic uncertainty is added in quadrature to form the total systematic uncertainty. These are included in the final result by scaling up the limits with the total systematic uncertainty.
The dominant theoretical systematic uncertainty is related to fitted parameters of the dark matter halo profiles. Considering the 1σ variation on both parameters for the individual models result in a 150-200% uncertainty on the sensitivity on Since this effect is theory-dependent, and may change as dark matter halo models evolve, it is not included in the total systematic uncertainty. Instead, the results are presented for both dark matter halo models.

Results
After the final event selection, 22,632 events were observed in 1005 days of IceCube data. The data are presented in Fig. 6 illustrating that the data are compatible with the background-only hypothesis. Since no significant excess has been observed, an upper limit on σ A v is determined. Figure 7 shows the 90% confidence upper limits (solid black line) for the W + W − -annihilation channel for the two dark matter halo profiles. The colored bands represent the range of expected outcomes of this measurement with no signal present. The result is very near the median sensitivity, and thus compatible with the background-only hypothesis, which is the case across all annihilation channels. Tables 2 and 3 show the final upper limits on σ A v for all annihilation channels and WIMP masses considered in this analysis after accounting for the systematic uncertainties.
IceCube has previously searched for a neutrino signal from annihilating dark matter in the center of the Milky Way, using a combined event selection at low and high energies. The low energy selection observed an underfluctuation that resulted in an enhanced limit on σ A v , while the high energy selection gave access to higher energies. This analysis improves on the previous result at most of the energies considered. In order to compare this work to previous results, Fig. 8 shows the upper limits on σ A v for the τ + τ − annihilation channel and NFW halo profile of this work to previous results from IceCube and other indirect dark matter detection experiments. It can   [28,[38][39][40] and ANTARES [41]. Also shown are upper limits from gamma-ray searches from the dwarf galaxy Segue 1 (Seg1) by FermiLAT+MAGIC [42] and from the galactic center by H.E.S.S. [43]. The 'natural scale' refers to the value of σ A v that is needed for WIMPs to be a thermal relic [44] be seen that the analysis presented in this paper sets the best limits of a neutrino experiment on WIMP self-annihilation in the galactic center for WIMPs with masses between 10 and 100 GeV annihilating to τ + τ − .

Conclusions
This analysis demonstrates the continued improvements in dark matter searches with neutrinos, providing a valuable complement to the bounds from Cherenkov telescopes and gamma-ray satellites. A more inclusive event selection and the use of an improved event reconstruction algorithm have increased the sensitivity of IceCube to the signal of dark matter self-annihilation. However, no significant excess above the expected background has been observed in 3 years of Icecube/DeepCore data. Upper limits have been put on σ A v providing the leading limits on WIMPs with a mass between 10 and 100 GeV for a neutrino observatory.