Improved limits on dark matter annihilation in the Sun with the 79-string IceCube detector and implications for supersymmetry

We present an improved event-level likelihood formalism for including neutrino telescope data in global fits to new physics. We derive limits on spin-dependent dark matter-proton scattering by employing the new formalism in a re-analysis of data from the 79-string IceCube search for dark matter annihilation in the Sun, including explicit energy information for each event. The new analysis excludes a number of models in the weak-scale minimal supersymmetric standard model (MSSM) for the first time. This work is accompanied by the public release of the 79-string IceCube data, as well as an associated computer code for applying the new likelihood to arbitrary dark matter models.


Introduction
Searches for high-energy neutrinos from the Sun are currently the most sensitive means of probing spin-dependent interactions between protons and most models for dark matter (DM) [1,2]. Most analyses take a semi-model-independent approach, assuming that capture and annihilation have reached equilibrium in the Sun, and that DM annihilates exclusively into a single final state. These assumptions are expressly violated in many concrete models for the identity of DM, including supersymmetry [3][4][5][6][7]. Resulting limits are often difficult to meaningfully connect to theoretical predictions [8][9][10][11][12][13][14][15], in part because the necessary data and likelihood functions for recasting limits to other theories are unavailable. The computational expense required to replicate the experimental analyses for millions of parameter combinations can also be prohibitive. All these issues arise in some form in direct detection, collider searches and other forms of indirect detection as well [16][17][18][19][20][21]. This paper provides a solution to these problems for the indirect dark matter search with neutrinos. We previously presented a 79-string search for dark matter annihilation in the Sun (IC79; [1]), deriving limits on single annihilation channels. We later developed a formalism (Paper I; [7]) that allows event-level neutrino telescope data to be used to constrain DM models with mixed annihilation final states, thereby allowing IceCube searches to be properly included in global fits to theories beyond the Standard Model of particle physics (BSM). Paper I provided methods applicable to neutrinos with high energies (50 GeV and above) that were observed with the 22-string configuration of IceCube. This paper (Paper II) revises this formalism to include the impact of non-negligible angles between the neutrino direction and the muon produced, extending the reach of the technique to neutrino energies as low as 10 GeV. We then apply the formalism to IC79 data and use it to rule out some example supersymmetric models. Compared to the original IC79 analysis [1,22], which was based solely on the observed arrival directions of events, here we also include event-level energy information and an explicit treatment of the total number of observed events within the signal region, leading to an improvement in limits at high DM masses. Extensive references on neutrino searches for dark matter and BSM global fits can be found in Paper I.
We publicly provide the fast likelihood code (nulike 1 ) that implements the improved analysis presented in this paper, using the public IC79 event information and detector response. Nulike also provides pre-computed, fully model-independent 'partial likelihoods' for every event observed by IC79, making new limits quick and easy to obtain for any annihilation final state or combination thereof. This is a distinct advantage over the standard IceCube analysis pathway, where full signal propagation and detector simulations are required for each model. While the approach in this paper relies on many results of the direct simulation method, such as effective areas and volumes, it provides a complete framework in which they can then be applied to essentially any neutrino annihilation signal that can be safely treated as a point source. The methods and the corresponding code are agnostic with respect to the details of the experiment and can be used to perform similar analyses for other neutrino telescopes, given appropriate input data in the form of event and detector response files.
In Section 2, we will provide details of the IC79 data that we use in the updated analysis, before describing the improved likelihood formalism in Section 3. We then show the impacts of the new analysis on generic weakly-interacting massive particle (WIMP) models in Section 4 and models in the minimal supersymmetric standard model (MSSM) in Section 5. We will conclude in Section 6.
2 The 79-string IceCube search for dark matter

The IceCube detector
Completed in December 2010, the IceCube neutrino observatory [23] is a neutrino telescope situated at the South Pole. IceCube is installed in the glacial ice at depths of between 1450 m and 2450 m, instrumenting a total volume of one cubic kilometre. Digital Optical Modules (DOMs) arranged on vertical strings deep in the ice sheet record the Cherenkov light induced by relativistic charged particles, including those created by neutrinos interacting with the ice. The detection of photon yields and arrival times in DOMs allows for the reconstruction of the directions and energies of the secondaries. In its 79-string configuration, 73 strings have a horizontal spacing of 125 m and a vertical spacing of 17 m between DOMs. The six remaining strings are located near the central string of IceCube and feature a reduced vertical spacing between DOMs of 7 m and higher quantum efficiency photomultiplier tubes. Along with the seven surrounding regular strings, they form the DeepCore subarray [24]. The horizontal distance between strings in DeepCore is less than 75 m. The higher sensor density in clear ice provides an order of magnitude lower energy threshold of O(10) GeV compared to the main IceCube array.

