---
author:
- 'P. Scott,'
- 'C. Savage,'
- 'J. Edsjö'
- 'and\'
- |
    The IceCube Collaboration:\
    R. Abbasi$^{3}$, Y. Abdou$^{4}$, M. Ackermann$^{5}$, J. Adams$^{6}$, J. A. Aguilar$^{7}$, M. Ahlers$^{3}$, D. Altmann$^{8}$, K. Andeen$^{3}$, J. Auffenberg$^{3}$, X. Bai$^{9,10}$, M. Baker$^{3}$, S. W. Barwick$^{11}$, V. Baum$^{12}$, R. Bay$^{13}$, K. Beattie$^{14}$, J. J. Beatty$^{15,16}$, S. Bechet$^{17}$, J. Becker Tjus$^{18}$, K.-H. Becker$^{19}$, M. Bell$^{20}$, M. L. Benabderrahmane$^{5}$, S. BenZvi$^{3}$, J. Berdermann$^{5}$, P. Berghaus$^{5}$, D. Berley$^{21}$, E. Bernardini$^{5}$, D. Bertrand$^{17}$, D. Z. Besson$^{22}$, D. Bindig$^{19}$, M. Bissok$^{23}$, E. Blaufuss$^{21}$, J. Blumenthal$^{23}$, D. J. Boersma$^{23}$, C. Bohm$^{2}$, D. Bose$^{24}$, S. Böser$^{25}$, O. Botner$^{26}$, L. Brayeur$^{24}$, A. M. Brown$^{6}$, R. Bruijn$^{27}$, J. Brunner$^{5}$, S. Buitink$^{24}$, K. S. Caballero-Mora$^{20}$, M. Carson$^{4}$, J. Casey$^{28}$, M. Casier$^{24}$, D. Chirkin$^{3}$, B. Christy$^{21}$, F. Clevermann$^{29}$, S. Cohen$^{27}$, D. F. Cowen$^{20,30}$, A. H. Cruz Silva$^{5}$, M. Danninger$^{2,a}$, J. Daughhetee$^{28}$, J. C. Davis$^{15}$, C. De Clercq$^{24}$, F. Descamps$^{3}$, P. Desiati$^{3}$, G. de Vries-Uiterweerd$^{4}$, T. DeYoung$^{20}$, J. C. D[í]{}az-Vélez$^{3}$, J. Dreyer$^{18}$, J. P. Dumm$^{3}$, M. Dunkman$^{20}$, R. Eagan$^{20}$, J. Eisch$^{3}$, R. W. Ellsworth$^{21}$, O. Engdeg[å]{}rd$^{26}$, S. Euler$^{23}$, P. A. Evenson$^{9}$, O. Fadiran$^{3}$, A. R. Fazely$^{31}$, A. Fedynitch$^{18}$, J. Feintzeig$^{3}$, T. Feusels$^{4}$, K. Filimonov$^{13}$, C. Finley$^{2}$, T. Fischer-Wasels$^{19}$, S. Flis$^{2}$, A. Franckowiak$^{25}$, R. Franke$^{5}$, K. Frantzen$^{29}$, T. Fuchs$^{29}$, T. K. Gaisser$^{9}$, J. Gallagher$^{32}$, L. Gerhardt$^{14,13}$, L. Gladstone$^{3}$, T. Glüsenkamp$^{5}$, A. Goldschmidt$^{14}$, J. A. Goodman$^{21}$, D. Góra$^{5}$, D. Grant$^{33}$, A. Gro[ß]{}$^{34}$, S. Grullon$^{3}$, M. Gurtner$^{19}$, C. Ha$^{14,13}$, A. Haj Ismail$^{4}$, A. Hallgren$^{26}$, F. Halzen$^{3}$, K. Hanson$^{17}$, D. Heereman$^{17}$, P. Heimann$^{23}$, D. Heinen$^{23}$, K. Helbing$^{19}$, R. Hellauer$^{21}$, S. Hickford$^{6}$, G. C. Hill$^{35}$, K. D. Hoffman$^{21}$, R. Hoffmann$^{19}$, A. Homeier$^{25}$, K. Hoshina$^{3}$, W. Huelsnitz$^{21,36}$, P. O. Hulth$^{2}$, K. Hultqvist$^{2}$, S. Hussain$^{9}$, A. Ishihara$^{37}$, E. Jacobi$^{5}$, J. Jacobsen$^{3}$, G. S. Japaridze$^{38}$, O. Jlelati$^{4}$, H. Johansson$^{2}$, A. Kappes$^{8}$, T. Karg$^{5}$, A. Karle$^{3}$, J. Kiryluk$^{39}$, F. Kislat$^{5}$, J. Kläs$^{19}$, S. R. Klein$^{14,13}$, J.-H. Köhne$^{29}$, G. Kohnen$^{40}$, H. Kolanoski$^{8}$, L. Köpke$^{12}$, C. Kopper$^{3}$, S. Kopper$^{19}$, D. J. Koskinen$^{20}$, M. Kowalski$^{25}$, M. Krasberg$^{3}$, G. Kroll$^{12}$, J. Kunnen$^{24}$, N. Kurahashi$^{3}$, T. Kuwabara$^{9}$, M. Labare$^{24}$, K. Laihem$^{23}$, H. Landsman$^{3}$, M. J. Larson$^{41}$, R. Lauer$^{5}$, M. Lesiak-Bzdak$^{39}$, J. Lünemann$^{12}$, J. Madsen$^{42}$, R. Maruyama$^{3}$, K. Mase$^{37}$, H. S. Matis$^{14}$, F. McNally$^{3}$, K. Meagher$^{21}$, M. Merck$^{3}$, P. Mészáros$^{30,20}$, T. Meures$^{17}$, S. Miarecki$^{14,13}$, E. Middell$^{5}$, N. Milke$^{29}$, J. Miller$^{24}$, L. Mohrmann$^{5}$, T. Montaruli$^{7,43}$, R. Morse$^{3}$, S. M. Movit$^{30}$, R. Nahnhauer$^{5}$, U. Naumann$^{19}$, S. C. Nowicki$^{33}$, D. R. Nygren$^{14}$, A. Obertacke$^{19}$, S. Odrowski$^{34}$, A. Olivas$^{21}$, M. Olivo$^{18}$, A. O’Murchadha$^{17}$, S. Panknin$^{25}$, L. Paul$^{23}$, J. A. Pepper$^{41}$, C. Pérez de los Heros$^{26}$, D. Pieloth$^{29}$, N. Pirk$^{5}$, J. Posselt$^{19}$, P. B. Price$^{13}$, G. T. Przybylski$^{14}$, L. Rädel$^{23}$, K. Rawlins$^{44}$, P. Redl$^{21}$, E. Resconi$^{34}$, W. Rhode$^{29}$, M. Ribordy$^{27}$, M. Richman$^{21}$, B. Riedel$^{3}$, J. P. Rodrigues$^{3}$, F. Rothmaier$^{12}$, C. Rott$^{15}$, T. Ruhe$^{29}$, D. Rutledge$^{20}$, B. Ruzybayev$^{9}$, D. Ryckbosch$^{4}$, S. M. Saba$^{18}$, T. Salameh$^{20}$, H.-G. Sander$^{12}$, M. Santander$^{3}$, S. Sarkar$^{45}$, K. Schatto$^{12}$, M. Scheel$^{23}$, F. Scheriau$^{29}$, T. Schmidt$^{21}$, M. Schmitz$^{29}$, S. Schoenen$^{23}$, S. Schöneberg$^{18}$, L. Schönherr$^{23}$, A. Schönwald$^{5}$, A. Schukraft$^{23}$, L. Schulte$^{25}$, O. Schulz$^{34}$, D. Seckel$^{9}$, S. H. Seo$^{2}$, Y. Sestayo$^{34}$, S. Seunarine$^{46}$, M. W. E. Smith$^{20}$, M. Soiron$^{23}$, D. Soldin$^{19}$, G. M. Spiczak$^{42}$, C. Spiering$^{5}$, M. Stamatikos$^{15,47}$, T. Stanev$^{9}$, A. Stasik$^{25}$, T. Stezelberger$^{14}$, R. G. Stokstad$^{14}$, A. Stö[ß]{}l$^{5}$, E. A. Strahler$^{24}$, R. Ström$^{26}$, G. W. Sullivan$^{21}$, H. Taavola$^{26}$, I. Taboada$^{28}$, A. Tamburro$^{9}$, S. Ter-Antonyan$^{31}$, S. Tilav$^{9}$, P. A. Toale$^{41}$, S. Toscano$^{3}$, M. Usner$^{25}$, N. van Eijndhoven$^{24}$, D. van der Drift$^{14,13}$, A. Van Overloop$^{4}$, J. van Santen$^{3}$, M. Vehring$^{23}$, M. Voge$^{25}$, C. Walck$^{2}$, T. Waldenmaier$^{8}$, M. Wallraff$^{23}$, M. Walter$^{5}$, R. Wasserman$^{20}$, Ch.  Weaver$^{3}$, C. Wendt$^{3}$, S. Westerhoff$^{3}$, N. Whitehorn$^{3}$, K. Wiebe$^{12}$, C. H. Wiebusch$^{23}$, D. R. Williams$^{41}$, H. Wissing$^{21}$, M. Wolf$^{2}$, T. R. Wood$^{33}$, K. Woschnagg$^{13}$, C. Xu$^{9}$, D. L. Xu$^{41}$, X. W. Xu$^{31}$, J. P. Yanez$^{5}$, G. Yodh$^{11}$, S. Yoshida$^{37}$, P. Zarzhitsky$^{41}$, J. Ziemann$^{29}$, A. Zilles$^{23}$, and M. Zoll$^{2}$
bibliography:
- 'DMbiblio.bib'
- 'SUSYbiblio.bib'
title: 'Use of event-level neutrino telescope data in global fits for theories of new physics'
---

