Update on decaying and annihilating heavy dark matter with the 6-year IceCube HESE data

TL;DR

This paper analyzes 6-year IceCube HESE data to explore if a combination of dark matter decay and astrophysical sources better explains high-energy neutrino observations, finding improved fits and considering dark matter annihilation scenarios.

Contribution

It introduces updated fits for dark matter decay and annihilation channels using 6-year IceCube data, comparing them to astrophysical models and previous 4-year analyses.

Findings

Combined dark matter decay and astrophysical flux improves data fit.

Dark matter decay alone can explain the neutrino data.

Dark matter annihilation scenarios are also consistent with observations.

Abstract

In view of the IceCube's 6-year high-energy starting events (HESE) sample, we revisit the possibility that the updated data may be better explained by a combination of neutrino fluxes from dark matter decay and an isotropic astrophysical power-law than purely by the latter. We find that the combined two-component flux qualitatively improves the fit to the observed data over a purely astrophysical one, and discuss how these updated fits compare against a similar analysis done with the 4-year HESE data. We also update fits involving dark matter decay via multiple channels, without any contribution from the astrophysical flux. We find that a DM-only explanation is not excluded by neutrino data alone. Finally, we also consider the possibility of a signal from dark matter annihilations and perform analogous analyses to the case of decays, commenting on its implications.

Tables

Table 1. Table 1: DM decays (single channel) plus astrophysical power-law flux: Best-fit values for 𝜽 = { τ DM ​ ( N DM ) , m DM , ϕ astro ​ ( N astro ) , γ } 𝜽 subscript 𝜏 DM subscript 𝑁 DM subscript 𝑚 DM subscript italic-ϕ astro subscript 𝑁 astro 𝛾 \bm{\theta}=\{\tau_{\text{DM}}(N_{\text{DM}}),m_{\text{DM}},\phi_{\rm astro}(N_{\rm astro}),\gamma\} , where τ DM subscript 𝜏 DM \tau_{\text{DM}} is expressed in units of 10 28 superscript 10 28 10^{28}  s, m DM subscript 𝑚 DM m_{\text{DM}} in TeV and ϕ astro subscript italic-ϕ astro \phi_{\rm astro} in units of 10 − 18 ​ GeV − 1 ​ cm − 2 ​ s − 1 ​ sr − 1 superscript 10 18 superscript GeV 1 superscript cm 2 superscript s 1 superscript sr 1 10^{-18}~{}{\rm GeV}^{-1}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm sr}^{-1} . The corresponding numbers of DM and astrophysical events are also indicated in parenthesis, as N DM subscript 𝑁 DM N_{\rm DM} and N astro subscript 𝑁 astro N_{\rm astro} . The overall best fit for all those channels is highlighted.

Decay channel	${\tau }_{\mathrm{DM}}\left[{10}^{28}s\right]\left(N_{\text{DM}}\right)$	$m_{\text{DM}}$ [TeV]	${\varphi }_{\mathrm{astro}}\left(N_{\mathrm{astro}}\right)$	$\gamma$
$u\overline{u}$	0.11       (28.4)	1761	0.52   (13.0)	2.34
$b\overline{b}$	0.07       (26.9)	1103	0.58   (14.3)	2.35
$t\overline{t}$	0.11       (28.7)	598	0.45   (12.5)	2.27
${𝑾}^{+}{𝑾}^{-}$	0.37       (28.5)	412	0.47   (12.6)	2.29
$ZZ$	0.43       (27.8)	407	0.52   (13.3)	2.32
$hh$	0.12       (28.8)	611	0.45   (12.6)	2.27
$e^{+}{e}^{-}$	2.20       ( 4.0)	4160	3.53   (37.3)	3.36
${\mu }^{+}{\mu }^{-}$	9.77       ( 4.9)	6583	3.51   (36.5)	3.39
${\tau }^{+}{\tau }^{-}$	0.89       (27.4)	472	0.59   (14.3)	2.36
${\nu }_e{\overline{\nu }}_e$	4.12       ( 3.6)	4062	3.52   (37.7)	3.33
${\nu }_{\mu }{\overline{\nu }}_{\mu }$	4.63       ( 5.0)	4196	3.52   (36.4)	3.41
${\nu }_{\tau }{\overline{\nu }}_{\tau }$	0.96       (16.6)	341	1.58   (24.9)	2.74

Table 2. Table 2: DM-only two-channel decays: Best-fit values for DM → p 1 ​ p ¯ 1 , p 2 ​ p ¯ 2 → DM subscript 𝑝 1 subscript ¯ 𝑝 1 subscript 𝑝 2 subscript ¯ 𝑝 2 \text{DM}\to p_{1}\,\bar{p}_{1},\,p_{2}\,\bar{p}_{2} defined by 𝜽 = { τ DM , m DM , BR } 𝜽 subscript 𝜏 DM subscript 𝑚 DM BR \bm{\theta}=\{\tau_{\text{DM}},m_{\text{DM}},\text{BR}\} , where τ DM subscript 𝜏 DM \tau_{\text{DM}} is expressed in units of 10 27 superscript 10 27 10^{27}  s, m DM subscript 𝑚 DM m_{\text{DM}} in TeV, and BR = Γ DM → p 1 ​ p ¯ 1 / ( Γ DM → p 1 ​ p ¯ 1 + Γ DM → p 2 ​ p ¯ 2 ) BR subscript Γ → DM subscript 𝑝 1 subscript ¯ 𝑝 1 subscript Γ → DM subscript 𝑝 1 subscript ¯ 𝑝 1 subscript Γ → DM subscript 𝑝 2 subscript ¯ 𝑝 2 {\rm BR}=\Gamma_{\text{DM}\to p_{1}\,\bar{p}_{1}}/\left(\Gamma_{\text{DM}\to p_{1}\,\bar{p}_{1}}+\Gamma_{\text{DM}\to p_{2}\,\bar{p}_{2}}\right) . The overall best fit for all those combinations of channels is highlighted.

Decay channels		${\tau }_{\mathrm{DM}}\left[{10}^{27}s\right]$	$m_{\text{DM}}$ [PeV]	BR
$b\overline{b}$	$e^{+}{e}^{-}$	1.82	4.001	0.91
$b\overline{b}$	${\mu }^{+}{\mu }^{-}$	2.00	4.517	0.97
$b\overline{b}$	${\nu }_e{\overline{\nu }}_e$	1.79	3.942	0.97
$b\overline{b}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	1.82	4.015	0.97
$u\overline{u}$	$e^{+}{e}^{-}$	1.76	3.898	0.92
$u\overline{u}$	${\mu }^{+}{\mu }^{-}$	1.85	4.148	0.98
$𝒖\overline{𝒖}$	${𝝂}_{𝒆}{\overline{𝝂}}_{𝒆}$	1.73	3.845	0.97
$u\overline{u}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	1.76	3.906	0.97
$hh$	$e^{+}{e}^{-}$	3.78	4.184	0.92
$hh$	${\mu }^{+}{\mu }^{-}$	4.54	6.132	1.00
$hh$	${\nu }_e{\overline{\nu }}_e$	3.72	4.076	0.97
$hh$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	4.54	6.132	1.00
$W^{+}{W}^{-}$	$e^{+}{e}^{-}$	6.32	4.498	1.00
$W^{+}{W}^{-}$	${\mu }^{+}{\mu }^{-}$	6.32	4.498	1.00
$W^{+}{W}^{-}$	${\nu }_e{\overline{\nu }}_e$	6.32	4.498	1.00
$W^{+}{W}^{-}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	6.32	4.498	1.00