Data samples
In the analysis described in this paper, we start with pre-selected data from a search for WIMP dark matter annihilation in the Sun with the IceCube 79-string configuration [1]. This analysis uses 317 live-days of data taken between May 2010 and May 2011. As described in Refs. [1,22], the DeepCore subarray is included for the first time in the analysis, lowering the energy threshold and extending the search to the austral summer (when neutrinos from the Sun pass downwards through the ice). In order to be sensitive to a wide range of potential WIMP masses, the analysis comprises three independent non-overlapping event selections. First, the full dataset is split into two seasonal streams, where September 22nd 2010 and March 22nd 2011 mark the beginning and end of the 'summer' dataset. The 'summer' sample ('summer low-energy' event selection, SL) is a dedicated low energy event sample that uses the surrounding IceCube strings as an instrumented muon veto in order to select neutrinoinduced events that start within DeepCore. The 'winter' dataset comprises two samples. The first sample ('winter high-energy' event selection, WH) has no particular track-containment requirement and aims to select upward-going muon tracks. The second sample ('winter lowenergy' event selection, WL) is a low energy sample, and focuses on neutrino-induced muon tracks that start or are fully contained in DeepCore. The event selection was carried out separately for each independent sample. By design, the uncorrelated nature of the three datasets makes it straightforward to combine them in a joint likelihood. The analysis in sections 4 and 5 uses the event-level data at final analysis level and corresponding signal simulations from [1] and [22].

Signal and background simulation
Solar WIMP signals are simulated using WIMPSim [25], which describes the annihilation of WIMPs inside the Sun. WIMPSim simulates the production, interaction, oscillation and propagation of all three flavours of neutrinos from the core of the Sun to the detector.
Muons arising in single or coincident air showers as well as atmospheric neutrinos form the background to this analysis. We did not simulate these contributions, as they can be estimated by scrambling real data at the final analysis level (detailed within section 2.7).

Calculation of detector efficiencies
The effective volume V eff (E µ ) of the detector for muon or anti-muon events produced through charged current interactions differs for each of the three event selections of Ref. [1]. V eff (E µ ) for the detection of muons from the Sun is a function of muon energy, averaged over the live-time of the respective event selections. It corresponds to an equivalent volume of 100% detection efficiency, and is identical for both muons and anti-muons. We also calculated the effective area A eff (E) for detection of muon neutrinos as a function of neutrino energy. We use A eff (E) later to compute 'bias factors', which account for selection effects in the analysis (see Section 3.6). The effective areas for muon neutrinos and muon anti-neutrinos differ due to the differences in the (anti-)neutrino cross-sections with hadrons. All effective volumes and areas for the 79-string analysis are available online [26].
We specify the total systematic uncertainties related to the detector response at the 1 σ confidence level within each energy bin, in a manner similar to how it was done in Paper I. These uncertainties come from simulation studies, where identified sources of uncertainty, e.g. absolute DOM efficiency, photon propagation in ice, or calibration constants, were individually varied within reasonable ranges of their original values. Similarly, the uncertainties arising from limited simulation statistics are also given for each energy bin of V eff , at the 1 σ confidence level. In the final analysis we combine these two errors in quadrature.

Angular response
The point spread function (PSF) describes the uncertainty in the reconstructed arrival direction of muons. Closely following Paper I, the reduced (one-dimensional) PSF for the angular deviation ∆ between the true arrival direction of a muon on the sky and its reconstructed direction is We extract the parameter σ µ , which we refer to as the 'mean angular error', directly from the one-dimensional PSF constructed from IceCube signal simulations. As in Paper I, we determine σ µ in the same energy bins that were used for calculating the detector efficiencies. For simplicity, we neglect the curvature of the PSF on the sky, owing to the fact that for dark matter signals detected with DeepCore, the muon production angle is typically expected to be the dominant source of angular deviation. We therefore restrict our analyses to signal regions of radii φ cut around the solar position on the sky so as to minimise the error induced by this approximation (and the fact that we include the entire sky in a data-driven estimation of the background; cf. Sec. 2.7). We determined that φ cut = 20 • provides satisfactory signal acceptance and background rejection for the WH sample, and φ cut = 40 • is appropriate for the WL and SL datasets. We associate angular uncertainties with real data events on an event-by-event basis, using the paraboloid method [27]. A paraboloid function is fitted to the muon track reconstruction likelihood space in the neighbourhood of the best fit. The resulting confidence ellipse on the sky is represented by the two principal axes, which correspond to the standard deviations of the likelihood function in each of two linearly-independent directions. The overall reconstructed likelihood track uncertainty, σ para (the 'paraboloid sigma'), is calculated as the mean in quadrature of the uncertainties along the two axes. Good track fits generally result in paraboloids that are narrow along both axes and therefore have small σ para values.