Introduction {#intro}
============

Despite ongoing efforts, we have yet to identify dark matter.  One of the most promising candidate classes is the weakly-interacting massive particle (WIMP) [@Bergstrom00; @Bertone05; @BertoneBook], so named because WIMPs interact with standard model (SM) particles only via the weak nuclear force.  WIMPs are appealing because the expected cosmological density of a thermal relic with a weak-scale annihilation cross-section is the same order of magnitude as the observed value [@WMAP7]. 

WIMPs are typically sought via three observational channels: direct WIMP-nucleon scattering (e.g. [@CerdenoGreen10; @Pato11; @CoGeNTAnnMod11; @CRESST11; @XENON100; @CDMSLowE12; @COUPP12]), production at accelerators (e.g. [@White07; @Bai10; @Goodman10; @Fox12a; @Fox12b; @CMS12a; @Frandsen12]) and indirect detection of SM products of WIMP self-annihilation (e.g. [@Pamelapositron; @CompositeAlt; @LATDwarfComposite; @MAGICSegue; @VERITASSegue; @IceCube09; @IceCube09_KK]).  A promising version of indirect detection is to use neutrino telescopes such as IceCube [@IC40DM], ANTARES [@Lim09] and SuperKamiokande [@SuperK11] to search for high-energy neutrinos produced by annihilation of WIMPs in the core of the Sun or Earth [@Press85; @Krauss:1985ApJ; @Freese86; @Krauss86; @Gaisser86].  Due to their weak interactions with nuclei, such WIMPs would have scattered on solar nuclei and become gravitationally bound to the solar system, eventually returning to scatter repeatedly and settle (in the case of the Sun) to the solar core [@Gould87b] (see also [@Scott09] for a review of this process and impacts of WIMPs on stellar structure, and [@Gould91; @Lundberg04; @Peter09a; @Peter09; @Sivertsson12] for additional corrections due to the influence of planets). 

WIMPs arise in many proposed extensions to the SM.  Standard neutrino telescope analyses [@IceCube09; @IceCube09_KK; @ICRC2011ic86; @IC40DM; @Lim09; @SuperK11; @Kumar12] take an effective view of WIMP interactions, placing limits on WIMP-nucleon scattering cross-sections as a function of the WIMP mass, by assuming a certain annihilation cross-section and final state.  These limits are difficult to translate into actual particle models, where the annihilation cross-section and branching fractions may take on a range of different values for any given WIMP mass and nuclear-scattering cross-section.  To properly interpret limits on neutrino fluxes in terms of the parameters (defined at e.g. the Lagrangian level) of a theory for new physics, it becomes necessary to compare the observed neutrino flux with the predicted neutrino signal for each individual point in the parameter space of the theory. 

Having translated the flux limits into direct constraints on the parameters of a theory, it is then also possible to compare and combine the sensitivities of multiple experiments, even if they probe entirely different sectors of the theory (e.g. neutrino and accelerator searches).  Within the realm of theories for new weak-scale physics, this ‘global fit’ approach has so far been applied mostly to supersymmetry (SUSY) [@Baltz04; @SFitter; @Allanach06; @Ruiz06; @Trotta08; @Fittino; @Abdussalam09a; @Abdussalam09b; @Mastercode12], in the form of the minimal supersymmetric standard model (MSSM).  Recent analyses have included new data from the Large Hadron Collider (LHC) [@MastercodeHiggs; @SuperBayeSHiggs; @Gfitter11], direct [@BertoneLHCDD; @Akrami11DD; @MastercodeXENON100; @SuperBayeSXENON100] and indirect detection [@Scott09c; @Ripken11; @Fittino12].  Whilst great care needs to be taken over the details of the statistical methods employed [@Akrami09; @SBSpike; @SBcoverage; @Akrami11coverage; @Strege12], these analyses have proven a clear success, pointing the way to a future of closer comparison between astronomical and terrestrial experiments. 

To date global fits have not included neutrino telescope data.  The ability of future incarnations of IceCube to detect WIMP annihilation in the constrained MSSM (CMSSM) has been studied [@Trotta09; @Nightmare; @Roszkowski12], based on the number of observed events only and simple estimates of the instrumental sensitivity.  Other authors have looked at predicted neutrino rates in 2D slices through MSSM parameter spaces (e.g. [@Feng01; @Barger02; @Baer04; @Ellis09; @Ellis11]), or from sets of models in non-statistical random scans of various incarnations of the MSSM (e.g. [@Bergstrom98b; @Wikstrom09; @Cotta12]). 

Here we present a fast likelihood method for including full event-level data from neutrino telescopes in global fits and related analyses.  In particular, our formulation allows the directions and energy estimators associated with each event to be included in the final unbinned likelihood calculation, which can then be employed as a likelihood component in a global fit.  The inclusion of spectral information in neutrino searches for WIMPs has been of particular interest recently [@Rott11; @Allahverdi12].  As a by-product of our approach, we also present a rigorous but simple exclusion measure for individual models, based just on the observed number of events in a neutrino telescope.  We give full details of the input data required, and explicit examples using real data from the IceCube neutrino telescope. 

The IceCube data and simulations we employ are described in Section \[icecube\].  These data have been made publicly available on the IceCube webserver [@filesloc] and in .[^1] The likelihood formalism, outlined in Section \[likelihood\], is implemented and also available in .  Our example global fits and MSSM scans are detailed in Section \[examples\].  We summarise our method and examples in Section \[conclusions\]. 

The IceCube Neutrino Telescope {#icecube}
==============================

Description {#icdescription}
-----------

Completed on December 18 2010, the IceCube neutrino observatory [@perfomancePaperIC] has $5160$ digital optical modules (DOMs) installed on $86$ strings between $1450$m and $2450$m below the surface in the glacial ice at the South Pole.  The detector consists of a hexagonal grid with horizontal spacing between strings of $125$m and a vertical spacing between DOMs of $17$m, leading to a total instrumented volume of $1$ km$^{3}$. Eight of the 86 strings were deployed in a more densely instrumented core in the middle of the array, with average inter-string separation of 42m and vertical DOM separation of 7m. Together with the 12 adjacent standard IceCube strings, this forms the DeepCore subarray.  DeepCore increases the sensitivity of IceCube at low energies, and substantially lowers the energy threshold.  IceCube records Čerenkov light in the ice from relativistic charged particles created in neutrino interactions in or near the detector.  By recording the arrival times and locations of these photons with DOMs, the direction and energy of the muon, and identity of its parent neutrino, can be reconstructed. 

Data samples {#icevent}
------------

For this paper, we use the analysis details of a search for WIMP dark matter annihilation in the Sun with the IceCube $22$-string configuration [@IceCube09; @filesloc].  This data set has $104.3$ days of live time, and was recorded between June 1 and September 23 $2007$. For the results in Sections \[cmssmglobal\] and \[cmssmrecon\] we use the event-level data at final analysis level and signal simulations from [@IceCube09].  As described in Ref.  [@IceCube09], final analysis level is reached through a series of increasingly stringent event selections that are applied to remove the cosmic-ray induced backgrounds.  Events are selected which appear to originate from the direction of the Sun and which exhibit low energy signatures. 

The results discussed in Section \[mssm7\] are based on a detailed study to determine the sensitivity of the $86$-string detector to signals originating from dark matter annihilation in the centre of the Sun [@ICRC2011ic86].  This study was performed as a full analysis in all details and gives a realistic expectation of the capabilities of IceCube to observe dark matter-induced signals, given the state of data extraction, reconstruction, and signal discrimination techniques available at the time of the study. 

Signal and background simulation {#icsim}
--------------------------------

Previous work [@IceCube09; @ICRC2011ic86] simulated solar WIMP signals using WimpSim [@Blennow08; @wimpsimweb], which describes the annihilation of WIMPs inside the Sun.  WimpSim simulates the production, interaction, oscillation and propagation of neutrinos from the core of the Sun to the detector, in a fully consistent three-flavour way.  For the previous IceCube solar WIMP analyses that we consider (IceCube $22$-string data analysis [@IceCube09] and $86$-string sensitivity analysis [@ICRC2011ic86]), two annihilation channels were simulated: $\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}\to W^{+}W^{-}$ (the ‘hard’ channel) and $\tilde{\chi}^{0}_{1}\tilde{\chi}^{0}_{1}\to b\bar{b}$ (the ‘soft’ channel).  These two channels were chosen in an attempt to approximately cover the range of possible SUSY models, by assuming 100% branching into two channels with very different characteristics.  For the $22$-string analysis, neutralino masses $m_{\chi^{0}_{1}}$ from $250$GeV to $5000$GeV were simulated.  For the $86$-string analysis, this range was extended down to $m_{\chi^{0}_{1}} = 50$GeV. A similar exercise was performed for Kaluza-Klein WIMPs in models of Universal Extra Dimensions [@IceCube09_KK]. 

In contrast, the analysis method we describe here (Section \[likelihood\]) is developed to employ the specific annihilation spectrum appropriate for each model in a given theory of new physics.  We use various parameterisations of the MSSM as examples, but the analysis technique we present is applicable to any theory containing a WIMP.  To calculate neutrino annihilation spectra for individual SUSY models, we create yield tables from a series of WimpSim simulations, with neutralino masses from 3GeV to 10TeV and many different annihilation channels.  For each SUSY model, we interpolate in these tables and add the contributions to the total neutrino spectrum from each annihilation channel, according to the partial annihilation cross section for the given channel.  We include more complicated annihilation channels involving model-dependent decays, such as those of Higgs bosons, by summing up the contributions from their decay in flight.  This procedure is included in DarkSUSY [@darksusy], which we use to calculate the detailed neutrino spectra model by model. 

Background contributions for this analysis are muon events from single and coincident air showers, as well as atmospheric neutrinos.  No dedicated background simulations are needed for the work we present here, as we estimate the expected background from scrambling real data at the final analysis level (detailed within Section \[icbackground\]). 

Effective area calculation {#iceffarea}
--------------------------

We derived detector responses for the 22-string and 86-string IceCube configurations mentioned in Section \[icdescription\] from the sets of signal simulations used in each analysis [@IceCube09; @ICRC2011ic86].  In each case, we calculated an effective area for detection of muon neutrinos by IceCube from the direction of the Sun as a function of neutrino energy.  These effective areas correspond to averages over the austral winter.  The 22-string area is available online [@filesloc].  For convenience, both datafiles are also redistributed in .  The effective area for muon neutrinos and muon anti-neutrinos are given separately. 

Detector systematics determined from the analyses are energy-dependent; the data files therefore include systematic uncertainties on the detector response within each energy bin, to which we assign a $1\sigma$ confidence level.  These uncertainties were determined within 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 the effective areas, at the $1\sigma$ confidence level. 

Angular response {#icangular}
----------------

The point spread function (PSF) describes the uncertainty in the reconstructed arrival directions of incoming neutrinos.  The full PSF is a 2D joint probability distribution for the angular separation between the true arrival direction and the reconstructed one, in each of two linearly-independent directions on the sky.  For the small angles under consideration and suitably chosen axes, this is well approximated by a 2D Gaussian distribution of the form $$\label{PSF2}
P(\theta_1,\theta_2) = \frac{1}{2\pi \sigma_1 \sigma_2} \, \exp\left[-\frac{\theta_1^2}{2\sigma_1^2}-\frac{\theta_2^2}{2\sigma_2^2}\right] \, ,$$ where $\theta_1$ and $\theta_2$ are the angular separation along each axis.  Given that the Sun is a circularly symmetric source, we work exclusively with a reduction of the PSF into the 1-dimensional probability distribution for the absolute angular separation on the sky, $\phi\equiv\sqrt{\theta_1^2+\theta_2^2}$, removing the azimuthal dependence.  For an azimuthally-symmetric PSF, $\sigma = \sigma_1 = \sigma_2$ and the 1D PSF is $$\label{PSF}
P(\phi) = \frac{\phi}{\sigma^2} \, \exp\left[-\frac{\phi^2}{2\sigma^2}\right].$$

Over a given energy range, it is possible to accurately construct the 1D PSF directly for neutrinos arriving from all directions, using IceCube data analysis signal simulations [@IceCube09; @ICRC2011ic86].  In that case, one need not assume that the 2D PSF is even Gaussian, let alone azimuthally symmetric.  However, the radial probability density function (PDF) of an azimuthally-symmetric 2D Gaussian is a reasonable approximation to the radial PDF of a moderately non-symmetric 2D Gaussian, and has the added advantage that it can be parameterised in terms of a single $\sigma$. As the 2D PSF in IceCube is not strongly asymmetric, we therefore model the 1D PSF with Eq.  \[PSF\], but extract the parameter $\sigma$, which we refer to as the “mean angular error”, directly from the 1D PSF constructed from IceCube signal simulations.  This can be done by integrating to the containment angle implied by e.g. the mean ($\langle\phi\rangle=\sqrt{\frac{\pi}{2}}\sigma$), second moment ($\langle\phi^2\rangle=2\sigma^2$), or median ($\sigma \sqrt{2\ln{2}}$; this is what we use here) of the $\phi$ distribution in the signal simulations.  As in previous work [@IceCube09; @ICRC2011ic86], we consider $\sigma$ over the same energy bins as used for calculating the effective area. 

We associate angular uncertainties with real data events on an event-by-event basis, using the paraboloid method [@paraboloid].  A paraboloid function is fitted to the muon track reconstruction likelihood function in the neighbourhood of the best fit.  The resulting confidence ellipse on the sky is represented by the axes $\sigma_{1}$ and $\sigma_{2}$, which correspond to the standard deviations of the likelihood function in each of two linearly-independent directions.  The overall reconstructed likelihood track uncertainty, $\sigma_{\mathrm{para}}$ (the “paraboloid sigma”), is calculated as the mean in quadrature of the uncertainties in the two axes, $\sigma_{\mathrm{para}}^2=(\sigma_{1}^{2} + \sigma_{2}^{2})/2$. Good track fits generally result in a paraboloid that is narrow along both axes, and therefore have small $\sigma_{\mathrm{para}}$ values. 

The simulated distribution of reconstructed angles, used for determining the mean angular error, contains contributions from both the intrinsic angular deviation of neutrino-induced muons from the incoming neutrino direction, and the error in the reconstructed muon direction itself.  For actual data events, the former cannot be known except in an average sense, whereas the latter can be estimated for each event with the paraboloid sigma.  For the 22-string configuration, the reconstruction uncertainty dominates over the intrinsic angular deviation of neutrino-induced muons, so we are able to employ the paraboloid sigma $\sigma_{\mathrm{para}}$ as an accurate approximation to the total angular error in Eq.  \[PSF\], for individual events.  For future IceCube configurations or other experiments where the muon angular reconstruction is of higher precision than the intrinsic muon angular spread, an additional component must be added to $\sigma_{\mathrm{para}}$ for each event, in order to properly account for the difference between the incoming neutrino and neutrino-induced muon directions. 

![ Predicted probability distributions of $N_\mathrm{chan}$ derived from high-statistics $\nu$ simulations generated in the initial $22$-string analysis [@IceCube09].  Each distribution is constructed for neutrinos having energies in a specific logarithmic energy interval of width 0.2.  The fitted Gaussian functions are to guide the eye only, and are not used in our calculations.[]{data-label="fig:nchan"}](Nchan_12plots){width="90.00000%"}

Energy estimator {#icenergy}
----------------

In the broader context of IceCube data analysis, the study we present here deals with comparatively low-energy neutrinos.  The corresponding muon events are associated with minimum-ionising tracks, where stochastic losses, which increase with muon energy, are not dominant.  Energy estimators, as generally used within IceCube analyses targeting higher muon energies, exploit the energy-dependence of such losses in order to estimate the energy of the muon and the parent neutrino.  The volume instrumented with the $22$-string configuration is too small to use the reconstructed track length of either fully or partially contained muon tracks to determine the energy.  Instead, we use the number of lit DOMs $N_\mathrm{chan}$, a measure of the amount of recorded light per event, as a suitable energy estimator.  We make no attempt to assign a specific energy and associated uncertainty to each event according to its $N_\mathrm{chan}$ value.  Instead, we calculate the expected distribution of observed $N_\mathrm{chan}$ values for a series of intervals in neutrino energy, and use these together with the predicted energy spectrum of the signal for a given SUSY model, to calculate the predicted distribution of $N_\mathrm{chan}$. Using an unbinned likelihood (Section \[likelihood\]), we then simply compare the observed value of $N_\mathrm{chan}$ for each event to this predicted distribution.  We derived the probability distributions per $\nu$-energy interval from high-statistics $\nu$ simulations generated in the initial $22$-string analysis [@IceCube09].  These distributions are shown for the energy range relevant to our analysis in Fig.  \[fig:nchan\].  Note that we do not use the fitted Gaussian functions in Fig.  \[fig:nchan\] for any purpose other than to guide the eye; we employ the actual distributions directly for our signal predictions and likelihood calculations. 

Background estimation {#icbackground}
---------------------

We estimate the necessary background distributions from actual data.  For the likelihood functions in Section \[likelihood\], we need two background distributions: the angular distribution of background events ${\mathrm{d}}P_\mathrm{BG}(\phi')/{\mathrm{d}}\phi'$ — given as a function of $\phi'$, the angle between the reconstructed track direction and the Sun — and the distribution of $N_\mathrm{chan}$ due to background events, ${\mathrm{d}}P_\mathrm{BG}(N_\mathrm{chan})/{\mathrm{d}}N_\mathrm{chan}$.

The angular distribution of the dominant cosmic-ray shower background contribution is expected to be azimuthally independent, depending only upon the angle from the horizon (zenith angle).  Thus, events incident in a horizontal band across the sky covering the same vertical range as the analysis (centred upon the solar zenith), but including all 360$^{\circ}$ around the horizon, are representative of the backgrounds in the small angular region around the Sun actually used in the analysis.  This allows for an estimate of the necessary background distributions driven purely by observation.[^2] The background is roughly constant per unit angular area in this band, with any mild zenith angle dependence softened by an averaging over the zenith angle of the Sun during the austral winter.  Hence $\frac{{\mathrm{d}}P_{\mathrm{BG}}}{{\mathrm{d}}\phi'} \sim \sin\phi'$, though the actual $\frac{{\mathrm{d}}P_{\mathrm{BG}}}{{\mathrm{d}}\phi'}$ we use is tabulated rather than fit to any functional form. 

Data format and availability {#icdata}
----------------------------

The likelihood formalism described in the following Section, as implemented within , requires four simple text files containing IceCube data for any given analysis.  The first data file gives the binned $\nu$ and $\bar{\nu}$ effective areas as a function of energy, along with the associated $1\sigma$ statistical and systematic errors in each bin.  This file also provides the mean angular error.  The second data file contains all events that a) passed the event selection corresponding to the effective areas, and b) have reconstructed directions that fall within some specified cut cone around the solar position.  The 22-string version of this file consists of 180 events within 10 degrees of the Sun.  This file includes the observed value of $N_\mathrm{chan}$ for each event, along with the reconstructed arrival angle relative to the Sun, and the corresponding paraboloid sigma.  This file also includes the total live (i.e. exposure) time of observations towards the Sun.  The third file gives the background distributions, with the angular distribution provided as a function of $\phi'$ and the spectral distribution as a function of $N_\mathrm{chan}$. The final file contains the $N_\mathrm{chan}$ response histograms presented in Fig.  \[fig:nchan\]. 

We provide such files in  for both the 22- and 86-string examples that we show in this paper.  The 22-string files (containing real data), are constructed from the original IceCube analysis [@IceCube09] and its associated data release [@filesloc].  The 86-string files (containing simulated data), are constructed from the IceCube sensitivity analysis [@ICRC2011ic86] and the expected number of all-sky background events [@DanningerLicThesis].[^3] It is anticipated that 79-string data will later also be made available on the web and in DarkSUSY, in connection with an analysis of those data using the techniques presented in this paper. 

Likelihood functions {#likelihood}
====================

Using the data, responses and simulations described in Section \[icecube\], our goal is to evaluate to what degree IceCube observations support or constrain a theory of new physics $\Psi$, given some particular values ${\psi}$ of the theory’s $m$ free parameters.  We refer to a specific new physics scenario $\Psi$ as a ‘theory’ and a specific choice of its parameters ${\psi}$ as a ‘model’.  We may be interested in

-   the absolute goodness of fit of the model with parameter set ${\psi}$, or

-   the goodness of fit of ${\psi}$ relative to other points in the parameter space of the theory. 

Case a) is model exclusion: we are interested in the maximal confidence level at which we can exclude the theory with specific parameter values ${\psi}$ as a true hypothesis.  Case b) is most often parameter estimation: we are interested in determining how well ${\psi}$ fits the data, compared to the best-fit point $\hat{{\psi}}$ within the parameter space of the theory.  Assuming the overall theory to be correct, we use this information to determine a region of parameter space that includes the values of the parameters, with some specific level of confidence. 

The absolute goodness of fit for a model is described by its $p$-value, where the confidence level with which the model can be excluded is $1-p$. The $p$-value can be obtained from the likelihood function ${\mathcal{L}}({D}|{\psi})$, where ${\mathcal{L}}({D}|{\psi})\,\mathrm{d}{D}$ describes the probability of observing data in the range $[{D}, {D}+\mathrm{d}{D}]$ if $\Psi$ is true and the true values of its parameters are ${\psi}$. The relative goodness of fit can be obtained from the ratio $$\Lambda({D}|{\psi})\equiv\frac{{\mathcal{L}}({D}|{\psi})}{{\mathcal{L}}_\mathrm{max}({D}|\hat{\psi})}.$$ Here ${\mathcal{L}}_\mathrm{max}$ is the maximum likelihood, found at the best-fit point $\hat{{\psi}}$. Calculation of $p$ or $\Lambda$ requires knowledge of the distribution of ${\mathcal{L}}({D}|{\psi})$ for repeated experiments (i.e. the distribution in data space).  The simplest approach is to assume that $-2\ln\Lambda$ follows a $\chi^2$ distribution with $m$ degrees of freedom (as predicted in the limit of infinite data by Wilks’ Theorem [@Wilks]), whilst the most rigorous approach is to explicitly construct the distribution by brute force simulation (the Feldman-Cousins approach [@FeldmanCousins]). 

For both model exclusion and parameter estimation, we require an expression for ${\mathcal{L}}({D}|{\psi})$. Here we develop a likelihood function for IceCube data, based on the total number of observed events and the individual properties of each event.  We go on to develop a $p$-value for model exclusion based on our likelihood construction.  Although we work exclusively at neutrino level in this paper (cf.  Section \[icevent\]), the treatment in this Section should be equally valid at muon level. 

General unbinned likelihood
---------------------------

Consider a set of ${{n_\mathrm{tot}}}$ events observed in IceCube.  Denote the true energy of the $i$th event as ${E_i}$, and its arrival angle relative to the position of the Sun on the sky as ${\phi_{i}}$. The $i$th event will be reconstructed with arrival angle ${\phi'_i}$, and will cause a number of DOMs ${N_i}$ to fire.  The unbinned likelihood for these data is $$\label{unbinned_like}
{\mathcal{L}}_\mathrm{unbin} \equiv {\mathcal{L}}({{n_\mathrm{tot}}}|\theta_\mathrm{tot})\prod_{i=1}^{{{n_\mathrm{tot}}}} \int_0^\pi\int_0^\infty Q({N_i}, {\phi'_i}| {E_i}, {\phi_{i}})\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}\,{\mathrm{d}}{\phi_{i}}}({E_i},{\phi_{i}},{\psi})\,{\mathrm{d}}{E_i}\, {\mathrm{d}}{\phi_{i}}.$$ Here $Q({N_i}, {\phi'_i}| {E_i}, {\phi_{i}})$ is the probability per unit angle and $N_\mathrm{chan}$ (i.e. probability density) for observing ${N_i}$ and ${\phi'_i}$ for the $i$th event when the true values of the energy and arrival angle are ${E_i}$ and ${\phi_{i}}$. The *a priori* probability density for the $i$th event to arrive with energy ${E_i}$ from angle ${\phi_{i}}$ is given by $\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}\,{\mathrm{d}}{\phi_{i}}}({E_i},{\phi_{i}},{\psi})$; this is a prediction of the model ${\psi}$.

Number likelihood
-----------------

Neglecting systematic issues, which we address below, the prefactor ${\mathcal{L}}({{n_\mathrm{tot}}}|\theta_\mathrm{tot})$ in Eq.  \[unbinned\_like\] is the standard Poissonian likelihood, equal to the probability of observing ${{n_\mathrm{tot}}}$ events produced by a Poisson process with mean $\theta_\mathrm{tot}$ $$\label{standard}
{\mathcal{L}}({{n_\mathrm{tot}}}|\theta_\mathrm{tot}) = \frac{\theta_\mathrm{tot}^{{n_\mathrm{tot}}}e^{-\theta_\mathrm{tot}}}{{{n_\mathrm{tot}}}!}.$$ The expectation value for the number of events is $\theta_\mathrm{tot}\equiv\theta_\mathrm{tot}({\psi})$, the total number of events predicted by model ${\psi}$. This is given by $$\theta_\mathrm{tot}({\psi}) \equiv \theta_\mathrm{BG} + \theta_\mathrm{S}({\psi}),$$ the sum of the predicted number of background events $\theta_\mathrm{BG}$ and signal events $\theta_\mathrm{S}({\psi})$. Whilst the predicted number of signal events depends on the model, the background does not. 

Eq.  \[standard\] accounts for statistical fluctuations in the number of observed events for a precisely-known mean $\theta_{\mathrm{tot}}$. However, the predicted mean *itself* may contain a systematic error, due to e.g. an error in the estimate of the effective area of the instrument.  Taking $\epsilon$ to be the (unknown) ratio of the true expected signal contribution $\theta_{\mathrm{S,true}}$ to the nominally-predicted signal contribution $\theta_{\mathrm{S}}$, we have $\theta_{\mathrm{S,true}} = \epsilon \theta_{\mathrm{S}}$. This means that the relative (fractional) error on $\theta_{\mathrm{S}}$ is $\epsilon-1$. To account for this potential systematic error in the number of predicted events, we marginalise over a probability distribution for the relative error $\epsilon-1$ in a semi-Bayesian manner [@Conrad03].  Assuming a log-normal probability distribution for $\epsilon-1$, the likelihood becomes $$\label{number}
{\mathcal{L}}_\mathrm{num}({{n_\mathrm{tot}}}|\theta_\mathrm{tot}) = \frac{1}{\sqrt{2\pi}\sigma_\epsilon}\int_0^\infty \frac{(\theta_\mathrm{BG}+\epsilon\theta_\mathrm{S})^{{{n_\mathrm{tot}}}}e^{-(\theta_\mathrm{BG}+\epsilon\theta_\mathrm{S})}}{{{n_\mathrm{tot}}}!}\frac1\epsilon\exp\left[-\frac{1}{2}\left(\frac{\ln\epsilon}{\sigma_\epsilon}\right)^2\right]\mathrm{d}\epsilon \, .$$ Here $\sigma_\epsilon$ is the fractional systematic error in the predicted number of events.  We take this term to be the sum in quadrature of the fractional uncertainty of the IceCube effective area, and a theoretical error $\tau$. The fractional error in the effective area is itself the sum in quadrature of the fractional statistical and systematic errors determined in deriving the effective area (both have the character of a systematic for the purposes of the likelihood calculation, if not for the effective area calculation itself).  For $\tau$ we use a conservative error of 5%, designed to account for neglected higher-order corrections and possible accumulated numerical round-off errors. 

We apply the systematic uncertainty only to the *signal* prediction in Eq.  \[number\], as neither an uncertainty in the signal prediction from theory, nor an error in the estimation of the effective area, would impact the number of predicted background events.  This is because $\theta_\mathrm{BG}$ is based on the total number of *observed* events away from the Sun.  Eq.  \[number\] gives the complete contribution of the number of events to the total likelihood; we hence refer to this as our adopted ‘number likelihood’. 

We could have adopted a Gaussian probability distribution for $\epsilon-1$ in Eq.  \[number\], but this leads to a non-zero probability that $\epsilon=0$, which is unrealistic.  We have therefore adopted a log-normal distribution for $\epsilon$ in Eq.  \[number\] and throughout this paper, although for small $\sigma_\epsilon$ the difference with the Gaussian case is minimal.  As the Gaussian calculation is significantly faster in some cases, we provide both Gaussian and log-normal routines in . 

Eq.  \[number\] assumes a single, constant systematic error on the effective area.  For real detectors, this error is typically highly energy-dependent.  In general there is no consistent way to allow the prior distribution for $\epsilon$ in Eq.  \[number\] to vary with energy (or any other observable) in an unbinned likelihood calculation.  One option is to model the variation of the systematic error with energy using some parametric form, and then marginalise over the parameters; this however requires assumptions as to the functional form and the degree to which the size of the systematic is correlated at different energies.  Another solution is to bin events into broad ‘energy’ ranges according to their observed $N_\mathrm{chan}$ values, using a Bayesian sorting technique that assumes a prior for the source spectral shape.  Each bin is then assigned a different systematic error for the effective area, based on the range of neutrino energies it covers, and an unbinned analysis is performed on the events *within each bin*.  Ideally, one keeps the number of bins in such a setup to a minimum, to minimise prior-dependence in the final results.  A third, more conservative option, is simply to take $\sigma_\epsilon$ in Eq.  \[number\] to be the largest percentage error seen on the effective area at any energy.  We follow the third strategy; we also attempted an analysis using the second, but found that the additional noise introduced by the need to sort events into the different bins grossly outweighed any advantage gleaned from having a more accurate effective area at higher neutrino energies. 

Spectral and angular likelihoods
--------------------------------

In general $Q({N_i}, {\phi'_i}| {E_i}, {\phi_{i}})$ in Eq.  \[unbinned\_like\] is a function of the spectral and angular resolution of the detector $$\label{sep_instrument}
Q({N_i}, {\phi'_i}| {E_i}, {\phi_{i}},{\psi}) = E_\mathrm{disp}({N_i}|{E_i})PSF({\phi'_i}|{\phi_{i}}), \\$$ where $E_\mathrm{disp}$ is the distribution of $N_\mathrm{chan}$ for incoming neutrinos of energy ${E_i}$ (the energy dispersion of the instrument), and $PSF$ is the detector point-spread function.  Here we have assumed that $E_\mathrm{disp}$ has no angular dependence, and that the energy-dependence of the PSF can be neglected at this stage.  These assumptions are certainly not true in general, given the large width of $E_\mathrm{disp}$ and the strong energy-dependence of the PSF.  Here we apply a hard angular cut around the solar position, reducing the possible range of ${\phi_{i}}$, and therefore the range of variation of $E_\mathrm{disp}$ with angle.  Although this does not make $E_\mathrm{disp}$ entirely isotropic even within the analysis cone, our use of such a cone makes neglecting this angular variation reasonable.  The second approximation is well justified because our per-event analysis already takes the energy-dependence of the PSF into account implicitly, because we apply a unique PSF for each event using the paraboloid sigma (cf.  Section \[icangular\]). 

Similarly, $\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}\,{\mathrm{d}}{\phi_{i}}}({E_i},{\phi_{i}},{\psi})$ in Eq.  \[unbinned\_like\] depends on the true energy spectrum $\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}}$ and spatial distribution $\frac{{\mathrm{d}}P}{{\mathrm{d}}{\phi_{i}}}$ of the signal and background $$\label{sep_theory}
\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}\,{\mathrm{d}}{\phi_{i}}}({E_i},{\phi_{i}}) = \frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}}({E_i},{\psi}) \frac{{\mathrm{d}}P}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}},{\psi}),$$ with $$\begin{aligned}
\label{dpde}
\frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}}({E_i},{\psi}) &= f_\mathrm{BG}\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{E_i}}({E_i}) + f_\mathrm{S}\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{E_i}}({E_i},{\psi})\\
\label{dpdphi}
\frac{{\mathrm{d}}P}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}},{\psi}) &= f_\mathrm{BG}\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}})+f_\mathrm{S}\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}},{\psi}).\end{aligned}$$ Here $\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{E_i}}$, $\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{\phi_{i}}}$, etc refer to the signal and background contributions to the overall predicted spectral and angular distributions, and $$f_\mathrm{S}=f_\mathrm{S}({\psi})\equiv\frac{\theta_\mathrm{S}({\psi})}{\theta_\mathrm{tot}({\psi})}, \hspace{5mm} f_\mathrm{BG}=f_\mathrm{BG}({\psi})\equiv\frac{\theta_\mathrm{BG}}{\theta_\mathrm{tot}({\psi})}$$ are the signal and background fractions, with $\theta_\mathrm{S}$ and $\theta_\mathrm{BG}$ the total predicted number of signal and background events within the analysis cone.  Like $\theta_\mathrm{tot}$ and $\theta_\mathrm{S}$, the fractions $f_\mathrm{S}$ and $f_\mathrm{BG}$ depend upon the model ${\psi}$ by definition, but we typically do not write this dependence out explicitly.[^4]

