All-flavour Search for Neutrinos from Dark Matter Annihilations in the Milky Way with IceCube/DeepCore

We present the first IceCube search for a signal of dark matter annihilations in the Milky Way using all-flavour neutrino-induced particle cascades. The analysis focuses on the DeepCore sub-detector of IceCube, and uses the surrounding IceCube strings as a veto region in order to select starting events in the DeepCore volume. We use 329 live-days of data from IceCube operating in its 86-string configuration during 2011-2012. No neutrino excess is found, the final result being compatible with the background-only hypothesis. From this null result, we derive upper limits on the velocity-averaged self-annihilation cross-section,<\sigma_A v>, for dark matter candidate masses ranging from 30 GeV up to 10 TeV, assuming both a cuspy and a flat-cored dark matter halo profile. For dark matter masses between 200 GeV and 10 TeV, the results improve on all previous IceCube results on<\sigma_A v>, reaching a level of 10^{-23} cm^3 s^-1, depending on the annihilation channel assumed, for a cusped NFW profile. The analysis demonstrates that all-flavour searches are competitive with muon channel searches despite the intrinsically worse angular resolution of cascades compared to muon tracks in IceCube.

Abstract We present the first IceCube search for a signal of dark matter annihilations in the Milky Way using all-flavour neutrino-induced particle cascades. The analysis focuses on the DeepCore sub-detector of IceCube, and uses the surrounding IceCube strings as a veto region in order to select starting events in the DeepCore volume. We use 329 live-days of data from IceCube operating in its 86-string configuration during [2011][2012]. No neutrino excess is found, the final result being compatible with the background-only hypothesis. From this null result, we derive upper limits on the velocity-averaged self-annihilation cross-section, σ A v , for dark matter candidate masses ranging from 30 GeV up to 10 TeV, assuming both a cuspy and a flat-cored dark matter halo profile. For dark matter masses between 200 GeV and 10 TeV, the results improve on all previous IceCube results on σ A v , reaching a level of 10 −23 cm 3 s −1 , depending on the annihilation channel assumed, for a cusped NFW profile. The analysis demonstrates that all-flavour searches are competitive with muon channel searches despite the intrinsically worse angular resolution of cascades compared to muon tracks in IceCube.

Introduction
There is strong evidence for extended halos of dark matter surrounding the visible component of galaxies. Independent indications of the existence of dark matter arise from gravitational effects at both galactic and galaxy-cluster scales, as well as from the growth of primordial density fluctuations which have left their imprint on the cosmic microwave background [1]. The nature of the dark matter is, however, still unknown. The most common assumption is that dark matter is composed of stable relic particles, whose present-day density is determined by freeze-out from thermal equilibrium as the universe expands and cools [2][3][4]. We focus here on a frequently considered candidate -a cosmologically stable massive particle having only weak interactions with baryonic matter, namely a Weakly Interacting Massive Particle (WIMP).
Within this particle dark matter paradigm, the Milky Way is expected to be embedded in a halo of WIMPs, which can annihilate and produce a flux of neutrinos detectable at Earth. The differential flux depends on the annihilation cross section of the WIMPs as where σ A v is the product of the self-annihilation cross section, σ A , and the WIMP velocity, v, averaged over the velocity distribution of WIMPS in the halo, which we assume to be spherical, m χ is the WIMP mass, dN ν /dE is the neutrino energy spectrum per annihilation and J a (ψ) is the integral of the squared of the dark matter density along the line of sight. Therefore, searches for the dark matter annihilation signal in the Galactic halo can probe the WIMP self-annihilation cross-section, given their spatial distribution. The expected signal is particularly sensitive to the adopted density profile of the dark matter halo, which determines the term J a (ψ) in Eq. (1), ψ being the angle between the direction to the Galactic Centre and the direction of observation [5,6]. The density profile of dark matter halos determined by numerical simulations of structure formation is still under debate [7][8][9][10][11][12]. To explicitly quantify the effect of the choice of the halo profile on the results of our analysis, we adopt two commonly used models: the Navarro-Frenk-White (NFW) cusped profile [9], and the Burkert cored profile [8,13]. We use the values for the parameters that characterize each profile from the Milky Way model presented in [14]. The difference between the two profiles is relevant only within the Solar circle, i.e., at radii less than 10 kpc.
In this paper we use data from the IceCube neutrino telescope to search for high energy neutrinos from the Galactic Centre and halo that may originate from dark matter annihilations. There have been several studies triggered by the observation of a electron and positron excess in the cosmic ray spectrum [15][16][17] which favour models in which WIMPs annihilate preferably to leptons [18][19][20][21][22][23][24]. We keep, though, the analysis agnostic in terms of the underlying specific particle physics model that could give rise to WIMP dark matter. In this sense it is a generic approach, and our results can be interpreted within any model that predicts a WIMP.
We use data collected in 329.1 live-days of detector operation between May 2011 and March 2012. The analysis focuses on identifying particle cascades produced by neutral or charged current neutrino interactions occurring inside the DeepCore sub-array of IceCube, being thus sensitive to all flavours. The analysis does not explicitly try to remove muon tracks from charged current ν μ interactions, but the event selection has been optimized to identify and select the more spherical light pattern produced in the detector by particle showers.