Energy estimator
Paper I used the number of lit DOMs (N chan ) as a suitable energy estimator. This definition worked well for a detector with a consistent density of optical modules, like the 22-string configuration of IceCube. This paper uses data recorded in the 79-string configuration of IceCube. This configuration includes the DeepCore subarray, which has a higher density of DOMs than the rest of the detector (Sec. 2.1). A simple count of lit DOMs would yield different results depending on whether the event crosses, partially crosses, or is contained within DeepCore. In an attempt to address this, we introduce a modified N chan value, N c chan , which corrects for the variation in DOM density across the detector. In this context, the corrected energy proxy N c chan is where N DC chan and N IC chan are the number of lit DOMs in DeepCore ('standard' IceCube and 'high quantum efficiency' DOMs) and the remainder of IceCube, respectively. The factor f DC = 0.28 is the the ratio of the number of 'standard' IceCube DOMs inside DeepCore to the total number of DOMs in DeepCore, multiplied by an additional correction factor.  Figure 1: Predicted probability distributions of N c chan for the WH event selection, derived from high-statistics simulations used in [1]. Each distribution is defined for muons having energies in a specific logarithmic energy interval of width 0.2. The fitted functions are to guide the eye only and are not used in our calculations. The lower plot compares the fitted functions, illustrating the ability to differentiate events between different energy intervals.
The ratio accounts for the higher density of DOMs in DeepCore compared to the rest of the detector, and the additional correction factor accounts for the higher quantum efficiency of  the photomultiplier tubes in the 'high quantum efficiency' DOMs.
We calculated the expected distributions of observed N c chan values for a series of intervals in muon energy, as we did in neutrino energy in Paper I. Figs. 1 -3 show these probability distributions for each event selection (WH, WL and SL) and muon energy range. The total interval in muon energy is different for each event selection due to the respective event selection criteria that are applied. We use these probability distributions, together with the  predicted energy spectrum of the signal from each WIMP model, to calculate the predicted distribution of N c chan . The fitted functions in Figs. 1 -3 are only to guide the eye; our signal predictions and likelihood calculations employ the actual distributions. The lower plots in Figs. 1 -3 compare the fitted functions for each event selection, illustrating the ability to differentiate events between different energy intervals. To reach energies as low as the first interval in Figs. 2 and 3, DeepCore uses an independent, low-threshold, simple majority trigger (SMT), with a 2.5 µs time window, applied to DOMs comprising the DeepCore fiducial volume [24]. This trigger requires that three or more DOMs satisfy the so-called hard local coincidence (HLC) condition (as opposed to the threshold of eight or more DOMs, more typically used in IceCube analyses). DOMs meet the HLC condition when two or more DOMs in close proximity to each other (nearest or next-to-nearest neighbours on the same string) register hits within a 1 µs time window. This trigger is 70% efficient for a simulated sample of atmospheric ν µ events of 10 GeV neutrino energy [24].
More advanced energy reconstruction methods other than N c chan are available in IceCube that are based on the reconstruction of charged-particle energies and topologies from the observed Cherenkov light yield [28]. Here we use N c chan for simplicity and robustness.

Background estimation
As in Paper I, the background distributions for each event selection come directly from data. The angular distribution of background events dP BG (φ )/dφ is a function of φ , the angle between the reconstructed track direction and the Sun. Muons produced in cosmic-ray showers are the dominant contributors to the background. Their angular distribution is observed to be largely independent of azimuth, so we estimated dP BG (φ )/dφ from real data events at the final selection level with scrambled azimuths. We used all observed events at the final selection level for this exercise. Given the tight upper limit on a signal contribution in the original analysis [1], including the nominal signal region does not bias the background estimate. We calculated the distribution of N c chan due to background events, dP BG (N c chan )/dN c chan , and observed no significant correlation between the arrival angles of events relative to the Sun and their measured N c chan values.

Data format, public code and availability
Full event data from the analysis of Ref. [1], including angles, N c chan values and paraboloid sigmas, can be found at http://icecube.wisc.edu/science/data/IC79 solarWIMP data release. Effective areas and volumes, along with N c chan and angular responses, can be found at the same location.
The nulike code can be downloaded from http://nulike.hepforge.org. The release of nulike coincides with the release of DarkSUSY v5.1.3. This release of DarkSUSY provides optimised interpolation routines for WIMPSim [25] outputs contained in DarkSUSY, and ensures that they are fully compatible with the parallel likelihood routines in nulike (i.e. the routines one would use together with nulike are threadsafe in the latest DarkSUSY release).