In Eqs.  \[sep\_theory\]–\[dpdphi\] we have again dropped any explicit energy dependence in the angular part, or angular dependence in the spectral part.  This is perfectly well justified for the background component, as the arrival directions and $N_\mathrm{chan}$ values (and therefore by inference, energies) of background events are observed to be essentially uncorrelated.  For the signal component alone, the assumption that the energy spectrum does not vary with arrival angle is well justified because the IceCube PSF is far larger than the the angular extent of the region where high-energy neutrinos are produced in the Sun, regardless of the energies of the neutrinos.  That the angular extent of the signal does not depend on the energy of the neutrinos from the solar core, nor therefore the model ${\psi}$, follows on exactly the same grounds. ( We can therefore model the angular distribution of the signal as simply a delta function at the solar position.) 

However, these assumptions do not rigorously hold for the *combined* signal and background prediction, as the spatial and spectral characteristics of the signal and background differ; neutrinos coming from the Sun should presumably exhibit a more signal-like spectrum on average than those arriving from elsewhere on the sky.  We deal with this complication implicitly by way of our finite analysis cone.  We evaluate the signal spectral distribution in Eq.  \[dpde\] at precisely the solar position, and then weight its contribution against that of the background by way of the total predicted signal and background fractions $f_\mathrm{S}$ and $f_\mathrm{BG}$. These fractions are the integrated predictions *over the full analysis cone*, meaning that in Eq.  \[dpde\] we essentially derive a mean predicted energy spectrum over the full cone.  Similarly, we implicitly account for the energy-dependence of the angular distribution by weighting the signal and background angular distributions by the total signal and background fractions. 