The IceCube neutrino observatory
The IceCube Neutrino Observatory [25] is a neutrino telescope located about one kilometer from the geographical South Pole and consisting of an in-ice array and a surface air shower array, IceTop [26]. The in-ice array utilizes one cubic kilometer of deep ultra-clear glacial ice as its detec- Fig. 1 Schematic overview of the IceCube string layout seen from above. Gray-filled markers indicate IceCube strings and black markers indicate the DeepCore strings with denser DOM spacing. All IceCube strings marked with a black border are included in the definition of the extended DeepCore volume used in the analysis tor medium. This volume is instrumented with 5160 Digital Optical Modules (DOMs) that register the Cherenkov photons emitted by the particles produced in neutrino interactions in the ice. The DOMs are distributed on 86 strings and are deployed between 1.5 km and 2.5 km below the surface. Out of the 86 strings, 78 are placed in a triangular grid of 125 m side, evenly spaced over the volume, and are referred to as IceCube strings. The remaining 8 strings are referred to as DeepCore strings. They are placed in between the central IceCube strings with a typical inter-string separation of 55 m. They have a denser DOM spacing and photomultiplier tubes with higher quantum efficiency. These strings, along with some of the surrounding IceCube strings, form the Deep-Core low-energy sub-array [27]. In the analysis described below, an extended definition of DeepCore was used, which includes one more layer of the surrounding IceCube strings, leaving a 3-string wide veto volume surrounding the fiducial volume used, see Fig. 1. While the original IceCube array has a neutrino energy threshold of about 100 GeV, the addition of the denser infill lowers the energy threshold to about 10 GeV.
The analysis presented in this paper uses a specific Deep-Core trigger, which requires that at least three hits are registered within 2.5 µs of each other in the nearest or next-tonearest neighboring DOMs in the DeepCore sub-array. When this condition is fulfilled, the trigger opens a ±6 µs readout window centered around the trigger time, where the full in-ice detector is read out. The average rate of this trigger is about 260 s −1 .