General form
The primary improvement in the likelihood treatment here compared to Paper I [7] is that we allow for differences between the arrival directions of neutrinos (φ) and the muons they produce (φ µ ). At neutrino energies above O(100 GeV), to a good approximation one can neglect the difference between φ and φ µ . This was the case for all data and calculations considered in Paper I. With the DeepCore infill array however, the actual 79-and 86-string IceCube configurations are sensitive to neutrino energies even below 10 GeV. For example, for a neutrino of energy 10 GeV producing a muon of 4 GeV, φ − φ µ can be as large as 30 degrees, and must therefore be explicitly included in all calculations.
The distribution of muon production angles introduces an explicit energy dependence to the detector PSF. This improves on our earlier approximation that the detector response factorises into separate functions of angle and energy (Eq. 3.6 in Paper I). In this paper we therefore work with the general form of the unbinned likelihood, The vector ξ ξ ξ refers to the parameters of a given BSM model. N c i and φ i are the actual observed event-level data for the ith event of n tot total events. N c i in this analysis is the generalised N chan , whereas φ i is the angle between the reconstructed muon track and the direction of the Sun. As in Paper I, is the probability density (in effective units of inverse angle and N c chan ) for observing N c i and φ i for the ith event when the true values of the incoming neutrino energy and angle relative to the Sun are E and φ, respectively.
The prefactor L num is the number likelihood for observing n tot events given a prediction θ tot , marginalised over the systematic error on the predicted number of events. This is where θ S is the predicted number of signal events, θ BG is the predicted number of background events, θ tot = θ S + θ BG , is the rescaling variable assumed to have a log-normal distribution, and σ is the fractional systematic error on the signal prediction (which sets the width of the distribution of ). The width σ is the sum in quadrature of a theoretical error τ and the fractional uncertainty on the detector response. This treatment requires the selection of a single indicative systematic error on the effective volume, which is then applied identically at all muon energies. When computing results, to be conservative we chose the largest systematic error on the effective volume over the entire range of detectable muon energies. For the theoretical error τ we adopted a minimum of 5% for WIMP masses m χ ≤ 100 GeV to account for neglected higher order corrections and round-off errors, increasing to 50% at m χ = 10 TeV as τ = 0.05 m χ 100 GeV This sliding scale is designed to encapsulate the increasing error with WIMP mass of predicted spectra from DarkSUSY, due to internal tables in which it interpolates results from WIMPSim. Paper I and Refs. [16,29] give further details and background on the number likelihood. The expected distribution of incident neutrino energies (E) and angles (φ) is given by d 2 P/dE dφ(E, φ, ξ ξ ξ), which is a prediction of the model parameters ξ ξ ξ. This separates into a weighted sum of the signal (S) and background (BG) contributions, so that Eq. 3.1 can be expressed as where f S ≡ θ S /θ tot and f BG ≡ θ BG /θ tot are the fractions of the total expected events from signal and background, respectively, and gives the signal (X = S) and background likelihoods (X = BG).

Background likelihood
The calculation of the background likelihood component follows the treatment in Paper I closely: the integral in Eq. 3.5 for X = BG is the actual observed background, which is independent of the model parameters ξ ξ ξ. Within the zenith angle range considered in this analysis, to a very good approximation, the background spectrum and angular distributions are not correlated. L BG,i can then be written as where dP BG /dN c i and dP BG /dφ i are the observed N c chan and angular distributions of the background, respectively (Sec. 2.7). The expected number of background events θ BG used to calculate the background fraction f BG refers to the events contained in the angular cut φ cut around the solar position.

Signal likelihood
In order to take into account the distribution of production angles in calculating the signal likelihood, the integral in Eq. 3.5 for X = S should be expressed in terms of the kinematics of the produced muons. In Eq. 3.5 this integrand is the product of the predicted arrival probability of a neutrino of a given energy and arrival angle, and the detector response to it. We express this as the product of the predicted differential flux of incoming neutrinos (d 2 Φ ν /dE dφ), the exposure time of the observation (t exp ), the effective differential crosssection for neutrino conversion into muons in the ice (d 2 Σ ν→µ /dE µ dφ µ ), and the response of the detector to muon-conversion events (Q µ ). We then integrate over the distribution of muon energies and angles that might be created in the interaction, so as to recover a pure function of the neutrino properties (as the theoretical predictions of different dark matter models ξ ξ ξ are given at neutrino level). We divide by the expected number of signal events θ S inside the angular cut cone, in order to normalise the integral of the resulting probability distribution to unity. We also multiply by a bias factor f b (E), which is an analysis-dependent function of the neutrino energy. θ S and f b (E) are discussed in detail in Sec. 3.5. Finally, we add the contributions of both incoming neutrinos and antineutrinos, giving: Here φ µ is the angle of the produced muon relative to the Sun, E µ is its energy, and barred quantities are the equivalent measures for anti-particles. The angular component of the signal prediction is a delta function at the solar position, so the integral of Eq. 3.7 over φ (required by Eq. 3.5 in order to obtain the signal likelihood) can be done analytically. We then find where the sum indicates that the corresponding antiparticle expression must also be included.
With φ = 0, the true muon arrival angle relative to the Sun, φ µ , becomes identical to the microscopic muon production angle in the frame where the target nucleus is at rest. The value of this angle depends on the incoming neutrino energy and the outgoing muon energy, as well as the momentum carried by the parton within the nucleon with which the neutrino interacts. It can be written as where m N refers to the mass of the nucleon involved. The Björken scaling variable x indicates the fraction of the nucleonic momentum carried by the parton involved in the interaction. By definition, x varies between 0 and 1, as does the other Björken variable y = 1 − E µ /E.
Together, x and y provide a convenient and well-bounded way to express the dependence of the neutrino interaction cross-sections on the outgoing muon energy and angle. We therefore trade E µ and φ µ for x and y, so that remembering that E µ = E µ (y, E) and φ µ = φ µ (x, y, E). For each observed event inside our analysis cone, we precompute the inner double integral of Eq. 3.11 for a set of 50 logarithmically-spaced neutrino energies per decade over the range 0.5 ≤ log 10 (E/GeV) ≤ 4.0. To obtain the contribution of the predicted signal to the total likelihood for that event, we re-weight these 'partial likelihoods' according to the predicted neutrino spectrum for each model ξ ξ ξ, as well as the bias factor f b . To allow for a fast and straightforward application to any theoretical neutrino spectrum, we provide the partial likelihoods for the 79-string IceCube analysis in nulike, precomputed, along with routines for computing the bias factors f b . We also provide the underlying event data online [26] and a utility within nulike that can precompute and save the partial likelihoods from any other neutrino telescope, provided the data are in the same format. The effective differential conversion cross-section is given by where it should again be understood that E µ = E µ (y, E). The replacement {ν, µ} → {ν,μ} provides the corresponding expression for d 2 Σν →μ /dx dy. This is the product of the number density n N of nucleon species N (proton or neutron) in the detector, the effective volume V eff (E µ ) of the detector for muon or anti-muon conversion events, and d 2 σ ν→µ,N /dx dy, the microscopic differential cross-section for muon production by charged-current interactions.
V eff is the same for both muons and anti-muons. In contrast, the differential cross-sections differ for particles and antiparticles. These are known from the theory of weak interactions, up to a dependence on the parton distributions for x. We obtain these from nusigma [30], which by default relies on the CTEQ6-DIS parton distribution functions [31]. Users of nulike who prefer other parton distributions can simply switch those employed by nusigma and recompute the partial likelihoods.