Collecting all the terms from Eq.  \[unbinned\_like\] where we have retained an explicit dependence on the energy of an incoming neutrino or the number of DOMs it triggers in the detector, we see that the contribution to the overall likelihood from the spectral properties (i.e. ${N_i}$) of event $i$ is $$\begin{aligned}
\label{energyintegral0}
{\mathcal{L}}_{\mathrm{spec},i}({N_i}|{\psi})&\equiv\int_0^\infty E_\mathrm{disp}({N_i}|{E_i}) \frac{{\mathrm{d}}P}{{\mathrm{d}}{E_i}}({E_i},{\psi}) \,{\mathrm{d}}{E_i}\\
\label{energyintegral}
&= \int_0^\infty \left[f_\mathrm{BG}E_\mathrm{disp}({N_i}|{E_i})\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{E_i}}({E_i})+ f_\mathrm{S}E_\mathrm{disp}({N_i}|{E_i})\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{E_i}}({E_i},{\psi}) \right]\,{\mathrm{d}}{E_i}.\end{aligned}$$ The true energy spectrum of the atmospheric background $\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{E_i}}$ would need to be estimated empirically using the observed background flux away from the Sun.  The distribution actually observed by IceCube (Section \[icbackground\]) is the true background distribution $\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{E_i}}$ convolved with $E_\mathrm{disp}$; this is simply the first term in Eq.  \[energyintegral\], so $$\label{speclike}
{\mathcal{L}}_{\mathrm{spec},i}({N_i}|{\psi}) = f_\mathrm{BG}\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{N_i}}({N_i}) +f_\mathrm{S}\int_0^\infty E_\mathrm{disp}({N_i}|{E_i})\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{E_i}}({E_i},{\psi}) \,{\mathrm{d}}{E_i}.$$