Table 3. Table 3: DM annihilations (single channel) plus astrophysical power-law flux: Best-fit values for 𝜽 = { ⟨ σ ​ v ⟩ , m DM , ξ , ϕ astro , γ } 𝜽 delimited-⟨⟩ 𝜎 𝑣 subscript 𝑚 DM 𝜉 subscript italic-ϕ astro 𝛾 \bm{\theta}=\{\langle\sigma v\rangle,m_{\text{DM}},\xi,\phi_{\rm astro},\gamma\} , where ⟨ σ ​ v ⟩ delimited-⟨⟩ 𝜎 𝑣 \langle\sigma v\rangle  is expressed in units of 10 − 22 ​ cm 3 ​ s − 1 superscript 10 22 superscript cm 3 superscript s 1 10^{-22}~{}{\rm cm}^{3}~{}{\rm s}^{-1} , m DM subscript 𝑚 DM m_{\text{DM}} in TeV and ϕ astro subscript italic-ϕ astro \phi_{\rm astro} in units of 10 − 18 ​ GeV − 1 ​ cm − 2 ​ s − 1 ​ sr − 1 superscript 10 18 superscript GeV 1 superscript cm 2 superscript s 1 superscript sr 1 10^{-18}~{}{\rm GeV}^{-1}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm sr}^{-1} . The corresponding numbers of galactic and extragalactic DM and astrophysical events are also indicated as N DM , G ann superscript subscript 𝑁 DM G ann N_{\rm DM,G}^{\rm ann} , N DM , EG ann superscript subscript 𝑁 DM EG ann N_{\rm DM,EG}^{\rm ann} and N astro subscript 𝑁 astro N_{\rm astro} . The overall best fit for all those channels is highlighted.

Ann. channel	${⟨\sigma v⟩}_{22}$	$m_{\text{DM}}$ [TeV]	$\xi$	${\varphi }_{\mathrm{astro}}$	$\gamma$	$N_{\mathrm{DM},G}^{\mathrm{ann}}$	$N_{\mathrm{DM},\mathrm{EG}}^{\mathrm{ann}}$	$N_{\mathrm{astro}}$
$u\overline{u}$	52.24	260	0.001	1.02	2.52	20.6	0.0	20.2
$b\overline{b}$	24.10	491	0.001	0.81	2.45	23.2	0.0	17.3
$t\overline{t}$	8.20	270	0.001	0.69	2.40	24.8	0.0	15.8
${𝑾}^{+}{𝑾}^{-}$	1.51	178	0.001	0.87	2.48	22.5	0.0	18.1
$ZZ$	1.27	177	0.001	0.91	2.50	22.2	0.0	18.4
$hh$	7.46	278	0.001	0.69	2.40	24.9	0.0	15.8
$e^{+}{e}^{-}$	1.03	159	0.635	1.65	2.75	13.5	1.3	25.8
${\mu }^{+}{\mu }^{-}$	0.63	205	0.001	0.71	2.41	24.6	0.0	15.9
${\tau }^{+}{\tau }^{-}$	0.96	218	0.001	0.66	2.39	25.5	0.0	15.4
${\nu }_e{\overline{\nu }}_e$	0.33	158	3.388	1.67	2.76	10.8	3.8	26.0
${\nu }_{\mu }{\overline{\nu }}_{\mu }$	0.70	159	1.791	0.96	2.52	19.0	3.1	18.9
${\nu }_{\tau }{\overline{\nu }}_{\tau }$	0.70	159	1.945	0.96	2.52	18.8	3.4	18.9

Table 4. Table 4: DM-only two-channel annihilations: Best-fit values for DM ​ DM → a ​ a ¯ , b ​ b ¯ → DM DM 𝑎 ¯ 𝑎 𝑏 ¯ 𝑏 \text{DM}\,\text{DM}\to a\bar{a},\,b\bar{b} defined by 𝜽 = { ⟨ σ ​ v ⟩ , m DM , ξ , BR } 𝜽 delimited-⟨⟩ 𝜎 𝑣 subscript 𝑚 DM 𝜉 BR \bm{\theta}=\{\langle\sigma v\rangle,m_{\text{DM}},\xi,\text{BR}\} , where ⟨ σ ​ v ⟩ delimited-⟨⟩ 𝜎 𝑣 \langle\sigma v\rangle  is expressed in units of 10 − 22 ​ cm 3 ​ s − 1 superscript 10 22 superscript cm 3 superscript s 1 10^{-22}~{}{\rm cm}^{3}~{}{\rm s}^{-1} and m DM subscript 𝑚 DM m_{\text{DM}} in PeV. Galactic and extragalactic DM event numbers corresponding to channel a are indicated as N G a subscript superscript 𝑁 𝑎 G N^{a}_{\rm G} , N EG a subscript superscript 𝑁 𝑎 EG N^{a}_{\rm EG} , respectively. The overall best fit for all those channels is highlighted.

Ann. channels		${⟨\sigma v⟩}_{22}$	$m_{\text{DM}}$ [PeV]	$\xi$	BR	$N_G^a$	$N_{\mathrm{EG}}^a$	$N_G^b$	$N_{\mathrm{EG}}^b$
$𝒖\overline{𝒖}$	$e^{+}{e}^{-}$	18.32	1.99	4.52	0.92	17.9	18.3	1.7	2.8
$u\overline{u}$	${\mu }^{+}{\mu }^{-}$	19.41	2.13	3.24	0.98	19.8	14.9	2.6	3.5
$u\overline{u}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	17.95	1.99	4.15	0.95	18.2	17.0	2.5	3.2
$u\overline{u}$	${\nu }_{\tau }{\overline{\nu }}_{\tau }$	17.42	1.99	4.33	0.95	17.6	17.2	2.6	3.5
$u\overline{u}$	${\tau }^{+}{\tau }^{-}$	17.86	2.16	3.64	0.97	18.1	15.2	3.1	4.3
$u\overline{u}$	$ZZ$	18.07	2.15	3.64	0.93	17.7	14.9	3.6	4.6
$b\overline{b}$	$e^{+}{e}^{-}$	17.14	2.07	5.02	0.92	16.6	19.4	1.7	3.0
$b\overline{b}$	${\mu }^{+}{\mu }^{-}$	19.23	2.40	3.47	0.97	19.1	15.8	2.3	3.3
$b\overline{b}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	17.10	2.06	4.44	0.95	17.2	17.7	2.4	3.4
$b\overline{b}$	${\nu }_{\tau }{\overline{\nu }}_{\tau }$	16.48	2.05	4.66	0.95	16.6	17.9	2.6	3.7
$b\overline{b}$	${\tau }^{+}{\tau }^{-}$	17.91	2.60	4.04	0.97	17.3	16.9	2.5	4.0
$b\overline{b}$	$ZZ$	17.74	2.29	3.77	0.93	17.0	15.2	3.6	4.9
${\mu }^{+}{\mu }^{-}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	0.36	1.06	26.91	0.65	3.0	28.6	0.8	5.7
${\mu }^{+}{\mu }^{-}$	${\nu }_{\tau }{\overline{\nu }}_{\tau }$	0.55	1.78	32.81	0.49	1.9	24.4	1.0	9.7
${\tau }^{+}{\tau }^{-}$	${\nu }_{\mu }{\overline{\nu }}_{\mu }$	0.70	1.90	20.28	1.00	4.5	34.3	0.0	0.0
${\tau }^{+}{\tau }^{-}$	${\nu }_{\tau }{\overline{\nu }}_{\tau }$	0.71	1.89	20.61	0.95	4.4	33.6	0.1	0.7
$W^{+}{W}^{-}$	$e^{+}{e}^{-}$	4.85	2.20	5.09	1.00	13.9	24.8	0.0	0.0
$ZZ$	$e^{+}{e}^{-}$	4.92	2.20	4.64	1.00	14.4	23.9	0.0	0.0