Detector response
Based on the observation that the angular and spectral (N c chan ) distributions of detected events are essentially uncorrelated across the sky (Sec. 2.7), we assume that the detector response to events producing muons factorises into the product is the energy dispersion of the detector and PSF(φ i |φ µ , E µ ) is its point spread function, assuming these to be identical for muons and antimuons. The energy dispersion is the N c chan response to events that produce muons of a given energy (in contrast to the neutrino N chan response that we employed in Paper I). We obtained this from IceCube detector Monte Carlo simulations (Sec. 2.6).
The uncertainty in the muon reconstruction direction is given on a per-event basis by the IceCube paraboloid sigma σ para,i for the ith event (Sec. 2.5), which accounts for the dependence of the PSF on the incoming muon energy. To obtain the PSF in terms of φ i and φ µ , we shift from the coordinate system centred on the true muon arrival direction (i.e. ∆ = 0 in Eq. 2.1) to the one with the Sun at the origin (φ = 0), integrating over all azimuths to obtain where I 0 is the lowest-order modified Bessel function of the first kind. 2

Predicted event rate
The total predicted number of signal events θ S follows similarly to Eq. 3.11 as the sum of the predicted number of neutrino-initiated signal events (3.15) and the corresponding quantity θ S,ν for anti-neutrinos. Again, we remind the reader that φ µ and E µ are functions of x, y and E. The only difference here with respect to what one would naively read off Eq. 3.11 is the factor L(φ µ , E µ , φ cut ), a dimensionless, energy-dependent angular loss factor that is independent of the muon charge. L corrects for neutrinos that originate from the direction of the Sun but produce muons that are ultimately reconstructed as arriving from outside the analysis cut cone (φ i > φ cut ). Similarly to Paper I, we use the mean angular error of IceCube (σ µ ; cf. Sec. 2.5) to calculate L, integrating the PSF over the analysis cut cone to give This is known as the Marcum P -function or Complementary Marcum Q-function, which we evaluate with the code of Ref. [32]. There are two crucial differences here as compared to Paper I. The first is that L is a muon-level correction factor, expressed in terms of the muon energy and the width of the muon-level angular uncertainty σ µ (E µ ), not the corresponding neutrino quantities. The other is that because of the non-zero muon production angle, the off-centre PSF (Eq. 3.14) must be used instead of the central distribution (Eq. 2.1).
The mean angular error is the correct PSF width to use in Eq. 3.16, because we are interested in determining a priori what fraction of incoming neutrinos with a given energy should be absent from the final set of observed events, due to the chosen angular cut. This is in contrast to the case of the contribution to the partial likelihood coming from the detector response (Eq. 3.14), where we are interested in the probability that a given event originated from the Sun, where the event-level paraboloid σ para,i should be preferred.