Signal and background simulations
In order to keep the analysis general we will assume that WIMPs annihilate with 100 % branching ratio into a few benchmark channels (bb, W + W − , νν, μ + μ − and τ + τ − ) and present results for these cases. Those channels effectively bracket the final particle spectra of realistic models with several final states. The neutrino spectra were calculated using PYTHIA [28] by producing a resonance at twice the mass under consideration and forcing it to decay to the desired channel. The program then takes care of the further hadronization and/or decays in the standard way. We ignore the possible WIMP spin in this approach, which can effect the final neutrino spectrum, mainly when considering annihilations through the W + W − channel [29]. We assume that the detected neutrinos have undergone full flavour mixing given the very long oscillation baseline from the source, so there are equal numbers of the three flavours. The expected angular distribution of signal events in the sky is obtained by reweighting the originally simulated isotropic distribution by J a (ψ).
There are two backgrounds to any search for neutrinos from the Galaxy: atmospheric neutrinos and atmospheric muons, both produced in cosmic-ray interactions in the atmosphere. To estimate the effect of these backgrounds on the analysis, a sample of atmospheric muons was generated with the CORSIKA package [30] and a sample of atmospheric neutrinos was simulated with GENIE [31] between 10 GeV and 200 GeV, and with NUGEN [32] from 200 GeV up to 10 9 GeV, adopting the spectrum in [33]. However, the analysis does not use background simulations to define the cuts, but instead relies on azimuth-scrambled data. This reduces the systematic uncertainties and automatically accounts for any unsimulated detector behavior. The background simulations were used to verify the overall validity of the analysis and the performance of the different cut levels. Since the majority of triggers in IceCube are due to atmospheric muons, the distributions of the variables used in the analysis must agree between data and the CORSIKA simulation at early selection levels, while at higher selection levels the data should show a significant fraction of atmospheric neutrinos. Atmospheric muons and particles resulting from neutrino interactions in or near the detector are propagated through the detector volume and their Cherenkov light emission simulated. Cherenkov photons are then propagated through the ice using the PPC package [34], and the response of the detector calculated. From this point onwards, simulations and data are treated identically through further filtering and data cleaning.