Equations

$${\frac{d{\Phi }_{{\nu }_{\alpha }}}{d{E}_{\nu }}|}_{\oplus }={\varphi }_{\mathrm{astro}}{\left(\frac{E_{\nu }}{100\mathrm{TeV}}\right)}^{-\gamma }+\sum _{\beta }\sum _i|U_{\alpha i}{|}^2|U_{\beta i}{|}^2\frac{d{\Phi }_{\mathrm{DM},{\nu }_{\beta }}}{d{E}_{\nu }},$$

$$\frac{\mathrm{d}{\Phi }_{\mathrm{DM},\mathrm{G}}^{\mathrm{ann}}}{\mathrm{d}{E}_{\nu }}(E_{\nu },b,l)=\frac{⟨\sigma v⟩}{2}\frac{1}{4\pi m_{\mathrm{DM}}^2}\frac{\mathrm{d}N}{\mathrm{d}{E}_{\nu }}{\int }_0^{\infty }{\rho }^2[r(s,b,l)]\mathrm{d}s,$$

$$\frac{\mathrm{d}{\Phi }_{\mathrm{DM},\mathrm{EG}}^{\mathrm{ann}}}{\mathrm{d}{E}_{\nu }}(E_{\nu })=\frac{⟨\sigma v⟩}{2}\frac{{\Omega }_{\mathrm{DM}}^2{\rho }_c^2}{4\pi m_{\text{DM}}^2}{\int }_0^{\infty }\frac{dz}{H(z)}(1+z{)}^3\zeta (z)\frac{\mathrm{d}N}{\mathrm{d}{E}_{\nu }}\left[\left(1+z\right)E_{\nu }\right],$$

$$\frac{\mathrm{d}{\Phi }_{\mathrm{DM}}^{\mathrm{ann}}}{\mathrm{d}{E}_{\nu }}(E_{\nu })=\frac{\mathrm{d}{\Phi }_{\mathrm{DM},\mathrm{G}}^{\mathrm{ann}}}{\mathrm{d}{E}_{\nu }}(E_{\nu })+\xi \frac{\mathrm{d}{\Phi }_{\mathrm{DM},\mathrm{EG}}^{\mathrm{ann}}}{\mathrm{d}{E}_{\nu }}(E_{\nu }).$$

$$⟨\sigma v⟩\lesssim \frac{4\pi }{m_{\text{DM}}^2{v}_{\mathrm{local}}}\simeq 1.5\times 10^{-23}{\mathrm{cm}}^3/\mathrm{s}{\left(\frac{100\mathrm{TeV}}{m_{\text{DM}}}\right)}^2\left(\frac{10^{-3}c}{v_{\mathrm{local}}}\right),$$

Full text

††preprint: IFIC/19-19

1 Introduction

Since the discovery of the first astrophysical high-energy neutrinos at the IceCube neutrino telescope Aartsen:2013jdh , more neutrinos have been steadily detected Aartsen:2014gkd ; Aartsen:2015rwa ; Kopper:2015vzf ; Aartsen:2016xlq ; Kopper:2017zzm ; Haack:2017dxi ; Wandkowsky:2018 . This has opened the era of neutrino astronomy and understanding the origin of this high-energy neutrino flux is currently one of the most important problems in astroparticle physics Laha:2013lka ; Halzen:2013dva ; Waxman:2013zda ; Anchordoqui:2013dnh ; Murase:2014tsa ; Meszaros:2017fcs ; Ahlers:2018fkn .

From the observational point of view, different types of event samples have been considered and a broad picture of the astrophysical neutrino flux has emerged. It can be broadly summarized as follows. The high-energy starting events (HESE), with electromagnetic (EM)-equivalent deposited energies above 10 TeV and produced inside the detector, can be explained as produced by a soft neutrino spectrum with a power-law index ($E_{\nu}^{-\gamma}$) close to $\gamma\simeq 2.9$ Kopper:2017zzm ; Wandkowsky:2018 . Neutrinos of all flavors contribute to this sample and two topologies have been discriminated so far: muon tracks, mainly produced from muon neutrinos, and showers, produced by all flavors. On the other hand, through-going muon events, produced by up-going muon neutrinos (i.e, from the Northern hemisphere), are sensitive to energies above $\sim$ 200 TeV and can be understood in terms of a harder spectrum, $\gamma\sim 2.2$ Aartsen:2015rwa ; Aartsen:2016xlq ; Haack:2017dxi . This tension has motivated the proposal of two-component fluxes Chen:2014gxa ; Aartsen:2015knd ; Palladino:2016zoe ; Vincent:2016nut ; Palladino:2016xsy ; Palladino:2017qda , not necessarily isotropic and with a softer spectrum for $E_{\nu}\lesssim 100$ TeV, or the possible existence of a break or softening in the spectrum Winter:2014pya ; Anchordoqui:2014hua ; Palomares-Ruiz:2015mka ; Aartsen:2015zva ; Anchordoqui:2016ewn .

On another hand, based on the near isotropy of the observed events Aartsen:2013jdh ; Ahlers:2013xia ; Aartsen:2014gkd ; Anchordoqui:2014rca ; Troitsky:2015cnk ; Kistler:2015oae ; Palladino:2016zoe ; Vincent:2016nut ; Denton:2017csz ; Albert:2017oba ; Aartsen:2017ujz ; Krings:2017pif , extragalactic sources are favored. Only a relatively small fraction of the neutrino flux seems to be compatible with a galactic disc origin Denton:2017csz ; Albert:2017oba ; Aartsen:2017ujz ; Krings:2017pif . Note, however, that the limits on the galactic contribution are usually obtained using spatial and energy templates for diffuse neutrino emission, so different conclusions might be reached for other distributions Pagliaroli:2017fse . Thus, currently a galactic contribution to the total flux cannot be excluded, especially if the galactic contribution comes from the halo. Interestingly, the recent observation of the hardening of the $\gamma$-ray flux above $\sim$ 300 GeV in the Fermi-LAT diffuse data at large galactic latitudes seems to point to the existence of a correlated galactic neutrino flux with a hard spectrum Neronov:2018ibl (see also Ref. Neronov:2016bnp ). An example of such a spectrum could be the one produced from the decays of heavy dark matter (DM) particles. In this case, a significant fraction of the IceCube neutrinos would come from DM decaying in the halo of our galaxy, while the total neutrino flux from DM decays would result from the sum of the galactic and extragalactic (isotropic) contributions. The consistency, or even mild preference, of the angular distribution of the IceCube events with such a scenario has already been discussed Bai:2013nga ; Esmaili:2014rma ; Troitsky:2015cnk ; Chianese:2016opp , without a conclusive answer yet.