Bias factor calculation
The inner double integral in Eq. 3.15 gives the unbiased neutrino effective area for this analysis. It differs from the effective area derived in the standard 79-string analysis [1] in two important ways. First, it includes the factor L, to account for the angular loss due to our analysis cut cone around the solar position. Second, it implicitly assumes that all muons of a given energy are equally likely to pass the original analysis cuts used in the 79-string analysis. In reality, low-energy muons created by high-energy neutrinos are, for example, far more likely to appear in the final event sample than muons of the same energy created by low-energy neutrinos. This is due to the additional light deposited in the detector from the hadronic recoil in the case of a higher-energy neutrino, and the analysis cuts placed on quantities such as the absolute number of activated DOMs.
This departure from a perfect mapping between the properties of a muon and its probability of ending up in the final event sample constitutes a bias that depends on the neutrino energy. This is precisely the reason for the bias factor f b (E) in the preceding expressions, which accounts for the departure of the the event sample from the minimum bias expectation. To quantify this effect, we take the ratio of the original 79-string effective area A eff (E) to the unbiased effective area calculated without the angular correction L. The final effective area in this paper is the product of the bias factor f b and the unbiased effective area with the angular correction L. In this way, our analysis is fully consistent with the original 79-string effective area by construction, and accounts for both the bias and the angular cut cone at the same time.
To facilitate the use of other neutrino spectra, we provide unbiased effective areas precomputed in nulike for the three 79-string IceCube event selections, both with and without the angular correction L. We also provide the routines necessary to repeat the computations  : Limits on dark matter annihilation in the Sun using an analysis that takes into account neutrino energy information. We show limits separately for the three different IC79 event samples SL (summer low) WL (winter low) and WH (winter high) and their combination. The difference between dashed and solid lines indicates the improvement gained by moving from a simple counts-based number likelihood to a full unbinned one, incorporating the number of events, their arrival directions and energies. The full limit is weaker than the WL sample taken alone at low masses, because the SL sample exhibits a weak excess (<2σ local significance) of events above background expectation, not borne out in the WL sample.
for any other dataset. In final likelihood mode, the user can choose to have nulike work with user-supplied bias factors, or use the unbiased effective areas to automatically determine the bias factors. Figures 4-7 show the 90% confidence level (CL) limits on simple effective WIMP DM models computed using IC79 data (Sec. 2) and the nulike 1.0.0 implementation of the likelihood described in Sec. 3. We use the ∆ ln L relative to the background-only prediction as the test statistic, summed over the three event selections, conditioning on all parameters except the cross-section to leave only a single degree of freedom. The distribution of this test statistic is very close to χ 2 , as shown in previous analyses by explicit Monte Carlo [1,22]; this allows CLs to be determined by standard ∆χ 2 methods.  Figure 5: Limits on the spin-dependent WIMP-proton cross-section from IC79 using the improved likelihood, for the canonical soft (bb) and hard (W + W − and τ + τ − ) annihilation channels often seen in SUSY models. Here we compare to the limits from the original IC79 analysis ('PRL'; [1]); note that the previous 'hard' channel limit is W + W − above the W mass, but τ + τ − below it. The addition of energy information provides an improvement of up to a factor of 4 at high WIMP masses over the previous analysis, whereas the limits are in excellent agreement for low WIMP masses. Here we have assumed an annihilation cross-section of σv 0 = 3 × 10 −26 cm 3 s −1 .