Collecting the remaining terms in Eq.  \[unbinned\_like\] with an explicit direction-dependence, the angular likelihood for event $i$ is $$\begin{aligned}
\label{angularintegral0}
{\mathcal{L}}_{\mathrm{ang},i}({\phi'_i}|{\psi})&\equiv\int^\pi_0 PSF({\phi'_i}|{\phi_{i}})\frac{{\mathrm{d}}P}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}},{\psi}) \,{\mathrm{d}}{\phi_{i}}\\
\label{angularintegral}
&=\int^\pi_0 f_\mathrm{BG}\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}})PSF({\phi'_i}|{\phi_{i}})+f_\mathrm{S}\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}},{\psi})PSF({\phi'_i}|{\phi_{i}})\,{\mathrm{d}}{\phi_{i}}.\end{aligned}$$ As for the spectral distribution of the background, we estimate the angular dependence of the background empirically, based on observations away from the Sun (cf.  Section \[icbackground\]).  This provides a direct estimate of the first term in Eq.  \[angularintegral\]. 

The predicted angular distribution of signal events $\frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{\phi_{i}}}$ depends on the actual angular dependence of the signal.  Given that the Sun is a point source to IceCube, we can simply write ${\mathrm{d}}P_\mathrm{S}/{\mathrm{d}}{\phi_{i}}= \delta({\phi_{i}})$.[^5] We can therefore easily solve the integral in the second term in Eq.  \[angularintegral\].  Combining this with the known angular distribution of background events, Eq.  \[angularintegral\] becomes $$\begin{aligned}
\label{anglike}
{\mathcal{L}}_{\mathrm{ang},i}({\phi'_i}|{\psi}) = f_\mathrm{S}PSF({\phi'_i}|{\phi_{i}}=0) + f_\mathrm{BG}\frac{{\mathrm{d}}P_\mathrm{BG}}{{\mathrm{d}}{\phi'_i}}({\phi'_i}).\end{aligned}$$

$PSF({\phi'_i}|{\phi_{i}})$ in Eqs.  \[sep\_instrument\], \[angularintegral0\], \[angularintegral\] and \[anglike\] refers to the unique 1D reduced PSF for each event $i$, as determined by its paraboloid sigma $\sigma_i$. With a full PSF Gaussian in 2D (cf.  Section \[icangular\]) and an analysis cone cut angle $\phi'_\mathrm{cut}$, the reduced PSF is given by $$PSF({\phi'_i}|{\phi_{i}}) = C({\phi_{i}},\sigma_i,\phi'_\mathrm{cut})\frac{|{\phi'_i}-{\phi_{i}}|}{\sigma_i^{\phantom{i}2}}\exp\left[-\frac12\left(\frac{{\phi'_i}-{\phi_{i}}}{\sigma_i}\right)^2\right],$$ where $$C({\phi_{i}},\sigma_i,\phi'_\mathrm{cut})^{-1} = 1 - \mathrm{exp}\left[-\frac12\left(\frac{\phi'_\mathrm{cut}-{\phi_{i}}}{\sigma_i}\right)^2\right]$$ is an angular correction factor that ensures the PSF is unitary within the analysis cone (i.e. the integrated probability for any given event to have come from anywhere at all is 1). 

Total likelihood function
-------------------------

To arrive at the final composite unbinned likelihood, we simply substitute Eqs.  \[number\], \[speclike\] and \[anglike\] into Eq.  \[unbinned\_like\], giving $$\label{full_final_unbinned_like}
{\mathcal{L}}_\mathrm{total}({{n_\mathrm{tot}}},\Xi|{\psi}) = {\mathcal{L}}_\mathrm{num}({{n_\mathrm{tot}}}|{\psi}) \prod_{i=1}^{{{n_\mathrm{tot}}}}{\mathcal{L}}_{\mathrm{ang},i}({\phi'_i}|{\psi})\,{\mathcal{L}}_{\mathrm{spec},i}({N_i}|{\psi}),$$ where $\Xi\equiv\{{N_i},{\phi'_i}\}_{i=1..{{n_\mathrm{tot}}}}$ is the set of all reconstructed event arrival directions and $N_\mathrm{chan}$ values.  This is the likelihood function we employ for all the examples of parameter estimation in this paper. 

Predicted distributions and event counts {#predicted}
----------------------------------------

So far we have focused on the construction of the likelihood function, assuming that one knows how to calculate the number of events and their distribution in energy predicted by a model ${\psi}$. The signal contains contributions from both neutrinos and anti-neutrinos, $$\theta_\mathrm{S}({\psi}) \equiv \theta_{\mathrm{S},\nu}({\psi}) + \theta_{\mathrm{S},{\bar{\nu}}}({\psi}).$$ The energy distribution of the predicted signal is the normalised sum of the neutrino and anti-neutrino flux spectra indicated by model ${\psi}$, each weighted by the respective effective area of the detector ($A_\nu$ or $A_{\bar{\nu}}$).  The sum is further weighted by an overall angular loss factor $L$, giving $$\label{spectral_shape}
\frac{{\mathrm{d}}P_{\mathrm{S}}}{{\mathrm{d}}E}(E,{\psi}) = \frac{t_\mathrm{exp}(\phi=0)}{\theta_{\mathrm{S}}}L(E,\phi'_\mathrm{cut})\left[A_\nu(E)\frac{{\mathrm{d}}\Phi_{\mathrm{S},\nu}}{{\mathrm{d}}E}(E,{\psi})+A_{\bar{\nu}}(E)\frac{{\mathrm{d}}\Phi_{\mathrm{S},\bar{\nu}}}{{\mathrm{d}}E}(E,{\psi})\right].$$ Here $\frac{{\mathrm{d}}\Phi_{\mathrm{S},\nu}}{{\mathrm{d}}E}(E,{\psi})$ and $\frac{{\mathrm{d}}\Phi_{\mathrm{S},\bar{\nu}}}{{\mathrm{d}}E}(E,{\psi})$ are the respective differential neutrino and anti-neutrino fluxes incident on the detector (due to neutrino production in the Sun), with units of events per unit area, energy and time.  The prefactor $t_\mathrm{exp}(\phi=0)/\theta_{\mathrm{S},\nu}$ ensures that Eq.  \[spectral\_shape\] is indeed normalised.  Both $A$ and $L$ are explicitly energy-dependent; $A$ has units of area and $L$ is dimensionless. 

The angular loss factor $L(E,\phi'_\mathrm{cut})$ accounts for events that originate from the Sun, but are so badly reconstructed that they fall outside the analysis cut cone (${\phi'_i}>\phi'_\mathrm{cut}$).  We use the (1D, reduced) mean PSF of the instrument to determine $L$, as defined by the energy-dependent mean angular error $\sigma(E)$ (cf.  Section \[icangular\]).  This gives $$\label{anglossfactor}
L(E,\phi'_\mathrm{cut}) = 1 - \mathrm{exp}\left[-\frac12\left(\frac{{\phi'_\mathrm{cut}}}{\sigma(E)}\right)^2\right].$$

We can simply read off the definition of $\theta_{\mathrm{S},\nu}$ from Eq.  \[spectral\_shape\] by virtue of it being a normalised distribution $$\label{thetaS}
\theta_{\mathrm{S},\nu}\equiv t_\mathrm{exp}(\phi=0)\int_0^\infty L(E)A_\nu(E)\frac{{\mathrm{d}}\Phi_{\mathrm{S},\nu}}{{\mathrm{d}}E}(E,{\psi}) {\mathrm{d}}E.$$ The corresponding expression for $\theta_{\mathrm{S},{\bar{\nu}}}$ is Eq.  \[thetaS\] with $\nu\rightarrow\bar{\nu}$. This integral has no cutoff, as $\lim_{E\to0}A(E)=\lim_{E\to\infty}A(E)=0$.

The total number of predicted background events in the analysis cone ($\theta_\mathrm{BG}$) can be determined by simply rescaling the background rate observed away from the Sun, to the desired exposure time and analysis cone sky fraction. 

$p$-value for model exclusion
-----------------------------

Unfortunately, there is no established method for determining goodness-of-fit using unbinned maximum likelihood estimators, short of a full Neyman construction using Monte Carlo simulations based on Eq.  \[full\_final\_unbinned\_like\].  Whilst this is possible for simple toy models, it is not a computationally feasible option for global fits in virtually any UV-complete theory; lengthy renormalisation-group and relic density calculations make it unrealistic for analyses of supersymmetry, for example. 

A more realistic option is to do model exclusion based entirely on the number of observed events, using Eq.  \[number\].  This has the obvious disadvantage of discarding the spectral and angular information that go into Eq.  \[full\_final\_unbinned\_like\], reducing our overall ability to exclude models compared to a full Neyman construction.  However, this strategy has the distinct advantage that we know exactly how Eq.  \[number\] is distributed in data space, as it simply describes a Poisson process.  This knowledge is what allows us to calculate an absolute $p$-value for any model, completely independent of the value of the likelihood function elsewhere in a theory’s parameter space, and exclude that model with confidence $1-p$.