The possibility to explain part or all of the neutrino event spectrum observed by IceCube in terms of heavy DM decays or annihilations has been extensively considered Feldstein:2013kka ; Esmaili:2013gha ; Bai:2013nga ; Ema:2013nda ; Bhattacharya:2014vwa ; Zavala:2014dla ; Higaki:2014dwa ; Ema:2014ufa ; Rott:2014kfa ; Esmaili:2014rma ; Fong:2014bsa ; Daikoku:2015vsa ; Murase:2015gea ; Esmaili:2015xpa ; Aisati:2015vma ; Roland:2015yoa ; Anchordoqui:2015lqa ; Boucenna:2015tra ; Ko:2015nma ; Troitsky:2015cnk ; EsmailiTaklimi:2016bbx ; Esmaili:2016swq ; Chianese:2016opp ; Dev:2016uxj ; Fiorentin:2016avj ; Dev:2016qbd ; DiBari:2016guw ; Chianese:2016smc ; Chianese:2016kpu ; Kuznetsov:2016fjt ; Cohen:2016uyg ; Borah:2017xgm ; Hiroshima:2017hmy ; Bhattacharya:2017jaw ; Chakravarty:2017hcy ; Chianese:2017nwe ; Dhuria:2017ihq ; Aartsen:2018mxl ; Sui:2018bbh ; Chianese:2018ijk ; Blanco:2018esa . Recently, the IceCube collaboration published two analyses using independent datasets Aartsen:2018mxl : 6-year muon data Aartsen:2016xlq and all-sky cascades from the 2-year medium energy starting events (MESE) Aartsen:2014muf . By only considering hard channels ($H+\nu$ and $Z+\nu$), a non-zero DM decay contribution was found for the best fits, including an astrophysical (isotropic) power-law flux. Nevertheless, this was not statistically significant and only limits on the lifetime were derived for these channels as well as for other ones, although in those other cases only using the cascade analysis. On the other hand, the first DM decay analysis using the 6-year HESE data Kopper:2017zzm was performed for fixed values of the spectral index of the astrophysical flux ($\gamma=2.0,\,2.2$) Chianese:2017nwe , and a best fit for DM masses of $\sim 400-500$ TeV was found. Another recent work Sui:2018bbh made use of the 6-year HESE data Kopper:2017zzm and additionally of the 8-year through-going muon data Haack:2017dxi , which dominates the statistics. This analysis, however, considered only a particular DM decay model, with a final hard neutrino spectrum. It was also found that adding a DM decay contribution, with a DM mass $\sim 300$ TeV, to an astrophysical power-law flux, provides a better fit to data.

In this work we consider just the 6-year HESE data above 60 TeV in EM-equivalent deposited energy,111Note that preliminary results using 7.5-year HESE data have already been presented Wandkowsky:2018 , including limits on the DM contribution Yuan:2018 , although details about the charge recalibration and the improved reconstructions are not public yet. and evaluate whether a DM contribution could provide a good fit to the HESE data and reduce by itself the tension with the through-going muon data. This is an update of our previous detailed analysis, which used the 4-year HESE data Bhattacharya:2017jaw . Thus, one of the main questions is whether our previous conclusions based on 4-year data are substantially altered when additional two years of data are included, which we will discuss thoroughly in this paper. In particular, we also consider the case of DM decaying into two channels (leaving the relative branching ratio as a free parameter), so that the entire HESE data could be explained without the need for a power-law astrophysical flux. As previously found Esmaili:2013gha ; Esmaili:2014rma ; Bhattacharya:2017jaw , we confirm that current HESE data could be interpreted in terms of this scenario.222Nevertheless, it could be in tension with gamma-ray data, as we mention below. This may seem to be at odds with the claim in Ref. Aartsen:2018mxl , which concluded that a DM-only scenario is disfavored, by neutrino data, with respect to a DM plus a power-law flux. However, apart from the fact that the referred tension is below the 1$\sigma$ confidence level (CL), that conclusion is based on a fixed combination of hard channels for DM decay.

Additionally, we also consider the possibility to explain the HESE data, in part or entirely, in terms of neutrinos from DM annihilations. Due to the dependence on the square of the DM density, this scenario could present a sharper enhancement towards the galactic center. To study this, we leave the normalization of the cosmological DM clumping factor as a free parameter, which modifies the relative importance of the galactic and extragalactic contributions and thus, both the spectral and the spatial distributions of the neutrino flux. Although it is well known that this requires very large annihilation cross sections Feldstein:2013kka ; Zavala:2014dla ; Esmaili:2014rma ; Chianese:2016opp ; Chianese:2016kpu , these data could, in principle, probe regions below the unitarity bound Griest:1989wd ; Hui:2001wy . We will see, however, that this is not the case with the current HESE data.

The main consideration imposed on our analysis due to the addition of two years worth of additional data, apart from the obvious improvement in statistics, comes from the proliferation of sub-PeV events without a balancing increase in the PeV events. The tension between data and uniform power-law fits, already exaggerated due to this tilting of the spectrum in favor of sub-PeV events, has been quantified as a $2.6\sigma\ (2.1\sigma)$ bump in events between EM-equivalent deposited energies of $60-100$ TeV when compared against a flat power-law spectrum with spectral index $\alpha=2.0\ (2.2)$ Chianese:2017nwe .

In view of these qualitative changes to the HESE data, we explore the nature of fits from combined DM (decays or annihilations) and power-law astrophysical neutrino fluxes, specifically looking at how these affect our previous 4-year HESE analysis.

The paper is organized as follows: in section 2 we briefly introduce the HESE sample we use and describe the main ingredients of the analysis we perform. In section 3 we present the results for the case of DM decays and in section 4 the corresponding ones for the annihilating scenario. Finally, in section 5 we draw our conclusions.

2 Analysis of the 6-year HESE data

The 6-year HESE data we use in this work correspond to 2078 days of runtime during which, 82 events have been observed (28 new events in the last two years) with EM-equivalent deposited energies above 10 TeV. Two of these events, however, were produced by a coincident pair of background muons from unrelated cosmic-ray air showers and are not included in our analysis. In this work we consider the EM-equivalent deposited energy interval [60 TeV–10 PeV], which contains 50 events (18 from the last two years), classified as either tracks or showers.