Improved limits on WIMP dark matter
For all limits in this section, we assume that DM annihilates exclusively to some specific final state, with a canonical thermal annihilation cross-section σv 0 = 3 × 10 −26 cm 3 s −1 . For all but the highest WIMP masses and lowest scattering cross-sections, these models have reached equilibrium between capture and annihilation in the Sun. We do not assume equilibrium in our calculations however, as is often done. We use DarkSUSY 5.1.3 to compute the predicted neutrino spectrum at the detector for each model, and to solve for the presentday DM population in the Sun. We adopt the standard halo model and default nuclear matrix elements as implemented in DarkSUSY; see discussions in Refs. [17,35]. Fig. 4 presents the limits on the spin-dependent WIMP-proton cross-section imposed by the three different IC79 event samples: WH, WL and SL individually, and in combination. As an example, here we show limits corresponding to annihilation solely to τ + τ − final states. As expected [1,22], the SL and WL samples dominate the sensitivity at low WIMP masses. For comparison, we also show limits based on the number likelihood (Eq. 3.2) alone, neglecting all event-level information. For the cut cone that we use (40 degrees for WL and SL, 20 degrees for WH), considering the arrival directions and energies of neutrino events provides up to a factor of 20 improvement in the resulting limits. Super-K bb IceCube bb IceCube W + W − IceCube τ + τ − PICO-60 (2015) Figure 6: Comparison of our limits with the latest constraints from Super-Kamiokande [2] and PICO [33,34]. Depending on the annihilation channel, IceCube provides the strongest limits above WIMP masses of ∼100-200 GeV. Super-K is more sensitive at the lowest masses. If the annihilation spectrum is soft or heavily suppressed, the PICO experiment provides stronger limits than neutrino telescopes; other direct limits are weaker. Here we have assumed an annihilation cross-section of σv 0 = 3 × 10 −26 cm 3 s −1 for deriving IceCube limits; Super-K limits assume complete equilibrium between capture and annihilation in the Sun.
At high masses, the combined limit in Fig. 4 essentially tracks the exclusion curve of the WH sample, which is orders of magnitude more sensitive than the WL and SL samples in this region of parameter space. At masses below 100 GeV however, where SL and WL both play significant roles, the combined limit is slightly weaker than the limit obtained by considering the WL sample alone. This is because the SL sample exhibits a weak excess above the background expectation inside the analysis cut cone that is not replicated in the WL sample: 819 observed events as compared to 770 predicted in the analysis cone from background alone.
In Fig. 5 we compare these new limits to the previous 79-string IceCube constraints on hard and soft annihilation channels. To allow a reasonable comparison, here we show limits for bb, W + W − and τ + τ − final states, matching what was used in the previous analysis ('soft channel' = bb, 'hard channel' = W + W − for m χ > m W and τ + τ − for m χ < m W ). The previous analysis used the same data as we use here, except that it did not include event energy information in the likelihood function. At low masses, the analysis agrees with the previous one, indicating that the energy information adds little information. Including the event-level energy information has the most impact at high WIMP mass, making use of the relatively good energy resolution of IceCube at high muon energies. The limits in   Figure 7: Limits on the spin-dependent WIMP-proton cross-section from IC79, for a range of different annihilation final states. The canonical hard (W + W − and τ + τ − ) and soft (bb) channels bracket the possible limits for different models reasonably well. More extreme channels (hardest: νν, softest: gg) less often found in SUSY can lead to even stronger or weaker constraints. For the νν channel we have assumed equal branching fractions for all three neutrino flavours. The ability to easily and quickly compute full limits for any combination of final states is a particular feature of the method and tools we present in this paper. As a convenience, datafiles for all curves in this figure are available precomputed in the nulike download 3 . are up to a factor of 4 stronger than the previous analysis at multi-TeV masses. The latest update of WIMPSim fixes an issue with propagation of neutrinos in the Sun that affected the version used to derive the original IC79 limits [1]. This resulted in conservative limits for WIMP masses above ∼500 GeV, ranging from a factor of 1.05 at 500 GeV to 1.2 at 1 TeV and up to 1.5 at 5 TeV for the W + W − and τ + τ − final states. Improvements beyond those factors are due to the improved analysis method in this paper. Fig. 6 compares these limits to other searches for spin-dependent DM-proton scattering, both from the Sun and direct detection experiments. The 79-string IceCube data provide the strongest limits of any search for all masses above ∼100-200 GeV (the exact value depends on the annihilation channel). Super-Kamiokande [2] is the most sensitive experiment at all lower masses. Limits from direct detection [33,34] are weaker, except in the case of DM with soft or suppressed annihilation spectra, in which case the PICO experiment [33,34] is the most constraining. Indirect DM searches by Antares [36] and Baksan [37] have set less stringent limits on the spin-dependent DM-proton scattering and are consequently not included in Fig. 6. Figure 7 shows new limits for all major two-body annihilation final states. Annihilation to either electroweak gauge boson final state is more or less equivalent, as W and Z have around the same mass and couplings to the rest of the SM, and consequently yield very similar neutrino spectra. We don't show hZ, but we have checked that it indeed lies mid-way between hh and ZZ, as expected.
As expected, most channels are indeed bracketed by the canonical 'hard' and 'soft' channels. The exceptions to this are gluon final states, where spectra are especially soft and limits particularly weak, and neutrino final states, which give very strong limits because they are monochromatic at the source. In the neutrino case, the monochromatic source spectrum means that of all final states, annihilation to neutrinos tracks the actual neutrino effective area most closely, with the only deviation from a monochromatic spectrum at the detector coming from reprocessing in the Sun following prompt production at the DM mass. This is also why the neutrino-channel limits at masses above one TeV become weaker than those from the τ + τ − channel: as a channel with an extremely hard annihilation spectrum, most of the neutrinos produced are close to the DM mass, and are therefore absorbed in the Sun. This is a general feature of all channels above one TeV: soft and hard channels begin to swap character in terms of the limits, as softer channels actually produce more neutrinos able to make it out of the Sun and to the detector. This effect can also be seen in the gluon channel limits, which become stronger as the mass increases past ∼7 TeV, as enough of the resulting very low-energy neutrinos are pulled into the observable energy window from below to counteract the slight increase in the number of neutrinos above one TeV that never make it out of the Sun.