We begin by identifying each $\psi$ as a specific signal hypothesis, and identifying $\psi$ plus the background distribution inferred from observations away from the Sun as a specific signal+background hypothesis $s+b$. The $p$-value for a hypothesis is defined [@Junk; @James] as the probability, in a repeat of the chosen experiment, of obtaining a test statistic at least as extreme as the one actually observed, if the hypothesis is true.  The test statistic with which we are probing the $s+b$ hypothesis is the likelihood ratio, $$\label{likelihoodratio}
X \equiv \frac{{\mathcal{L}}_\mathrm{num}({{n_\mathrm{tot}}}|\theta_\mathrm{S}+\theta_\mathrm{BG})}{{\mathcal{L}}_\mathrm{num}({{n_\mathrm{tot}}}|\theta_\mathrm{BG})}.$$ $X$ is the ratio of the number likelihood (Eq.  \[number\]) for the signal+background prediction to the number likelihood for a background-only prediction.  A more extreme result is one that gives a lower likelihood ratio (i.e. is less probable if $s+b$ is true).  We do not probe the hypothesis $\psi$ directly with Eqs.  \[number\] and \[likelihoodratio\], as $\psi$ alone does not give a specific prediction for the total number of observed events. 

Referring to Eqs.  \[standard\], \[number\] and \[likelihoodratio\] we see that $X$ is in fact a monotonically increasing function of the number of observed events $n_\mathrm{tot}$. This means that the $p$-value for a hypothesis, as tested using the number likelihood ratio $X$, is the total probability of observing $n_\mathrm{tot}$ or less events if the hypothesis is true.  This is precisely the sum of the likelihood over all possible numbers of observed events less than or equal to $n_\mathrm{tot}$. For the signal plus background hypothesis, this is $$\label{sbpval}
p_{s+b} = \sum_{n=0}^{n_\mathrm{tot}}{\mathcal{L}}_\mathrm{num}\left[n|\theta_\mathrm{tot}\left(\psi\right)\right] \,,$$ and for the background-only hypothesis $$\label{bgpval}
p_{b} = \sum_{n=0}^{n_\mathrm{tot}}{\mathcal{L}}_\mathrm{num}\left(n|\theta_\mathrm{BG}\right) \,.$$ Note that there is a small conceptual complication in the case of the background-only hypothesis: direct substitution of $\theta_\mathrm{S}=0$ into Eq.  \[likelihoodratio\] causes the distribution of the test statistic to collapse to a delta function at $X=1$. Eq.  \[bgpval\] is therefore strictly valid only when Eq.  \[likelihoodratio\] is taken in the limit $\theta_\mathrm{S}\rightarrow 0$, rather than evaluated at exactly $\theta_\mathrm{S}=0$. For further explanation, we refer the interested reader to the original description [@CLs1; @CLs2] of the $CL_s$ technique (of which our $p$-value is essentially an example), and the conditioned confidence interval technique from which the $CL_s$ method was generalised [@Zech]. 

For a true Poisson likelihood given by Eq.  \[standard\], the sums in Eqs.  \[sbpval\] and \[bgpval\] can be performed relatively painlessly, but for the “smeared” Poisson likelihood given by Eq.  \[number\] (the form we actually use[^6]), the sum is over terms each requiring a numerical integration, a potentially slow process.  We provide optimised routines for handling the numerical integration and summing that ensure calculations of the $p$-values are relatively quick. 

In a classical sense $p_{s+b}$ can be regarded as the $p$-value for testing the hypothesis $\psi$ when the background hypothesis is known to be correct.  However, it is sensitive to statistical fluctuations in the actual observed background; a downward fluctuation in the background can cause $\psi$ to be excluded with an unreasonably high confidence level, for example.  To correct for such problems, we use the modified $p$-value as our estimate for the hypothesis $\psi$ alone, defined as [@Junk] $$\label{finalpval}
p_{\psi} = \frac{p_{s+b}}{p_{b}} \, .$$

Example SUSY scans and global fits {#examples}
==================================

As concrete examples of our likelihood formalism in action, we perform three different analyses of the impact of IceCube data on theories of new physics.  We examine two variants of the minimal supersymmetric standard model (MSSM): the constrained MSSM (CMSSM), and a 7-parameter, low-energy phenomenological parameterisation (the MSSM-7).  In the first analysis, we consider the impact of including real IceCube 22-string data on a CMSSM global fit, using the full likelihood function (Eq.  \[full\_final\_unbinned\_like\]).  In this first analysis, as a rough proxy for the full detector we also rescale the effective area of the 22-string analysis, and consider the resulting impact on the CMSSM.  In the second analysis, we assess the ability of our method to recover a hypothetical WIMP annihilation signal from the Sun in terms of the CMSSM, using the same rescaled 22-string detector.  In the third analysis, we give an example of a simple model exclusion exercise in the MSSM-7, using Eq.  \[finalpval\] with the 86-string detector simulation. 

CMSSM global fit with 22-string data {#cmssmglobal}
------------------------------------

The CMSSM is a high-energy parameterisation of the MSSM.  It is defined in terms of the sign of the Higgs mixing parameter $\mu$, the ratio of up-type to down-type Higgs VEVs $\tan\beta$, the universal trilinear coupling $A_0$ and the universal mass parameters $m_0$ and $m_{1/2}$. $\tan\beta$ is defined at the electroweak scale, whereas $A_0$, $m_0$ and $m_{1/2}$ are defined at the scale of grand unification.  The WIMP mass, given by the mass of the lightest neutralino $m_{\chi_1^0}$, is a derived quantity.  We refer the reader to Ref.  [@BaerTata] for full details. 

In order to include IceCube 22-string data in a CMSSM global fit, we coupled  to [@Ruiz06; @Trotta08; @Trotta09] and included Eq.  \[full\_final\_unbinned\_like\] as an additional likelihood component in the global fit.  The total likelihood function we employ for scans in this section includes all components available in the public release of : LEP limits on sparticle and Higgs masses, the anomalous magnetic moment of the muon, limits on rare processes like $b\to s\gamma$, additional $B$-physics observables, the dark matter relic density, and measurements of SM nuisance parameters (refer to [@Trotta08; @Trotta09] for details).  We do not include more recent constraints from the LHC [@Allanach11a; @Allanach11b; @Mastercode12; @MastercodeHiggs; @SuperBayeSHiggs; @Gfitter11], XENON-100 [@MastercodeXENON100; @SuperBayeSXENON100; @Fittino12] or indirect detection [@Scott09c; @Ripken11; @Fittino12], as they are not necessary for our simple illustrative purposes here.  Future analyses using the likelihood formalism we describe, where the focus is on physical interpretation rather than methods, should include such data. 

We performed a scan of the CMSSM parameter space using MultiNest [@MultiNest1; @MultiNest] with 4000 live points and a tolerance of 0.5.  These are settings appropriate for mapping the posterior PDF in a CMSSM scan, but not the profile likelihood [@Akrami09; @SBSpike].  The likelihood formalism we present here is of course equally well adapted for use in frequentist as Bayesian analyses, but obtaining correctly converged profile likelihoods takes roughly an order of magnitude more computing power than posterior PDFs.  We thus rely almost exclusively on posterior maps in this paper, as we are essentially only interested in the results of the scans for the sake of illustration, so the choice of posterior PDF or profile likelihood is fairly arbitrary.  One must also be concerned with correct coverage in presenting profile likelihood results [@Akrami11coverage; @SBcoverage; @Strege12], an additional unnecessary complication for our current purposes. 

An important point however is that the full IceCube likelihood calculation takes a fraction of a second in nearly all cases, so it is not a significant bottleneck in most scans when compared to e.g. a relic density calculation. 

We parsed and plotted the results of our runs using [@pippi].  The parsing step requires binning samples and marginalising posterior PDFs or profiling the likelihood function.  In most cases, we sorted samples into 100 bins in each parameter/observable, and interpolated between those bins for display, at a resolution of 500 points per parameter/observable.  Plots of the predicted number of signal events in IceCube are the exception: in this case we used 130 bins, and interpolated to a resolution of 3000 points per observable direction, to better resolve the region of interest. 

![ Joint 2D posterior probability distributions from a CMSSM global fit including IceCube 22-string event data.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible regions.  Shading and black contours give the result for the fit with IceCube data, whilst grey contours correspond to an identical global fit, but with no IceCube data included.  Existing 22-string data has a minimal impact on the parameter space of the CMSSM.[]{data-label="fig:cmssmIC22"}](IC22_marg1){height="1.05\linewidth"}

![ Joint 2D posterior probability distributions from a CMSSM global fit including IceCube 22-string event data.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible regions.  Shading and black contours give the result for the fit with IceCube data, whilst grey contours correspond to an identical global fit, but with no IceCube data included.  Existing 22-string data has a minimal impact on the parameter space of the CMSSM.[]{data-label="fig:cmssmIC22"}](IC22_marg2){height="1.05\linewidth"}

![ Joint 2D posterior probability distributions from a CMSSM global fit including IceCube 22-string event data.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible regions.  Shading and black contours give the result for the fit with IceCube data, whilst grey contours correspond to an identical global fit, but with no IceCube data included.  Existing 22-string data has a minimal impact on the parameter space of the CMSSM.[]{data-label="fig:cmssmIC22"}](IC22_marg3){height="1.05\linewidth"}

Our first global fit employed the observed events, detector response and effective area of the published 22-string (IC22) analysis [@IceCube09].  For this analysis we used an angular cut of $\phi_\mathrm{cut}=10^\circ$ around the solar position.  The resulting marginalised posterior PDFs are are shown in colour in Fig.  \[fig:cmssmIC22\], along with contours indicating 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible intervals.  Here we show results in the standard $m_0$, $m_{1/2}$ parameter plane, as well as in terms of the neutralino mass – nuclear scattering cross-section planes, for both spin-dependent and spin-independent interactions.  We also plot corresponding credible contours in grey for an identical scan run *without* the inclusion of any IceCube data. 

As expected from the limits produced in the original 22-string analysis [@IceCube09], this data has very little impact upon the preferred regions in a CMSSM fit; the stau co-annihilation region (the small region at low $m_0$ in the leftmost panel and the curved, lower cross-section regions in the other two panels) is unaffected, and any impact on the focus point region (the larger region at high $m_0$ in the leftmost panel and the nearly horizontal, high cross-section bands in the other panels) is difficult to see.  On very close comparison of the black and grey contours in the right panel of Fig.  \[fig:cmssmIC22\] (somewhat easier if one actually refers to the grey contours in the right panel of the *next* figure, Fig.  \[fig:cmssmIC22x100\], which are identical to the grey contours plotted in Fig.  \[fig:cmssmIC22\]), a very small impact can perhaps be seen.  A tiny part of the focus point region, at large spin-dependent scattering cross-section and $m_{\chi_1^0}\sim200$GeV, appears to have been disfavoured by the IC22 likelihood at the level of moving it from the $1\sigma$ credible region to the $2\sigma$ region.  Whether this is simply scanning noise is difficult to say, but the fact that the region in question has a very high posterior PDF in the scan without IceCube data, and is therefore very well-sampled, would seem to argue against such an interpretation. 

To illustrate the approximate effect on the CMSSM of an eventual non-detection of neutrinos from WIMP annihilation in the Sun by the full IceCube detector, using our likelihood formalism, we performed a second global fit with a rescaled IceCube effective area.  We multiplied the IC22 effective area by a factor of 100, and kept all other aspects of the detector as in the 22-string analysis (angular errors, event sample, backgrounds and spectral response); we refer to this as the ‘IC$22\times100$’ configuration.  Although we employ the actual simulated 86-string analysis [@ICRC2011ic86] later for model exclusion, using it for a study such as we describe in this section is not possible, as the 86-string analysis does not contain the requisite $N_\mathrm{chan}$ information to include spectral information in the likelihood function. 

![ As per Fig.  \[fig:cmssmIC22\], but assuming that the detector effective area is boosted by a factor of 100, in order to make a rough estimate of the impact of the complete IceCube detector on the CMSSM.  Grey contours again correspond to the fit with no IceCube data.  Upcoming searches for dark matter with IceCube will robustly exclude the majority of the focus point region of the CMSSM.[]{data-label="fig:cmssmIC22x100"}](IC22x100_marg1){height="1.05\linewidth"}

![ As per Fig.  \[fig:cmssmIC22\], but assuming that the detector effective area is boosted by a factor of 100, in order to make a rough estimate of the impact of the complete IceCube detector on the CMSSM.  Grey contours again correspond to the fit with no IceCube data.  Upcoming searches for dark matter with IceCube will robustly exclude the majority of the focus point region of the CMSSM.[]{data-label="fig:cmssmIC22x100"}](IC22x100_marg2){height="1.05\linewidth"}

![ As per Fig.  \[fig:cmssmIC22\], but assuming that the detector effective area is boosted by a factor of 100, in order to make a rough estimate of the impact of the complete IceCube detector on the CMSSM.  Grey contours again correspond to the fit with no IceCube data.  Upcoming searches for dark matter with IceCube will robustly exclude the majority of the focus point region of the CMSSM.[]{data-label="fig:cmssmIC22x100"}](IC22x100_marg3){height="1.05\linewidth"}

The results of the IC$22\times100$ global fit are shown in Fig.  \[fig:cmssmIC22x100\].  Grey contours again refer to an identical scan performed without the inclusion of any IceCube data.  As is clear from these results, the sensitivity of something resembling the full IceCube detector to both spin-independent and spin-dependent WIMP-nucleon interactions should place very strong constraints on the CMSSM, all but ruling out the majority of the focus-point region.  Whilst this has been shown to already be the case when XENON-100 data is included in the global fit [@MastercodeXENON100; @SuperBayeSXENON100], those constraints are based only on spin-independent scattering, and are therefore particularly sensitive to the adopted prior on the hadronic matrix elements [@Ellis08].  If the value of the pion nuclear sigma term is taken from lattice calculations instead of from experiment, XENON-100 provides almost no constraint on the focus point region.  IceCube should provide a more complete exclusion of most of the focus-point region than XENON-100, because both spin-independent and spin-dependent couplings are expected to contribute to the solar capture rate, and spin-dependent scattering is far less sensitive to hadronic uncertainties than spin-independent scattering.  The full IceCube result will therefore constitute an important independent verification of the XENON-100 exclusion; if IceCube sees a signal in its 86-string configuration, this will be a strong indication that the CMSSM is not responsible for dark matter, or that there is an error in the experimental determination of the pion nuclear sigma term.  Even if there is no signal, IceCube will have a major impact on the parameter space of less constrained versions of supersymmetry, where the spin-independent and spin-dependent nuclear scattering cross-sections are not so tightly coupled as in the CMSSM. 

The appearance of the narrow region at low $m_{1/2}$ and $m_{\chi_1^0}$ (the well-known CMSSM light Higgs funnel region) only in the IC$22\times100$ scan is not surprising, and can be understood via a combination of prior and scanning effects.  This region ‘exists’ as a locale of reasonable fits in all scans, but only appears in the IC$22\times100$ posterior plots because in this case, given the linear prior we have employed on $m_{1/2}$ and $m_0$, removing most of the focus point increases the relative contribution of the funnel to the final model evidence.  These nuances are not directly related to the inclusion of IceCube data, and have been examined extensively already in the literature [@Trotta08; @Akrami09; @SBSpike]. 

a)

b)

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_1){width="\linewidth"}

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_3){width="\linewidth"}