At these energies, cosmic-ray secondaries produced in the atmosphere, muons and neutrinos, are the main sources of background. As done in our previous analysis Bhattacharya:2017jaw , here we do not consider the potential contribution from prompt atmospheric neutrinos from charm decays. For energies above 60 TeV, the number of background atmospheric muon and neutrino events for 2078 days is taken as:333For the muon background we use the IceCube numbers, whereas for the neutrino background we scale it with respect to the 4-year data Kopper:2015vzf . We do this because in the 6-year analysis slightly larger cross sections have been used Kopper:2017zzm and we have not corrected for this. The difference, though, is only about 0.6 events. $N_{\mu}=1.3$ and $N_{\nu}=5.1$.

The two additional contributions we consider in this work for the neutrino flux at Earth are: an isotropic power-law flux and the flux from DM decays or annihilations,

[TABLE]

where $\gamma$ is the spectral index and $\phi_{\rm astro}$ is the astrophysical flux normalization (per flavor). Here we assume the canonical flavor composition from hadronic sources after averaged oscillations, $(1:1:1)_{\oplus}$, for both neutrinos and antineutrinos, and hence, $\phi_{\rm astro}$ is flavor-independent. $U$ is the neutrino mixing matrix and the above product represents the oscillation probability of flavor neutrino $\nu_{\alpha}$ into $\nu_{\beta}$, after all terms dependent on mass-squared differences are averaged out. The term $\Phi_{{\rm DM}}$ corresponds to the flux of neutrinos from DM decays or annihilations at production. It depends on the DM mass, $m_{\text{DM}}$, on the decay/annihilation channel(s), and for decays, on the lifetime, $\tau_{\text{DM}}$, and for annihilations, on the annihilation cross section, $\langle\sigma v\rangle$, and on the clumping factor (see below). In both cases, decays and annihilations, the neutrino flux has two contributions: an extragalactic component, which originates from DM particles in halos at all redshifts, and thus it is isotropic; and a galactic contribution, which comes from DM in our galactic halo and it is anisotropic. Whereas the DM decay flux depends linearly on the DM density, the flux from DM annihilations depends on its square. Therefore, the galactic contribution is more sharply peaked towards the galactic center in the case of annihilations. For the details of our computation of the neutrino flux from DM decays we refer the reader to our previous work Bhattacharya:2017jaw . The flux from DM annihilations is described below.

Analogously to Ref. Bhattacharya:2017jaw , to which we refer the reader for details, we perform unbinned extended maximum likelihood analyses using the EM-equivalent deposited energy, the event topology, and the hemisphere of origin of each of the 50 events in the 6-year HESE sample within [60 TeV–10 PeV]. For the description of the computation of the event spectra of both signal and background, we refer the reader to Refs. Palomares-Ruiz:2015mka ; Vincent:2016nut .

3 DM decays

The neutrino flux from DM decays depends on the DM lifetime, $\tau_{\text{DM}}$, the mass, $m_{\text{DM}}$, and the decay channels. In this work we only consider two-body final state channels and perform two types of fits:

1) involving a single channel DM decay plus an isotropic power-law flux and

2) where the entire spectrum is explained purely by neutrinos from DM decay via two channels.

In the former case, the total neutrino flux for a given decay channel is described in terms of the set of free parameters $\bm{\theta}=\{\tau_{\text{DM}},m_{\text{DM}},\phi_{\rm astro},\gamma\}$, while in the latter scenario, for each pair of channels, it is described in terms of $\bm{\theta}=\{\tau_{\text{DM}},m_{\text{DM}},{\rm BR}\}$, where ${\rm BR}$ is the branching ratio for decay into one of the two channels.

In the case of decays, the relative contribution between the galactic and the extragalactic components is determined by the spatial distribution of the galactic density. The extragalactic component is isotropic and does not depend on the DM clustering properties across redshifts.

3.1 Results: DM decays plus isotropic astrophysical power-law flux

The best-fit parameters corresponding to a total signal flux comprising neutrinos from a power-law astrophysical spectrum determined by $\phi_{\rm astro}$ and $\gamma$, and from DM decays, as a function of $\tau_{\text{DM}}$ and $m_{\text{DM}}$, for different two-body decay channels, are indicated in Table 1. The comparison of the best-fit likelihoods for all channels is shown in Figure 1. On top of each bar, we indicate the best fit obtained for the astrophysical index $\gamma$. In Figure 2 we show the event spectra for two channels: in one, DM decay explains the low-energy data ($\text{DM}\to W^{+}W^{-}$) while in the other it explains the PeV events ($\text{DM}\to\nu_{e}\bar{\nu}_{e}$).

The overall best fit comes from decays to bosons $W,\ Z,\text{ and }h$ (also top quarks), which prefer $m_{\text{DM}}\sim(400-650)$ TeV, thereby contributing only to the sub-PeV spectrum. In these cases, a relatively flat astrophysical flux, with $\gamma\sim 2.3$, accounts for the high-energy events (left panel of Figure 2) and supplements the contribution from DM decays to sub-PeV energies therefore, making the total sub-PeV spectrum steeply falling.

Fits to hard channels with DM mass in the high PeV (e.g., decays to charged lepton pairs and neutrinos) do comparatively worse, yet the statistical significance is weak. In these cases, for the best fits (except $\text{DM}\to\tau^{+}\tau^{-}$), events from DM decays explain the higher end of the spectrum, while a soft astrophysical power-law flux fills in the sub-PeV part of the event spectrum, falling quickly to have no impact on the high energy events (right panel of Figure 2).444From the comparison of the two panels of Figure 2, it might seem that the fit for $\text{DM}\to\nu_{e}\bar{\nu}_{e}$ is slightly better. However, a note of caution is in order. The event spectra is shown in bins, so a particular choice of binning can have a misleading impact. We avoid this problem we perform an unbinned extended maximum likelihood analysis to determine the best fits. Moreover, important parameters that influence the fits, such as the type of event topology or the neutrino direction, cannot be illustrated in this event-rate plot. As it turns out, the full likelihood analysis reveals the overall best fit comes from decays to $W^{+}W^{-}$ and not from those to $\nu_{e}\bar{\nu}_{e}$ (Figure 1). As in the 4-year HESE analysis Bhattacharya:2017jaw , the lack of high-energy tracks in the data influences the fit, and the flux from $\text{DM}\to\mu^{+}\mu^{-}$ turns out to be slightly disfavored compared to cases with other hard channels.

Given that the best-fit spectral index for several channels — indeed for those giving better fits — naturally turns out to be within $1\sigma$ of the best fit obtained by the latest IceCube through-going muon analysis, $\gamma_{\rm TG}=2.19$ Haack:2017dxi , we forego the prior-added analysis we carried out in our previous work. This is one of the major takeaways from the current analysis: some specific combinations of DM decay and astrophysical flux naturally lead to best-fit spectral indices consistent with that from the IceCube through-going muon analysis, sensitive only to the highest energy neutrinos. This is in contrast to what occurred with the analysis of the 4-year HESE data Bhattacharya:2017jaw , where we obtained a best fit for a hard channel ($\text{DM}\to e^{+}e^{-}$) with PeV mass, explaining the high-energy part of the event spectrum and thus, with a very soft astrophysical spectrum to describe the low-energy side.

3.1.1 Parameter correlations and preferred regions