Implications for MSSM benchmarks
In this section we use the new IceCube 79-string likelihood to test a number of models of weak-scale supersymmetry, employing the same test statistic as in Sec. 4. Here we focus on the MSSM-25, a 25-parameter, weak-scale parameterisation of the minimal supersymmetric standard model (MSSM; see Ref. [40] for details). This contains the MSSM-19, otherwise known to as the 'phenomenological' (p)MSSM, as a subspace. Fig. 8 shows some MSSM-25 benchmark models from the study of Ref. [40], selected by requiring models with large spin-dependent scattering cross-sections. To give a broader indication of the possibilities in the MSSM, Fig. 8 also shows all models from the MSSM-19 benchmarking exercise of the Snowmass 2013 review [9], except for the Bino-stop coannihilation benchmark, which is very similar to the Bino-squark benchmark in this plane. 4 Except for the models that we show with faded symbols (which we return to later) these models are all consistent with constraints from the LHC, flavour physics and the relic density of dark matter, as well as direct and indirect searches for dark matter. The Snowmass 2013 benchmarks include a 'spoke' of models extending along a single direction in parameter space from one specific benchmark, shown as a vertical line in Fig. 8. We also show shaded bands between the strongest (τ + τ − ) and weakest (bb) limits for channels typically seen in the MSSM. This gives some idea of where essentially all MSSM models are excluded regardless of annihilation channel (above the bb limit), and where only some models are excluded (between Well-tempered neutralinos  [38,39]. These all correspond to so-called 'well-tempered' neutralinos, which exhibit a mixed gaugino-Higgsino character. Solid orange crosses indicate models in tension with IC79 data at more than 1σ (but excluded at less than 90% CL). Green plus symbols indicate models not constrained by IC79, labelled according to the dominant characteristic determining their relic density. The vertical green line corresponds to a benchmark 'spoke' of models [9], where the correct relic density is obtained by bino-squark co-annihilation. Benchmarks are from the MSSM-25 and MSSM-19 ('pMSSM'; a subset of the MSSM-25) scans of Refs. [9,40], and correspond to models allowed by LHC, relic density and other direct and indirect constraints. Benchmark scattering cross-sections are rescaled for the neutralino relic density, and the shaded regions are indicative only; these assume pure spin-dependent scattering and annihilation to the canonical 'hard' and 'soft' channels often seen in the MSSM (even though harder and softer spectra are also possible within the MSSM).
bb and τ + τ − ), depending on their specific annihilation branching fractions to different final states.
We have colour-coded the individual models in Fig. 8 by the extent to which they are excluded by the new IceCube limits, taking into account both spin-dependent and spinindependent scattering. We have also labelled different benchmark groups according to the means by which the neutralino achieves the appropriate relic density in the early Universe. Many neutralino models are excluded for the first time by the new limits we present here (bright red crosses). Other models exhibit a tension with data at the 68-90% confidence level (orange crosses). These are 'well-tempered' neutralino models, which exhibit a roughly even mixture of gaugino and Higgsino weak eigenstates, boosting their spin-dependent scattering cross-section without contributing too strongly to the spin-independent one. Other benchmarks (green plus symbols), where the relic density is achieved by squark or chargino co-annihilation with the neutralino, resonant annihilation via the CP-odd Higgs, or by virtue of the large annihilation cross-section exhibited by pure Higgsinos, remain unconstrained by spin-dependent searches of any kind.
We also show a number of well-tempered neutralino benchmarks with faded symbols in Fig. 8, indicating that although they were consistent with all earlier data, they have since been excluded by LUX [39]. One of these examples (the well-tempered neutralino MSSM-19 benchmark from Ref. [9]) was already strongly excluded by the original LUX spin-independent limits. The others satisfy the spin-independent limit, but are excluded by the recent LUXCalc [38] application of the LUX data to spin-dependent neutron scattering. 5 All of these models are strongly excluded by IceCube.
In Table 1, we give further details of all the benchmark models shown in Fig. 8. These include cross-sections for annihilation and nuclear scattering ( σv , σ SD , σ SI ), relic densities (Ωh 2 ), capture and annihilation rates (C, A), and dominant annihilation branching fractions (necessary to understand differences between the various well-tempered models).
The benchmark models we show here, whilst illustrative, are only isolated samples from the vast range of possible models in the MSSM. A full statistical analysis of MSSM theories in the context of these data awaits their inclusion in large-scale global fits, as expected shortly from the GAMBIT Collaboration [42].

Conclusions
We have presented a new analysis of data collected in the 79-string IceCube search for dark matter, taking into account energies of individual neutrino events. This resulted in stronger spin-dependent limits on WIMP dark matter, in particular for high WIMP masses, and allowed us to rule out a number of MSSM models for the first time. In the process, we developed an updated fast likelihood pipeline for event-level neutrino telescope DM search data, allowing it to be quickly and accurately applied to constrain essentially any dark matter model. We have also provided a public code implementing the new likelihood (nulike), and made data from the 79-string IceCube DM search publicly available in a format compatible with its use. Full details of the SUSY benchmarks and generic WIMP results presented in this paper are available as example programs in the public distribution of nulike. Future improvements can be expected from applications of nulike to other models, and from the 86string IceCube search for dark matter, which will include additional data and an improved energy proxy. Table 1: Properties of the benchmark models shown in Fig. 8.