c)

d)

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_2){width="\linewidth"}

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_4){width="\linewidth"}

e)

f)

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_5){width="\linewidth"}

![ 2D posterior probability distributions (a – e) and one comparison 2D profile likelihood (f), from a CMSSM signal recovery based on IceCube 22-string event data only.  We have again assumed that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, to approximate the final detector configuration of IceCube.  Here we have injected a simulated signal corresponding to a 500GeV WIMP annihilating into $W^+W^-$, at a rate that would give 60 signal events on average.  Contours indicate 68.3% ($1\sigma$) and 95.4% ($2\sigma$) credible/confidence regions.  Shading and black contours indicate the recovery using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]), containing number, angular and spectral information.  Grey contours give the recovery achieved with only number and angular information.  The recovery is quite tight, especially with the inclusion of spectral information.[]{data-label="fig:cmssmrecon"}](recon2D_6){width="\linewidth"}

CMSSM recovery validation with mock IceCube data {#cmssmrecon}
------------------------------------------------

To explicitly test the validity and performance of the likelihood construction outlined in Sec.  \[likelihood\], we performed a signal recovery within the CMSSM, using the IC$22\times100$ detector configuration.  We simulated a WIMP annihilation signal in the Sun expected to produce, on average, 60 signal events within the 10 degree angular cut.  We did this as a single realisation of the angular and spectral distribution expected from a 500GeV WIMP annihilating exclusively into $W^+W^-$ pairs, taking into account the angular resolution of the 22-string detector and the expected distribution of $N_\mathrm{chan}$ values from such a signal.  We added this signal to a realisation of the expected background, which we obtained by bootstrap Monte Carlo simulation using the actual all-sky events observed in the 22-string analysis (as described in Sec.  \[icbackground\]).  The particular realisation we employed contained 53 signal and 164 background events within 10 degrees of the solar position. 

We then ran a full parameter scan based on the simulated event list, using  and our modified version of .  For this recovery we included *only* the simulated IceCube data in the likelihood function (i.e. no relic density or other experimental constraints), as we are specifically interested in the ability of IceCube to pin down parameters of a theory of new physics using our likelihood formalism. 

In Fig.  \[fig:cmssmrecon\] we show the results of the recovery in terms of various two-dimensional marginalised posterior PDFs.  Here we give not only joint distributions for $m_0$ and $m_{1/2}$, and for $m_{\chi_1^0}$ and the two nuclear-scattering cross-sections, but also for $m_{\chi_1^0}$ and the velocity-averaged annihilation cross-section $\langle \sigma v\rangle_0$, and $m_{\chi_1^0}$ and the number of predicted signal events in IC$22\times100$. We also give in Fig.  \[fig:cmssmrecon1D\] the 1D marginalised posterior probability distributions of the WIMP mass $m_{\chi_1^0}$, and the three cross-sections.  For comparison, we give the results of a similar recovery *not including the spectral component of the likelihood* in grey. 

Figs.  \[fig:cmssmrecon\] and \[fig:cmssmrecon1D\] show that that our technique accurately recovers both the WIMP mass and expected number of events, and results in rather accurate measurements of the various cross-sections.  Fig.  \[fig:cmssmrecon\]a demonstrates that in this particular example, the fit hones in quite effectively on a part of the focus point region, corresponding to quite small regions in the mass–cross-section planes (Figs.  \[fig:cmssmrecon\]b,c,d).  The performance in the cross-section directions (Figs.  \[fig:cmssmrecon\]b,c,d and \[fig:cmssmrecon1D\]b,c,d) is particularly impressive given that our benchmark point does not specify any cross-sections, only an expected number of events.  A number of different combinations of cross-sections result in an appropriate number of events, as capture and annihilation need not be in equilibrium in the Sun in our calculation, and the relative contributions of spin-independent and spin-dependent scattering to the total capture rate is not specified. 

Although the benchmark point lies just outside the edge of the $1\sigma$ posterior credible region in $m_{\chi_1^0}$ and the predicted number of signal events (Fig.  \[fig:cmssmrecon\]e), this is entirely consistent with realisation noise and the impact of priors.  The profiled likelihood function (Fig.  \[fig:cmssmrecon\]f) in this plane easily encompasses the true value within its $1\sigma$ confidence interval, indicating that the slight preference towards lower masses in the posterior PDF is a prior-driven effect.  As expected, Fig.  \[fig:cmssmrecon\]f is somewhat noisy because the scan was not optimised for profile likelihoods, but the convergence towards a roughly Gaussian likelihood in this plane is obvious.  Recall also that the benchmark point refers to the expected number of total signal events (60), but that the random realisation that we happen to employ contains 53 signal events inside the analysis cone.  It is therefore not at all surprising to find that the best-fit and posterior mean values returned by the scan lie in the range 45–50, given the size of the credible/confidence regions. 

a)

b)

![ 1D posterior probability distributions for the mass of the lightest neutralino and its most important cross-sections, from a CMSSM signal recovery.  This is the same recovery shown in Fig.  \[fig:cmssmrecon\], based on IceCube 22-string event data only, but assuming a factor of 100 boost in the effective area in order to approximate the final detector configuration.  Blue curves give the probability distributions obtained using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]).  Grey curves show the distributions achieved without the use of spectral information.  The use of spectral information significantly improves the accuracy with which the mass and cross-sections can be determined.[]{data-label="fig:cmssmrecon1D"}](recon1D_1){width="\linewidth"}

![ 1D posterior probability distributions for the mass of the lightest neutralino and its most important cross-sections, from a CMSSM signal recovery.  This is the same recovery shown in Fig.  \[fig:cmssmrecon\], based on IceCube 22-string event data only, but assuming a factor of 100 boost in the effective area in order to approximate the final detector configuration.  Blue curves give the probability distributions obtained using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]).  Grey curves show the distributions achieved without the use of spectral information.  The use of spectral information significantly improves the accuracy with which the mass and cross-sections can be determined.[]{data-label="fig:cmssmrecon1D"}](recon1D_3){width="\linewidth"}

c)

d)

![ 1D posterior probability distributions for the mass of the lightest neutralino and its most important cross-sections, from a CMSSM signal recovery.  This is the same recovery shown in Fig.  \[fig:cmssmrecon\], based on IceCube 22-string event data only, but assuming a factor of 100 boost in the effective area in order to approximate the final detector configuration.  Blue curves give the probability distributions obtained using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]).  Grey curves show the distributions achieved without the use of spectral information.  The use of spectral information significantly improves the accuracy with which the mass and cross-sections can be determined.[]{data-label="fig:cmssmrecon1D"}](recon1D_2){width="\linewidth"}

![ 1D posterior probability distributions for the mass of the lightest neutralino and its most important cross-sections, from a CMSSM signal recovery.  This is the same recovery shown in Fig.  \[fig:cmssmrecon\], based on IceCube 22-string event data only, but assuming a factor of 100 boost in the effective area in order to approximate the final detector configuration.  Blue curves give the probability distributions obtained using the full likelihood (Eq.  \[full\_final\_unbinned\_like\]).  Grey curves show the distributions achieved without the use of spectral information.  The use of spectral information significantly improves the accuracy with which the mass and cross-sections can be determined.[]{data-label="fig:cmssmrecon1D"}](recon1D_4){width="\linewidth"}