Data selection
The triggered data are first cleaned of potential noise hits that could effect the performance of the track and cascade reconstruction algorithms. Hits that lie outside a predetermined time window around the trigger time or which do not have another causally connected hit within a predefined radius, are removed. The data is then filtered by a fast algorithm that selects events starting in the DeepCore fiducial volume, in order to remove events triggered by through-going atmospheric muons. The IceCube strings surrounding DeepCore are used as a veto for incoming tracks. The algorithm selects events with the "amplitude-weighted" centre of gravity of all hits inside the DeepCore volume, 1 and no more than one hit in the surrounding IceCube strings causally connected with that point. This filter reduces the passing event rate by nearly a factor of 10.
The event sample is further reduced by requiring a minimum number of eight hits in the event distributed in at least four strings. This ensures that the remaining events can be well reconstructed. The events are then processed through a series of reconstructions aimed at determining their type (cascade or track), arrival direction and energy. In a first stage, two first-guess reconstructions are applied; fits for a track hypothesis and for a cascade hypothesis are performed in order to obtain a quick characterization of the events and perform a first event selection. These fits are based on the position and time of the hits in the detector, but do not include information about the optical properties of the ice, in order to speed up the computation. The track hypothesis performs a χ 2 fit of a straight line to the hit pattern of the event, returning a vertex and a velocity [35]. The cascade hypothesis is based on determining the amplitude-weighted centre of gravity of the hits in the event and its associated time. The algorithm calculates the three principal axes of the ellipsoid spanned by the spacial distribution of hits, and the longest principal axis is selected to determine the generic direction of the event. Since the specific incoming direction along the selected axis is ambiguous, the hit times are projected onto this axis, from latest to earliest, to characterize the time-development of the track so that it points towards where the incident particle originated. The tensor of inertia reconstruction is generally only suitable as a first guess of the direction for track-like events, since for cascade-like events the three principal axes of the ellipsoid will be close to equal in size. This property, however, can be used to discriminate between tracks and cascades. Additionally, a series of cuts based on variables derived from the geometrical distribution of hits, as well as from information 1 The amplitude-weighted centre of gravity of an event is defined as r COG = a i r i / a i , where a i and r i are the amplitude and position of the ith hit. The sum runs over all the hits in the event (after hit cleaning). from the first guess reconstructions, are applied. These cuts bring the experimental data rate down by a factor of about 3000 with respect to trigger level, while keeping about 50 % of the signal, depending on the WIMP mass and annihilation channel considered.
At this point three sophisticated likelihood-based reconstructions are applied on all the remaining events. The likelihood reconstructions aim at determining a set of parameters a = (x 0 , t 0 , ξ , E 0 ) given a set of measured data points d i (e.g. time and spatial coordinates of every hit in an event). Here x 0 is an arbitrary point along the track, t 0 is the event time at position x 0 , ξ is the direction of the incoming particle and E 0 is the deposited energy of the event. The reconstructions attempt to find the value of a that maximizes the likelihood function, which is based on the Probability Density Function (PDF) of measuring the data point d i given the set of parameters a. For a cascade reconstruction there are seven degrees of freedom, while an infinite track reconstruction has only six since the point x 0 can be chosen arbitrarily along the track. The first reconstruction is based on an infinite track hypothesis, fitting only direction, not energy. The second reconstruction uses a cascade hypothesis, and it fits for the vertex position, direction and energy of the cascade. These two reconstructions use an analytic approximation for the expected hit times in the DOMs given a track or cascade hypothesis [36], rather than a full description of the optical properties of the ice. Since the focus of this analysis is to identify cascades, an additional, more advanced, cascade reconstruction is also performed, using the previous one as a seed. This second cascade reconstruction uses the full description of the optical properties of the Antarctic ice, as well as information of the position of non-hit DOMs through a term added to the energy likelihood. The three likelihood reconstructions return the best fit values of the variables of the vector a they fit for, as well as a likelihood value of their respective hypothesis, which is used in a further selection of events using linear cuts.
The final selection of events uses Boosted Decision Trees (BDT) [37] to classify events as signal-like or backgroundlike. Two BDTs were trained using data as background and a different benchmark reference signal each. One of the BDTs (BDT LE ) was trained using the neutrino spectrum from a 100 GeV WIMP annihilating into bb, while the other, BDT HE , was trained on the neutrino spectrum of a 300 GeV WIMP annihilating into W + W − . These two spectra were chosen to represent a soft and hard neutrino spectrum respectively, so the sensitivity of the analysis to other WIMP masses and/or annihilation channels with similar spectra can be evaluated with the same cuts on the BDT output scores. This removes the need to train a different BDT specifically for each mass and annihilation channel. Since no variables depending on the arrival direction of the events are used in the BDT Fig. 2 The score distribution for the BDT trained on the bb 100 GeV signal channel (left) and for the BDT trained on the W + W − 300 GeV signal channel (right). The plot shows the passing event rate (in Hz) for simulated atmospheric muons (blue line) and atmospheric neutrinos (green lines), as well as for the sum of these two components (total MC, purple line), compared with the data passing rate. The passing rate for each signal channel the BDT was trained for is also shown (grey lines), normalized to the experimental data rate. The final cuts on the score are marked with vertical lines. Events were kept if any of the scores were above the cut values. The lower panel in each plot shows the ratio of the data passing rate to the total expected background training, the event sample is kept blind with respect to the position of the Galactic Centre in the sky.
Seven variables that showed a good separation power between signal and background, selected among an initial larger set of variables that were tried, were used to train the BDTs. The variables are based on the different geometrical patterns that tracks and cascades leave in the detector, as well as on their different time development. The whole data set was classified by the two BDTs so each event was assigned two BDT scores. In order to decide on the best cut value on each BDT output, the range of BDT score values was scanned and the sensitivity of the analysis was calculated for each of them. The scores producing the best sensitivity for each of the two signal channels for which the BDTs were trained were selected. Events with a BDT LE score above the optimal value are referred to as the "low-energy" (LE) sample, and events with a BDT HE score above the corresponding cut value are referred to as the "high-energy" (HE) sample. The remaining number of events in each sample is 5892 events in the LE sample and 2178 events in the HE sample. The overlap between the two samples (events which have both BDT scores above the respective cut values) is 664 events. The final BDT score distributions for the 100 GeV bb and the 300 GeV W + W channels are presented in Fig. 2, with the vertical lines marking the optimal cut values used to select the final event sample.
After the BDT classification, the data has been reduced by a factor of about 1(3) × 10 6 for the LE(HE) sample, but still contains about 20 % of atmospheric muon contamination. The remaining signal in the two benchmark scenarios considered amounts to about 6 % (8 %) respectively. A summary of the event selection rates, as well as signal efficiency, is given in Table 1. The effective area for the two event selections, a measure of how efficient the detector is for the present Table 1 Data rates at different cut levels. For the two benchmark signal channels, 100 GeV bb and 300 GeV W + W − , values are given as percentage of signal retention relative to the DeepCore and event-quality filter level. The livetime of the experimental data set is 329.1 days   analysis, is shown in the left plot of Fig. 3. The right plot in the same figure shows the cumulative angular resolution (space angle between the reconstructed and true direction of the incoming neutrino) for the two benchmark channels used in training the BDTs.