Here we discuss the correlations between the parameters and compute the preferred region of parameter space to the statistical 1$\sigma$ and 2$\sigma$ CL against the corresponding best-fit point (Figures 3 and 4). Correlations between two parameters corresponding to different flux components (e.g., between $\tau_{\text{DM}}$ and $\gamma$) demonstrate the careful balancing act performed by these two complementary fluxes while fitting to the data. On the other hand, correlations between parameters representing the same flux (e.g., $\tau_{\text{DM}}$ and $m_{\text{DM}}$) reflect the sensitivity of the corresponding flux component to the best fit.

As noted above, for a decay channel like $\text{DM}\to W^{+}W^{-}$, the majority of the total signal events comes from DM decays, with all of these events populating the sub-PeV energies. The flatter astrophysical spectrum then fills in the missing events, with most of them lying at high energies. Correlation plots in Figure 3 are indicative of this nature. In particular, the correlation between $m_{\text{DM}}$ and $\tau_{\text{DM}}$ on one hand and the astrophysical spectral index $\gamma$ on the other shows the preference for a relatively flat spectrum ($1\sigma$ CL region lies below $\gamma=3$) and a large DM contribution (close to the lower edge of the illustrated lifetime range) at low energies (narrow $1\sigma$ CL region around low $m_{\text{DM}}$). It is interesting to compare these results to those obtained with the 4-year HESE data Bhattacharya:2017jaw . Then, there were two $1\sigma$ CL regions in the $(m_{\text{DM}},\tau_{\text{DM}})$ parameter space, but the best fit was in the high-mass one. Now, with the additional events detected around $\sim 100$ TeV, only the region around $m_{\text{DM}}\approx 400$ TeV remains, which substantially changes the event spectrum corresponding to the best fit for this channel.

For channels with harder spectra (e.g., $\text{DM}\to\nu_{e}\bar{\nu}_{e}$), multiple near-degenerate $1\sigma$ regions in terms of $m_{\text{DM}}$ open up (Figure 4). While the best fit for this channel lies in the high-mass region, comparative $1\sigma$ CL regions are also allowed in the low-$m_{\text{DM}}$ region. This is in contrast to what happened with the 4-year HESE data Bhattacharya:2017jaw , where only the high-mass region was preferred at $1\sigma$ CL. Again, this is a consequence of the increase in low-energy events in the 6-year HESE data vis-à-vis the 4-year sample, while leaving the PeV spectrum unchanged. Clearly, the $1\sigma$ CL contour corresponding to low-$m_{\text{DM}}$ allows for a flatter astrophysical flux, while the allowed $1\sigma$ CL contour corresponding to high-$m_{\text{DM}}$ extends to the high-$\gamma$ region, indicating a steeply falling astrophysical flux. It should be noted that for these hard channels, the constraints on the astrophysical flux parameters are more restrictive than in the case of $\text{DM}\to W^{+}W^{-}$.

3.1.2 Limits on the DM lifetime

For each decay channel, we also estimate the maximum allowed contribution from DM decays and thus, obtain limits on the corresponding values of $\tau_{\text{DM}}$ as a function to $m_{\text{DM}}$. The 95% CL limits on the DM lifetime and on the number of DM events ($N_{\text{DM}}$), for all single-channel DM decays in Table 1, are shown in Figure 5 as a function of the $m_{\text{DM}}$. These results are depicted along the bounds obtained from $\gamma$-ray observations in Ref. Cohen:2016uyg .555Note that a factor of a few stronger bounds have been recently obtained Blanco:2018esa .

In comparison with analogous limits obtained in our previous analysis using 4-year HESE data Bhattacharya:2017jaw , they are not significantly altered by the two years of additional data except at the low-$m_{\text{DM}}$ region for neutrino and charged lepton channels, where the bounds strengthen by about an order of magnitude at $m_{\text{DM}}$ values between 100–200 TeV.

Since the fits to the gauge boson channels shift toward lower values of $m_{\text{DM}}$, the corresponding best fit is somewhat in tension with the gamma-ray bounds on DM lifetime obtained in Ref. Cohen:2016uyg . On the other hand, both because these bounds are weaker for leptons and because they weaken with increasing $m_{\text{DM}}$, best fits for decays to charged leptons and neutrinos with $m_{\text{DM}}$ in the PeV region can evade these bounds. As already happened with the 4-year HESE data Bhattacharya:2017jaw , decays to light quark pairs are in strong tension with gamma-ray data: neutrino limits are weaker than the corresponding gamma-ray bounds, and the best fits run afoul of them by an order of magnitude. In all other cases, neutrino limits are stronger than the corresponding gamma-ray bounds in the high-$m_{\text{DM}}$ regions.

When understanding the tension between the best-fit DM parameters and gamma-ray constraints, a note of caution is in order. The gamma-ray bounds on DM decays obtained in Ref. Cohen:2016uyg , and recently reevaluated and slightly strengthened in Ref. Blanco:2018esa , are obtained from the DM contribution to the diffuse gamma-ray background666Also known as isotropic gamma-ray background (IGRB). (DGRB) observed by Fermi-LAT Ackermann:2014usa . A heavy decaying DM produces high-energy gamma-rays that initiate an electromagnetic cascade en route to the Earth and contributes to the DGRB in the energy range $\lesssim$ a few 100 GeV. On the other hand, low-energy astrophysical contributions to the DGRB could in fact entirely explain the observed DGRB Fornasa:2015qua .777The theoretical expectation for the contribution from these sources at some energies can completely explain the DGRB data, or even overshoot it. Obviously, by accounting for these contributions from astrophysical sources to the DGRB, as is done in Refs. Cohen:2016uyg ; Blanco:2018esa , there remains very little room for any contribution from the high-energy cascaded flux. Nevertheless, the same approach would strongly constrain any type of source of high-energy neutrinos, including the conventional $pp$ and $p\gamma$ astrophysical sources, since in these scenarios, an accompanying high-energy gamma-ray flux would also be produced at the source and initiate an electromagnetic cascade (except if the source is opaque to high-energy gamma-rays). In particular, the interpretation of the HESE data in terms of a single power-law flux with spectral index $\gamma\simeq 2.9$ is in strong tension with the DGRB (see Ref. Murase:2013rfa , where an upper limit of $\gamma\lesssim 2.2$ was derived for $pp$ sources by considering the whole DGRB data; using the DGRB data after subtracting low-energy astrophysical source contributions would lead to even stronger limits). When taking the more conservative gamma-ray limits derived from the whole DGRB data Murase:2015gea , and not the low-energy source-subtracted one, all the best fit points in Figure 5 are actually allowed.

3.2 Results: DM decays via multiple channels

We also analyze the consequences of turning off the astrophysical flux entirely and instead allowing DM to decay via multiple channels. This analysis proceeds in the same vein as the previous analysis with the 4-year HESE data (see Sec. 6 in Ref. Bhattacharya:2017jaw ). Allowing the DM to decay via two distinct channels allows us to parameterize the resulting flux in terms of the following physical quantities:

a) the DM mass, $m_{\text{DM}}$,

b) the DM lifetime $\tau_{\text{DM}}$, and