The grey contours in Fig.  \[fig:cmssmrecon\], and corresponding grey curves in Fig.  \[fig:cmssmrecon1D\], illustrate the dramatic improvement in the mass recovery achieved by including even a crude energy estimator like $N_\mathrm{chan}$. We see that the full likelihood consistently disfavours low masses compared to the recovery without spectral information.  This results in smaller credible regions in all planes of Fig.  \[fig:cmssmrecon\], particularly those involving cross-sections, and tighter credible intervals in Fig.  \[fig:cmssmrecon1D\].  Unsurprisingly, excluding the angular part of the likelihood also has a dramatic effect; using just the number likelihood results in a $2\sigma$ credible region that covers most of the parameter space. 

![ A model exclusion analysis in the MSSM-7, based on Eq.  \[finalpval\] and the consideration of all events seen to arrive from a direction within 20 degrees of the solar position.  Here we use a simulated 86-string IceCube configuration, and simulated data consistent with background.  A number of models will be directly excludable by the full IceCube detector, without any recourse to the rest of the model parameter space, nor any spectral or angular information beyond the number of events within the angular cut cone.  The strongest exclusions will be for models with mixed gaugino/Higgsino dark matter, giving rise to large spin-dependent nuclear scattering rates.  In particular, many of these models have spin-independent cross-sections too low to be detected by large terrestrial experiments like XENON-100.[]{data-label="fig:mssm7"}](MSSM7_gauginofrac){height="1.1\linewidth"}

![ A model exclusion analysis in the MSSM-7, based on Eq.  \[finalpval\] and the consideration of all events seen to arrive from a direction within 20 degrees of the solar position.  Here we use a simulated 86-string IceCube configuration, and simulated data consistent with background.  A number of models will be directly excludable by the full IceCube detector, without any recourse to the rest of the model parameter space, nor any spectral or angular information beyond the number of events within the angular cut cone.  The strongest exclusions will be for models with mixed gaugino/Higgsino dark matter, giving rise to large spin-dependent nuclear scattering rates.  In particular, many of these models have spin-independent cross-sections too low to be detected by large terrestrial experiments like XENON-100.[]{data-label="fig:mssm7"}](MSSM7_SI){height="1.1\linewidth"}

![ A model exclusion analysis in the MSSM-7, based on Eq.  \[finalpval\] and the consideration of all events seen to arrive from a direction within 20 degrees of the solar position.  Here we use a simulated 86-string IceCube configuration, and simulated data consistent with background.  A number of models will be directly excludable by the full IceCube detector, without any recourse to the rest of the model parameter space, nor any spectral or angular information beyond the number of events within the angular cut cone.  The strongest exclusions will be for models with mixed gaugino/Higgsino dark matter, giving rise to large spin-dependent nuclear scattering rates.  In particular, many of these models have spin-independent cross-sections too low to be detected by large terrestrial experiments like XENON-100.[]{data-label="fig:mssm7"}](MSSM7_SD){height="1.1\linewidth"}

MSSM-7 model exclusion with 86-string simulation {#mssm7}
------------------------------------------------

We have also performed scans in the more phenomenological MSSM-7 with seven free parameters at the electroweak scale: the Higgsino mass parameter $\mu$, the gaugino mass parameter $M_2$, the ratio of the Higgs vacuum expectation values $\tan \beta$, a common sfermion mass parameter $m_0$, the CP-odd Higgs boson mass $m_A$ and the trilinear couplings for the third generation squarks $A_b$ and $A_t$. See e.g. [@BergstromGondolo96] for a more detailed description of this model. 

The database we use here is made from several scans focusing on different regions of the parameter space.  The core of the database consists of general random scans allowing the mass parameters to take on values up to 10TeV, $A_t$ and $A_b$ in the range between $-3m_0$ and $3m_0$ and $\tan \beta$ between 1.2 and 60.  We have also performed importance sampling scans with ADSCAN [@Brein:2004kh].  ADSCAN is setup to use an adaptive integration routine (VEGAS) based on importance sampling.  We have then chosen to run ADSCAN with different “importance functions” $G$ focusing on different phenomenologically interesting areas in the parameter space.  The importance function we have mostly used is $$G\left(\Omega_\chi h^2\right) = \exp\left[-\frac12\left(\frac{\Omega_\chi h^2-\Omega_{\rm WMAP}h^2}{\sigma_{\Omega h^2}}\right)^2\right],$$ with $\Omega_{\rm WMAP}h^2 = 0.1099$ from WMAP data [@wmap5].  For the error we have chosen $\sigma_{\Omega h^2}=0.01$ to allow for some theoretical error on top of the experimental one.  We have also chosen importance functions where we add factors to drive the scans into regions of high gaugino fractions, high scattering cross sections or low neutralino masses.  To get a well-sampled database we have also forced ADSCAN into certain mass ranges with importance functions of the type $$G\left(\Omega_\chi h^2, m_{\chi_1^0}\right) = \exp\left[-\frac12\left(\frac{\Omega_\chi h^2-\Omega_{\rm WMAP}h^2}{\sigma_{\Omega h^2}}\right)^2 - \frac12\left(\frac{m_{\chi_1^0}-m_{\rm target}}{\Delta m}\right)^2 \right],$$ varying $m_{\rm target}$ between 10GeV and severalTeV, and commonly using $\Delta m=20$GeV.

In Fig.  \[fig:mssm7\] we show the results of a model exclusion analysis performed on this database for an 86-string IceCube configuration (IC86). 

The detector will be most sensitive to the region where the neutralino is an almost equal Higgsino-gaugino mixture, which is where the spin-dependent scattering cross-section is usually the largest (as can be seen especially in the right panel in Fig.  \[fig:mssm7\]).  This occurs because the spin-dependent cross-section is dominated by $t$-channel $Z$ exchange, which couples to the difference between the two Higgsino contributions; if the neutralino is strongly gaugino, this coupling is absent, whereas if it is strongly Higgsino it becomes equal parts $\tilde{H}_1^0$ and $\tilde{H}_2^0$ and the two contributions cancel. 

As many of these models are captured in the Sun mainly via spin-dependent scattering, the current direct detection experiments (which are mostly sensitive to the spin-independent scattering cross section) do not probe most of the models to which IC86 will be most sensitive.  This can be seen in the middle panel, where many of the models IceCube can probe are far below the sensitivities of both current and near-future direct detection experiments. 

Compared to the standard IceCube ‘hard’ and ‘soft’-channel analyses [@ICRC2011ic86], which are valid for just two specific annihilation final states, the (mixed final state) MSSM-7 models excluded in the analysis here lie (predictably) between the hard and soft sensitivity curves.  The implied sensitivity is closer to the soft-channel curve than it would have been if we had employed a similarly aggressive angular cut as in [@ICRC2011ic86].  For this analysis we used an angular cut of $\phi_\mathrm{cut}=20^\circ$ around the solar position.  In general we have chosen angular cuts so as to include the majority of a potential signal from the Sun.  The cut angle may be optimised for different models however, as higher-mass WIMPs would certainly produce more centrally-localised, high-energy neutrino events.  This is particularly important when doing model exclusion, as the $p$ value Eq.  \[finalpval\] depends only on the expected number of signal and background events, not their angular distributions.  The exclusion power of the full likelihood Eq.  \[full\_final\_unbinned\_like\] is less sensitive to the cut angle, but is still not independent of it, because of the necessary approximations we had to make about the separability of the spectral and angular components.  In principle the cut angle could be profiled or marginalised over as a nuisance parameter in model scans, but this is probably overly conservative; it can be optimised and chosen in advance for specific models based on simulated data, and an approximate scheme (depending, say, on WIMP mass and annihilation branching fractions) constructed for choosing the approximate optimal cut angle for any arbitrary model.  Such an extension may be implemented in future updates to the formalism we have presented in this paper. 

Conclusions
===========

We have constructed an unbinned likelihood formalism for including full event-level IceCube data in parameter explorations of theories for new physics.  The likelihood function includes information about the number, direction and spectral characteristics of neutrino events, and is fast to calculate.  We also constructed a simple associated measure for model exclusion, which is even faster to calculate and can be used for single models, without any reference to other parts of a parameter space. 

We performed a number of example model scans using our likelihood construction within the MSSM.  We carried out global fits to the CMSSM with actual 22-string data, and with the 22-string effective area rescaled to represent the final detector configuration.  Existing 22-string data has little impact on the CMSSM, but the final detector configuration will robustly exclude the majority of the focus point region if IceCube finds no evidence for WIMP annihilation.  We carried out a mock signal recovery in the CMSSM, showing that our method accurately recovers a benchmark point, and constrains model parameters very well.  In the process, we showed the utility of including spectral information in IceCube searches for dark matter.  Finally, we carried out a simple example model-exclusion analysis in the MSSM-7, showing that the 86-string configuration of IceCube will test a number of models that cannot be probed by direct detection experiments in the near future. 

Our likelihood construction is implemented in , and freely accessible to the community.  It will also be made available in a future release of .  The data and simulations we have used for this study, including 22-string IceCube event lists [@IceCube09], are freely available on the web [@filesloc], and in . 

[^1]: [www.darksusy.org](www.darksusy.org)

[^2]: Technically, the 360$^{\circ}$ band includes the direction of the Sun, so the events in this band may include some signal as well as background.  However, the background overwhelmingly dominates any potential signal, so background estimates would be only negligibly affected.  If this was a concern, a section of the band including the Sun can be removed and only an e.g. 300$^{\circ}$ band used. 

[^3]: In both cases the effective areas have been updated slightly with respect to the published values, to properly account for neutrino-antineutrino asymmetries; these changes do not affect previous published IceCube limits. 

[^4]: Note that although not all the actual background is genuinely due to neutrinos (there is a large atmospheric muon component), we can simply interpret the parts of Eqs.  \[sep\_theory\]–\[dpdphi\] due to non-neutrino events as effective neutrino angles and energies corresponding to those tracks.  As we only ever deal with the observed background distributions, not the predicted backgrounds (cf.  Eq.  \[anglike\]), this is entirely valid. 

[^5]: For an extended source, one also needs to take into account any anisotropy in the exposure time $t_\mathrm{exp}(\phi)$ of the detector, defined as the total time over which the detector is sensitive to neutrinos arriving from true angle $\phi$. The effective exposure may vary across the sky due to the seasonal variation of the detector orientation, and choice of analysis cuts.  In this case, the predicted distribution would be $$\label{signal_angular_dist}
    \frac{{\mathrm{d}}P_\mathrm{S}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}}) = \left.\frac{{\mathrm{d}}P_\mathrm{S,i}}{{\mathrm{d}}{\phi_{i}}}({\phi_{i}})t_\mathrm{exp}({\phi_{i}})\middle/\int_0^\pi\frac{{\mathrm{d}}P_\mathrm{S,i}}{{\mathrm{d}}\tilde{\phi_{i}}}(\tilde\phi) t_\mathrm{exp}(\tilde\phi)\,{\mathrm{d}}\tilde\phi\right. ,$$ where ${\mathrm{d}}P_\mathrm{S,i}/{\mathrm{d}}{\phi_{i}}$ is the intrinsic source distribution, and the denominator simply ensures correct normalisation. 

[^6]: Strictly speaking, it is Eq.  \[standard\] that we know the distribution of exactly, not Eq.  \[number\].  Given that Eq.  \[number\] is simply a smeared form of Eq.  \[standard\] however, the two distributions should be very similar.