Systematic uncertainties
In order to estimate the effect of experimental systematic uncertainties on the final sensitivity, Monte Carlo simulation studies were done, where the parameters defining a given input were varied within their estimated uncertainty. The main source of systematic uncertainties is the limited knowledge of the optical properties of the ice, both the bulk ice between 1450 m and 2500 m, as well as the "hole ice", i.e. the ice that forms as the water in the hole drilled for the string deployment refreezes. The scattering and absorption coefficients of the ice as a function of depth have been determined by in-situ flash measurements, and a standard "ice model" for IceCube has been derived [38]. The effect on the uncertainty of the estimated absorption and scattering length was investigated by varying the baseline settings by ±10 % individually. Their contribution to the uncertainty on the sensitivity lies in the range 8 %-12 %. Furthermore, there are indications that the hole ice contains residual air bubbles that result in a shorter scattering length in this ice compared to the ancient glacial bulk ice surrounding it. In the baseline simulation data sets the scattering length of the hole ice is set to 50 cm. Varying this parameter between 30 cm and 100 cm yields a 10 %-24 % change on the sensitivity. Recently, a more detailed modeling of the bulk ice has been developed [39]. It includes anisotropic scattering and accounts for the tilt of the different ice layers across the IceCube volume. Preliminary studies indicate that the effect on the sensitivity of this model is negligible for high-energy events, but it can be sizable for the lowest-energy events, reducing the sensitivity for low WIMP masses up to 25 %. These effects have not been included in this analysis. The overall efficiency of the process of converting the Cherenkov light into a detectable electrical signal by the DOM is another source of uncertainty. This effect was investigated by changing the DOM efficiency in the signal simulation by ±10 %, according to measurements of the performance of the DOMs in laboratory tests before deployment, as well as in in-situ calibration measurements after deployment. This uncertainty translates into an uncertainty on the final sensitivity of 10 %-35 %, depending on event selection. The effect is stronger for low-energy events that can fall under the detector threshold if less light is being captured. Additional, but minor, effects arise from the implementation of the photomultiplier dark noise in the simulation, the timing and geometry calibration of the detector and from the intrinsic randomness of several steps of the analysis, like time-scrambling of the data or the many pseudo-experiments performed to calculate the sensitivity.
All systematic uncertainties considered are summarized in Table 2 together with the total (quadratic sum) for the low and high-energy selections for both halo profiles. In order to be conservative, the limits presented in Sect. 6 for each WIMP mass and annihilation channel were rescaled by the corresponding total systematic uncertainty shown in Table 2.