c) the branching ratio in favor of the (arbitrarily ordered) first channel: ${\rm BR}=\Gamma_{\text{DM}\to p_{1}\,\bar{p}_{1}}/\left(\Gamma_{\text{DM}\to p_{1}\,\bar{p}_{1}}+\Gamma_{\text{DM}\to p_{2}\,\bar{p}_{2}}\right)$.

In this scenario, the DM mass is necessarily pushed to PeV values to accommodate both PeV and sub-PeV events. Results for a few selected combinations of decay channels are shown in Table 2 and the $\Delta\chi^{2}$ for each combination with respect to the best fit is shown in Figure 6. The overall best fits come from the combination of $u\bar{u}$ and $e^{+}e^{-}$ or $\nu_{e}\bar{\nu}_{e}$ channels with the DM mass around 4 PeV and a 92–97% branching ratio in favor of decays to the quark (soft) channel. The hard-channel decay, either to $\nu_{e}\bar{\nu}_{e}$ or to $e^{+}e^{-}$, explains the PeV events while the softer secondary neutrino flux from decays to $u\bar{u}$ explains the sub-PeV events. Other combinations of soft and hard channels provide similar results. Note, however, that if both channels are hard (e.g., $W^{+}W^{-}$ and $\nu_{e}\bar{\nu}_{e}$), the fit worsens and prefers the spectrum to be dominated by the softer of the two spectra.

Qualitatively, these patterns are similar to the results obtained with the 4-year HESE data Bhattacharya:2017jaw , albeit with slightly larger branching ratios into soft channels. In relative terms, when compared to the DM single channel decay plus astrophysical flux, we find that the overall best-fit two-channel combination does somewhat worse with the updated 6-year HESE data. Yet, this is not statistically significant.888The DM-only scenario here considered fits the data using three degrees of freedom, as opposed to four in the scenario with DM plus astrophysical fluxes, considered in the previous section. Consequently, a straight comparison to determine which of these two scenarios provides a comparatively better fit is not straightforward and is out of the scope of this work. The main reason for this slight swing is due to the fact that the 6-year HESE data strongly prefers DM decay fluxes with low-$m_{\text{DM}}$, and the requirement in this scenario of having to explain PeV events also from DM decay necessitates the addition of at least one hard component. This ends up worsening the fit slightly in comparison to the two-channel best fit obtained in the previous work. Nevertheless, unlike Ref. Aartsen:2018mxl , we cannot conclude that the DM-only scenario is disfavored by neutrino data. This is in agreement with previous findings using earlier HESE data sets Esmaili:2013gha ; Esmaili:2014rma ; Bhattacharya:2017jaw . Moreover, as already happened for the 4-year HESE data, for all these soft-hard combinations with a large branching ratio into the soft channel, the best fit for the DM lifetime is in tension with the gamma-ray limits from Refs. Cohen:2016uyg ; Blanco:2018esa , while compatible with the limits from Ref. Murase:2015gea (see the comment at the end of section 3.1.2).

As an illustration, in Figure 7 we show the best-fit event spectrum (among the combinations in Table 2), where we clearly see the interplay between the spectra of the two decay channels. In comparison with results obtained with the 4-year HESE data, the increase in events at EM-equivalent deposited energies of $\sim 100$ TeV makes it more difficult for the DM-only scenario to explain the PeV events simultaneously, making the overall fit slightly worse.

4 DM annihilations

The majority of the studies on the possible DM contribution to the IceCube high-energy neutrino flux have focused on the decay scenario, mainly because the required annihilation cross section, $\langle\sigma v\rangle$, for heavy DM ($m_{\text{DM}}\gtrsim 100$ TeV) Feldstein:2013kka ; Zavala:2014dla ; Esmaili:2014rma ; Chianese:2016opp ; Chianese:2016kpu violates the unitary bound $\langle\sigma v\rangle\leq 4\pi/(m_{\text{DM}}^{2}\,v)$ for thermal production Griest:1989wd ; Hui:2001wy . Yet, even if not produced thermally, the IceCube data could, in principle, test values below the local unitarity bound Esmaili:2014rma ; Chianese:2016kpu .999Note that the unitarity bound is less stringent locally, as the relative velocity of DM particles is smaller than in the early Universe. Moreover, the effect of DM annihilation in substructures, which would boost the signal with respect to the smooth contribution, along with a potential dependence of the relative velocity on negative powers, as in Sommerfeld-enhanced models Hisano:2004ds ; Profumo:2005xd ; Lattanzi:2008qa ; Feng:2010zp , could give rise to DM fluxes that can account for the observed number of high-energy neutrino-induced events Zavala:2014dla . In this section, however, we only consider a constant DM annihilation cross section and study the values that would give a potential contribution to the observed IceCube neutrino flux, in analogy to the case of decay presented above.

While in the case of the decaying DM scenario the resulting neutrino flux is linearly proportional to the (integral of the) DM number density, the neutrino flux from DM annihilations is proportional to the square of the DM density. Since the number density decreases with increasing $m_{\text{DM}}$, obtaining a neutrino flux comparable to that required for the IceCube event rate, requires a rather large annihilation cross section. As in the decaying DM case, the neutrino flux from DM annihilations at Earth receives both galactic and extra-galactic contributions. Nevertheless, there are important differences between these two cases.

The (smooth) galactic contribution of the neutrino flux from DM annihilation is given by

[TABLE]

where $\rho$ is the galactic DM density, $r(s,b,l)=\sqrt{s^{2}+R_{\odot}^{2}-2sR_{\odot}\cos b\cos l}$, is the distance to the Galactic center with $(b,l)$ being the Galactic coordinates and $R_{\odot}=8.5$ kpc. For the DM density we use an NFW profile Navarro:1995iw ; Navarro:1996gj , with the same parameters as in the decay case Bhattacharya:2017jaw . Given that the closer to the galactic center, the larger the DM density, this flux is more sharply peaked in that direction than in the case of DM decays.

In addition to the smooth galactic contribution, hierarchical structure formation in the standard cold DM scenario implies that DM halos must contain a large amount of smaller subhalos, which would further boost the annihilation flux. In this work, though, we do not account for this potential overall enhancement of the galactic DM annihilation signal. This boost factor depends on the direction of observation, as subhalos are not distributed homogeneously Springel:2008cc ; Hellwing:2015upa ; Rodriguez-Puebla:2016ofw and their properties also depend on the position within the main halo Springel:2008cc ; Diemand:2008in ; Pieri:2009je ; Moline:2016pbm . Nevertheless, we do not expect this to have an important impact on our results, as this enhancement is partly degenerated with the annihilation cross section and the ratio of the extragalactic to galactic contributions, which we keep as free parameters (see below).

The isotropic extragalactic contribution from DM annihilations is given by

[TABLE]