Analysis method
We use the distribution of the space angle ψ between event directions and the Galactic Centre to construct a likelihood function and test the signal hypothesis (excess of events at small ψ values) against the background-only hypothesis (an event distribution isotropic in the sky). The signal and background hypotheses are represented by probability density functions of the ψ distributions, where the subscripts S and B denote signal and background respectively and μ is the number of signal events present among the total number of observed events, n obs . The angle ψ is allowed to be in the full range [0 • , 180 • ], therefore covering the full sky, as shown in Fig. 4. This allows the analysis to be sensitive to the whole halo instead of just to the Galactic Centre. However, if the signal is allowed to come from anywhere in the halo, the background distribution, taken from data, is necessarily contaminated by a potential signal: thereby the dependence of f B (ψ | μ) on μ and not only on ψ. In particular the background distribution is constructed as where f ss and f sd are the PDF of the scrambled arrival directions of signal simulation and data events respectively. The likelihood that the data sample contains μ signal events is defined as where n obs is the number of observed events and f (ψ i | μ) is given in Eq. (2). We follow the method described in [40] to calculate a 90 % confidence level upper limit on μ, μ 90 , which gives an upper limit on the flux of neutrinos from the halo as defined in Eq. (1). This limit can, in turn, be translated into a limit on σ A v for any given WIMP mass, annihilation channel and halo profile.
The final limits are shown in the next section, for the event selection that showed the best sensitivity in each case.

Results and conclusion
At final selection level, a total of 5892 (2178) events were observed in the full sky for the low-energy (high-energy) samples respectively. Figure 5 shows the angular distribution of the two event samples at final cut level. The distributions are compatible with 0 signal events for all WIMP masses and annihilation channels tested. Tables 3, 4, 5 and 6 show the results for the best fit on the number of signal events,μ, together with the 90 % upper limits on the number of signal events, μ 90 , and the corresponding limit on the thermally-averaged WIMP annihilation cross section, σ A v 90 . Corresponding quantities with a tilde denote median upper limits (i.e., sensitivities). Each table corresponds to a given benchmark annihilation channel and it shows different WIMP masses for the two halo models considered. The available statistics at final level in the case of direct annihilation of 700 GeV WIMPs to neutrinos using the Burkert profile were not sufficient to define an angular distribution which was smooth enough to perform the shape analysis, so we choose not to quote results for this mass and channel in Table 6. Figures 6 and 7 show the results graphically for the NFW and Burkert dark matter profiles respectively. The plots show the 90 % C.L. upper limits (solid black line) on the velocity-averaged WIMP self-annihilation cross section, σ A v , together with the corresponding sensitivities (dashed black line) and the 1σ (green) and 2σ (yellow) statistical uncertainties.
In order to put the results of this analysis in perspective, Fig. 8 shows a comparison with results from previous Ice-Cube analyses and other experiments, for the τ τ annihilation channel and the NFW profile. Also shown is the allowed area in the ( σ A v , m χ ) parameter space if the e + + e − flux excess seen by Fermi-LAT and H.E.S.S. and the positron excess seen by PAMELA are interpreted as originating from dark matter annihilations [41]. There exist, however, conventional explanations based on local astrophysical sources [42,43] that, along with current limits on σ A v , disfavour such explanation. The figure shows that the analysis presented in this paper improves on previous IceCube analyses [44][45][46][47] for WIMP masses above about 200 GeV, as well as on the ANTARES [48] result for WIMP masses below ∼1 TeV. This demonstrates that particle cascades can be reconstructed with a good enough angular resolution in IceCube to make this channel competitive in searches for dark matter signals with neutrinos from the Galactic Centre and halo. Even if Cherenkov telescopes and gamma-ray satellites can reach stricter bounds on σ A v due to their better angular resolution and, depending on the source under consideration, low background, there is a much-needed complementarity in the field of dark matter searches, where neutrino telescopes can play a valuable role.   [44][45][46][47]. Also shown are the latest upper limits from gamma-ray searches obtained from the combination of FermiLAT and MAGIC results [49]. The three shaded areas indicate allowed regions if the e + + e − flux excess seen by Fermi-LAT, H.E.S.S. and the positron excess seen by PAMELA (3σ in dark green, 5σ in light green and gray area, respectively) would be interpreted as originating from dark-matter annihilations. The data for the shaded regions are taken from [41]. The natural scale denotes the required value of σ A v for a thermal-relic to constitute the dark matter [50]