where $H(z)=H_{0}\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}$, and the cosmological parameters $H_{0}$, $\rho_{c}$, $\Omega_{m}$ and $\Omega_{\Lambda}$ are the same used for the decay scenario Bhattacharya:2017jaw . The dimensionless quantity $\zeta(z)$ accounts for the non-homogeneous DM clustering. It represents the variance of the DM density fluctuations at redshift $z$, so can be expressed in terms of the non-linear matter power spectrum Serpico:2011in ; Sefusatti:2014vha . Equivalently, within the halo model approach Cooray:2002dia , it is proportional to the cumulative sum of the enhancement in each individual halo at redshift $z$ Bergstrom:2001jj ; Ullio:2002pj ; Taylor:2002zd , and thus it depends on the abundance and internal properties of DM halos. Both approaches require extrapolations to small scales and high redshifts of the findings of N-body numerical simulations. This results in relatively large uncertainties in the computation of the clumping factor $\zeta(z)$, which can amount up to two or three orders of magnitude Sefusatti:2014vha ; Moline:2014xua . In this work we consider as our default $\zeta(z)$ the revised halofit plus stable clustering (RHF-SC) prescription of Ref. Sefusatti:2014vha , such that $(1+z)^{3}\,\zeta(z)\,H_{0}/H(z)\sim 10^{5}$. To take into account uncertainties, we introduce a constant normalization, $\xi$, where $\xi=1$ represents the default case and we vary it within $[10^{-3},10^{3}]$. Note, however, that the estimated uncertainty is likely smaller Sefusatti:2014vha ; Moline:2014xua . The total flux is given by

[TABLE]

Unlike what happens for the case of DM decays, for which there is no enhancement of the flux due to the linear dependence on the density, in the DM annihilating scenario, the enhancement could be significant and determines the importance of the extragalactic contribution: the larger $\xi\,(1+z)^{3}\,\zeta(z)\,H_{0}/H(z)$, the more isotropic the flux. Therefore, the value of $\xi$ serves as a measure for the anisotropy level in the total flux. As a reference, for $\xi=1$, the number of galactic events is about a factor of ten larger than the number of extragalactic ones, although the actual relative factor depends on the DM mass and annihilation channel.

As with DM decays, we consider two distinct scenarios here too:

1) neutrinos from DM annihilations into a single channel plus an isotropic astrophysical power-law flux, and

2) neutrinos solely from DM annihilations into two channels.

The flux in the former case is defined in terms of the set of free parameters $\bm{\theta}=\{\langle\sigma v\rangle,m_{\text{DM}},\xi,\phi_{\rm astro},\gamma\}$, while in the latter scenario, for each pair of channels, it is defined in terms of $\bm{\theta}=\{\langle\sigma v\rangle,m_{\text{DM}},\xi,{\rm BR}\}$, where ${\rm BR}$ is the branching ratio for annihilations into the first of the two channels.

4.1 Results: DM annihilations plus isotropic astrophysical power-law flux

First we consider the annihilations of DM particles into a single channel and their possible contribution to the observed neutrino flux by IceCube, in addition to an astrophysical power-law flux. The set of free parameters in the fit is $\bm{\theta}=\{\langle\sigma v\rangle,m_{\text{DM}},\xi,\phi_{\rm astro},\gamma\}$. The best-fit parameters corresponding to a total signal flux comprising neutrinos from a power-law astrophysical spectrum determined by $\phi_{\rm astro}$ and $\gamma$, and from DM annihilations, as a function of $\langle\sigma v\rangle$, $m_{\text{DM}}$ and $\xi$, for different two-body annihilation channels, are indicated in Table 3. The comparison of the best-fit likelihoods for all channels is shown in Figure 8. As in the case of DM decays, on top of each bar, we indicate the best fit obtained for the astrophysical index $\gamma$. In Figure 9 we show the event spectra for the best-fit channel, $\text{DM}\,\text{DM}\to W^{+}W^{-}$.

Qualitatively similar results are obtained for all channels. The preferred value for the DM mass lies in the range $\sim$ (160–500) TeV, and with an annihilation cross section of the order of $\sim(3\times 10^{-23}-5\times 10^{-21})~{}{\rm cm}^{3}~{}{\rm s}^{-1}$. Thus, the DM annihilation signal tends to explain the low-energy part of the event spectrum, whereas a relatively hard astrophysical flux (although softer than the through-going muon best fit Haack:2017dxi ) explains the highest energy events (Figure 9). In Table 3 we also show the number of DM events from the galactic and extragalactic contributions separately, whose relative importance is governed by the parameter $\xi$. For most channels, especially for the soft channels, the preferred value of $\xi$ reaches our lower boundary (i.e, $10^{-3}$), which results in a negligible amount of events produced from the DM extragalactic neutrino flux. In these cases, the astrophysical power-law flux is the main contribution to events from the Northern hemisphere. Note, however, that in the analysis we perform here the angular information is only taken into account at the hemisphere level, which could reduce the sensitivity to the $\xi$ parameter. In any case, the upper limit on the anisotropy parameter $\xi$ can be explained by data preferring the isotropic component to be decoupled from the DM contribution, even if this is not statistically significant yet. This is related to the combination of a preference for a mild anisotropy, as shown for the 3-year HESE data Esmaili:2014rma , and the typically hard spectrum of the DM signal. Therefore, below some value of $\xi$, the number of events from extragalactic DM annihilations would be very small and data is not sensitive anymore to this parameter.

4.1.1 Parameter correlations and preferred regions

The preference for small values of the parameter $\xi$ can be seen from the plots in the bottom row of Figure 10, which show different correlations for the $W^{+}W^{-}$ annihilation channel plus and astrophysical isotropic power-law flux. The vertical area in those plots visualizes the lack of sensitivity for small values of $\xi$.101010Note that for DM decays, the relative contribution from the galactic and extragalactic components only depends on the galactic DM profile, so there is no freedom analogous to the $\xi$ parameter. We can also clearly see the preference — although at less than 2$\sigma$ CL — for $m_{\text{DM}}\sim 150$ TeV and annihilations cross section orders of magnitude larger than the common thermal freeze-out value, and for a relatively hard astrophysical power-law flux, in acceptable agreement with the through-going muon spectrum Haack:2017dxi .

In summary, a common feature in all cases is the preference for a negligible or subdominant flux of extragalactic neutrinos from DM annihilations (i.e., small $\xi$) and with galactic DM annihilations only contributing to the low-energy part of the spectrum.

4.1.2 Limits on the DM annihilation cross section

As done for the case of DM decays, we can also use the 6-year HESE data to evaluate the maximum contribution to the event spectrum that could come from DM annihilations. Thus, we also compute the limits on the annihilation cross section as a function of the DM mass, for the two-body annihilation channels indicated in Table 3. All these results are shown in Figure 11 and can be easily understood by comparison with Figure 5. On the other hand, the unitarity bound in the local halo Griest:1989wd ; Hui:2001wy is approximately given by

[TABLE]

where $v_{\rm local}$ is the typical local relative velocity of DM particles.

For all the channels we consider, this bound (not shown in the figure) is below the best fit point. In particular, soft channels (upper panels in Figure 11) are in very strong tension with unitarity constraints, with differences of about two orders of magnitude. And even in the case of the hardest channels, data prefers values of the annihilation cross section a factor of a few above the unitarity limit. Note, however, that this bound should be taken as an order of magnitude estimate. There are different effects that could modify it, as the average over the galactic distribution or the redshift dependence when applied to the cosmological contribution.

4.2 Results: DM annihilations via multiple channels