Decaying dark matter at IceCube and its signature on High Energy gamma experiments

TL;DR

This paper investigates whether decaying dark matter could explain IceCube neutrino observations by analyzing multimessenger data from neutrinos and gamma rays, and assesses future gamma-ray experiments' potential to test this hypothesis.

Contribution

It provides a comprehensive multimessenger analysis of decaying dark matter models using IceCube and Fermi data, and explores future gamma-ray experiments' capabilities.

Findings

Decaying dark matter models are consistent with IceCube data within certain parameter ranges.

Fermi gamma-ray data constrains some dark matter decay channels.

Future gamma-ray experiments could further test the dark matter hypothesis.

Abstract

The origin of neutrino flux observed in IceCube is still mainly unknown. Typically two flux components are assumed, namely: atmospheric neutrinos and an unknown astrophysical term. In principle the latter could also contain a top-down contribution coming for example from decaying dark matter. In this case one should also expect prompt and secondary gammas as well. This leads to the possibility of a multimessenger analysis based on the simultaneous comparison of the Dark Matter hypothesis both with neutrino and high energy gamma rays data. In this paper, we analyze, for different decaying Dark Matter channels, the 7.5 years IceCube HESE data, and compare the results with previous exclusion limits coming from Fermi data. Finally, we test whether the Dark Matter hypothesis could be further scrutinised by using forthcoming high energy gamma rays experiments.

Figures

Tables

Table 1. Table 1 : Contributions to the gamma spectrum: Bremsstrahlung (B), Synchrotron (S), Inverse Compton (IC) on Cosmic Microwave Background (CMB), Inverse Compton on Star Light (SL) and Infrared (IR) photons. The components which we take into account are denoted by a star.

	Primary	Secondary
Galactic	Prompt*	B, S, IC by SL+IR, IC by CMB*
Extragalactic	Prompt with attenuation*	Extra-Galactic IC by CMB with attenuation*

Table 2. Table 2 : Relation between deposited energies and expected true energies of the incident neutrino for some benchmark energy values in the three topologies of Track (T), Shower (S) and Neutral Current (NC) events.

$E_{dep}\left(\mathrm{TeV}\right)$	$E_{\nu }^T\left(\mathrm{TeV}\right)$	$E_{\nu }^S\left(\mathrm{TeV}\right)$	$E_{\nu }^{NC}\left(\mathrm{TeV}\right)$
100	278	139	294
200	570	271	531
300	830	389	736
400	1040	515	912

Table 3. Table 3 : Summary of the neutrino analysis: the best fit parameters are given for each analyzed channel. The ϕ 0 subscript italic-ϕ 0 \phi_{0} is expressed in flux units of f.u.  ≡ TeV − 1 ​ cm − 2 ​ s − 1 ​ sr − 1 absent superscript TeV 1 superscript cm 2 superscript s 1 superscript sr 1 \equiv\rm TeV^{-1}cm^{-2}s^{-1}sr^{-1} .

Channel	${\varphi }_0^{\mathrm{best}}(\times {10}^{-15}f.u.)$	${\gamma }^{\mathrm{best}}$	${\tau }^{\mathrm{best}}(\times {10}^{28}s)$	$m_{\mathrm{DM}}^{\mathrm{best}}\left(\mathrm{TeV}\right)$
${\nu }_e{\nu }_e$	2.24	3.33	19.10	4017.35
${\nu }_{\mu }{\nu }_{\mu }$	2.24	3.33	19.10	4017.35
${\nu }_{\tau }{\nu }_{\tau }$	2.24	3.33	19.10	4017.35
$e^{+}{e}^{-}$	2.14	3.86	2.09	3846.63
${\mu }^{+}{\mu }^{-}$	0.66	2.64	1.91	569.17
${\tau }^{+}{\tau }^{-}$	0.74	2.69	1.59	570.00
$W^{+}{W}^{-}$	0.68	2.67	0.53	620.81
$ZZ$	0.72	2.69	0.63	621.00
$hh$	0.67	2.66	0.39	645.65
$b\overline{b}$	1.15	3.19	0.83	9168.11
$c\overline{c}$	0.78	2.78	0.40	3376.76
$t\overline{t}$	0.73	2.69	0.25	776.47
$q\overline{q}$	0.88	2.81	0.44	3233.26
$gg$	0.77	2.74	0.40	3526.63

Table 4. Table 4 : Summary of the neutrino analysis with a prior distribution for the through-going muons: the best fit parameters are given for each analyzed channel. The ϕ 0 subscript italic-ϕ 0 \phi_{0} is expressed in flux units of f.u.  ≡ TeV − 1 ​ cm − 2 ​ s − 1 ​ sr − 1 absent superscript TeV 1 superscript cm 2 superscript s 1 superscript sr 1 \equiv\rm TeV^{-1}cm^{-2}s^{-1}sr^{-1} .

Channel	${\varphi }_0^{\mathrm{best}}(\times {10}^{-15}f.u.)$	${\gamma }^{\mathrm{best}}$	${\tau }^{\mathrm{best}}(\times {10}^{28}s)$	$m_{\mathrm{DM}}^{\mathrm{best}}\left(\mathrm{TeV}\right)$
${\nu }_e{\nu }_e$	0.88	2.43	6.92	501.19
${\tau }^{+}{\tau }^{-}$	0.77	2.42	1.91	569.17
$W^{+}{W}^{-}$	0.77	2.42	0.63	575.44
$b\overline{b}$	0.86	2.43	0.14	776.25
$t\overline{t}$	0.83	2.43	0.30	676.08

Table 5. Table 5 : Number of standard deviations with which we expect CTA to detect decaying Dark Matter signals for all channels: the DM parameters are fixed at their IceCube best fit values reported in Table  3 .

Channel	$N_{\sigma }$ CTA sensitivity
$\nu \nu$	0.00066
$e^{+}{e}^{-}$	0.012
${\mu }^{+}{\mu }^{-}$	0.038
${\tau }^{+}{\tau }^{-}$	0.046
$W^{+}{W}^{-}$	0.23
$ZZ$	0.20
$hh$	0.35
$b\overline{b}$	0.11
$c\overline{c}$	0.28
$t\overline{t}$	0.56
$q\overline{q}$	0.25
$gg$	0.27

Equations

$$\chi$$

$${\varphi }_{\chi }={\varphi }_{\mathrm{G}}+{\varphi }_{\mathrm{EG}}$$

$$\frac{d^2{\varphi }_G}{dEd\Omega }=\frac{1}{4\pi m_{\mathrm{DM}}\tau }\frac{d{N}_{\nu +\bar{\nu }}}{dE}\int ds\rho (s,l,b)$$

$$\rho (r(s,l,b))={\rho }_0{\left[\frac{r}{r_c}{\left(1+\frac{r}{r_c}\right)}^2\right]}^{-1}$$

$$r(s,l,b)=\sqrt{s^2+R_s^2-2s{R}_s\cos b\cos l}$$

$$\frac{d{\varphi }_{\mathrm{EG}}}{dEd\Omega }=\frac{{\Omega }_{\chi }{\rho }_{cr}}{4\pi m_{\mathrm{DM}}\tau }\int \frac{dz}{H(z)}\frac{d{N}_{\nu +\bar{\nu }}(E(1+z))}{dE}$$

$$\frac{d{N}_i}{dE}(m_{\mathrm{DM}},E)=\frac{d{N}_i}{dE}\left(M_{\Lambda },\frac{E}{m_{\mathrm{DM}}}{M}_{\Lambda }\right)\frac{M_{\Lambda }}{m_{\mathrm{DM}}}$$

$$\frac{d^2{\varphi }_{\mathrm{G}}}{dEd\Omega }=\frac{1}{4\pi m_{\mathrm{DM}}\tau }\frac{d{N}_{\gamma }}{dE}\int ds\rho (s,l,b)$$

$$\frac{d^2{\varphi }_{\mathrm{EG}}}{dEd\Omega }=\frac{{\Omega }_{\chi }{\rho }_{cr}}{4\pi m_{\mathrm{DM}}\tau }\int \frac{dz}{H(z)}\mathrm{Att}(E,z)\frac{d{N}_{\gamma }(E(1+z))}{dE}$$

$$\nabla \cdot \left(D\nabla n(\mathbf{x},E)\right)-\frac{\partial }{\partial E}\left(b(E)n(\mathbf{x},E)\right)=q$$

$$n(\mathbf{x},E)=\frac{1}{b(E)}{\int }_E^{m_{\mathrm{DM}}/2}q(E^{\prime })d{E}^{\prime }$$

$$b_{\mathrm{IC}}(E)=3{\sigma }_T{\int }_0^{+\infty }ϵdϵ{\int }_{\frac{1}{4{\gamma }^2}}^{1}dqn(ϵ,\mathbf{x})\frac{(4{\gamma }^2-{\Gamma }_{ϵ})q-1}{(q+{\Gamma }_{ϵ}q{)}^3}$$

$$\left[2q\ln q+q+1-2q^2+\frac{({\Gamma }_{ϵ}q{)}^2(1-q)}{2(1+{\Gamma }_{ϵ}q)}\right]$$

$$b_{\mathrm{SYN}}(E)=\frac{4{\sigma }_T{E}^2}{3m_e^2}\frac{{\mathbf{B}}^2(\mathbf{x})}{2}$$

$$B(r,z)=B_0\exp \left\{-\frac{|r-R_{sun}|}{r_B}-\frac{|z|}{z_B}\right\}$$

$$P_{\mathrm{IC}}(E_{\gamma },E_e)=\frac{3{\sigma }_T{E}_{\gamma }}{4{\gamma }^2}{\int }_{1/4{\gamma }^2}^{1}dq\left[1-\frac{1}{4q{\gamma }^2(1-ϵ)}\right]\frac{n(E_{\gamma }^0)}{q}$$

$$\left[2q\ln q+q+1-2q^2+\frac{{ϵ}^2(1-q)}{2(1-ϵ)}\right]$$

$$\frac{d^2{\varphi }_{\mathrm{IC}}^{\mathrm{G}}}{dEd\Omega }=\frac{1}{2\pi E_{\gamma }{m}_{\mathrm{DM}}\tau }{\int }_0^{+\infty }ds\rho (r(s,l,b)){\int }_{E_{\gamma }}^{m_{\mathrm{DM}}/2}d{E}_e\frac{P_{IC}(E_{\gamma },E_e)}{b(E_e)}{\int }_{E_e}^{m_{\mathrm{DM}}/2}\frac{d{N}_e}{d{E}^{\prime }}d{E}^{\prime }$$

$$\frac{d^2{\varphi }_{\mathrm{IC}}^{\mathrm{E}}}{dEd\Omega }=\frac{{\rho }_{cr}{\Omega }_{\chi }}{2\pi E_{\gamma }{m}_{\mathrm{DM}}\tau }{\int }_0^{+\infty }\frac{dz}{H(z)(1+z)}{\int }_{E_{\gamma }(1+z)}^{m_{\mathrm{DM}}/2}d{E}_e\frac{P_{IC}(E_{\gamma }(1+z),E_e,z)}{b(E_e)}$$

$${\int }_{E_e}^{m_{\mathrm{DM}}/2}\frac{d{N}_e}{d{E}^{\prime }}d{E}^{\prime }$$

$$\varphi ={\varphi }_{\mathrm{Astro}}+{\varphi }_{\chi }$$

$$\frac{d^2{\varphi }_{\mathrm{Astro}}}{dEd\Omega }={\varphi }_0{\left(\frac{E}{100\mathrm{TeV}}\right)}^{-\gamma }$$

$$p_S^{e,\tau }$$

$$p_{NC}^{e,\tau }$$

$$N_T$$

$$N_S$$

$$N_{NC}$$

$$\mathcal{L}=\prod _i({\mu }_i+b_i{)}^{n_i}\frac{e^{-({\mu }_i+b_i)}}{n_i!}$$

$$\Lambda =2\ln \frac{\mathcal{L}(\hat{\theta })}{\mathcal{L}({\theta }_0)}$$

$${\mathcal{L}}_{\mathrm{TG}}=\exp \left[-\frac{(\gamma -\bar{\gamma }{)}^2}{2{\sigma }_{\gamma }^2}-\frac{(\varphi -\bar{\varphi }{)}^2}{2{\sigma }_{\varphi }^2}\right]$$

$$n=\Delta t{\int }_{\Delta E}dE\int d{E}^{\prime }R(E,E^{\prime })A_{eff}(E^{\prime })\frac{d\Phi }{d{E}^{\prime }}(E^{\prime })$$

$$\bar{\Lambda }=\sum _i\left[(a_i+b_i)\ln \left(\frac{a_i+b_i}{b_i}\right)-a_i\right]$$

$$\sigma =\sqrt{\sum _i\left[(a_i+b_i){\ln }^2\left(\frac{a_i+b_i}{b_i}\right)\right]}$$

Full text

I Introduction

After the IceCube recent measurements of ultra-high energy (UHE) neutrinos at TeV-PeV Aartsen et al. (2013, 2014, 2017a) novel possibilities have opened up for astroparticle physics. We are specially interested for the present work in the High-Energy Starting Events (HESE) data set. While it is clear that these neutrinos cannot be explained in terms of the atmospheric neutrino background, their astrophysical origin is still not fully understood. The astrophysical processes which are expected to be the source of these highly energetic neutrinos occur through the decay of charged pions, in turn originated in hadronic Loeb and Waxman (2006); Murase et al. (2013); Tamborra et al. (2014); Bechtol et al. (2017) and photo-hadronic Winter (2013); Murase et al. (2016) interactions. Together with neutrinos, gamma rays are produced as well. This simultaneous production has recently been confirmed by the multi-messenger observation of the coincident gamma-ray and neutrino emission from the flaring blazar TXS 0506+056 Aartsen et al. (2018a, b); Ansoldi et al. (2018). It has been nevertheless pointed out that blazar flares can only contribute up to about 10% of the total observed neutrino flux Murase et al. (2018). Moreover, other analyses of spatial and temporal correlations with gamma-rays Adrian-Martinez et al. (2016); Aartsen et al. (2017b, 2019a, 2019b) and for angular clustering Murase and Waxman (2016); Aartsen et al. (2015a); Ando et al. (2017); Mertsch et al. (2017); Dekker and Ando (2019) have placed strong constraints on the contribution of several extragalactic astrophysical sources. Allowed astrophysical candidates are therefore hidden cosmic-ray accelerators for which the gamma-ray flux is highly suppressed Murase et al. (2016); Senno et al. (2016). The neutrino flux, which is therefore expected to be the superposition of many astrophysical sources, can be parameterized as a diffuse power-law spectrum $E_{\nu}^{-\gamma}$.

At the same time there are however tensions among different IceCube data samples. Under the assumption of a single power-law flux, the latest 10-year through-going (TG) muon events lead to a spectral index of $2.28^{+0.08}_{-0.09}$ as best-fit Stettner (2019). This is mildly compatible with the theoretical expectations from a neutrino flux originating by photohadronic production from charged particles accelerated through the Fermi mechanism Waxman and Bahcall (1999) and with the measurements of the blazar TXS 0506+056, which suggest a hard power-law flux with a spectral index in the range $2.0\div 2.3$ Aartsen et al. (2018b). However, the latest $7.5$ years High Energy Starting Events (HESE) have the best-fit $\gamma=2.89^{+0.20}_{-0.19}$ Schneider (2019).

A spectral index of about $2.3$ is in tension with this result. This might suggest the presence of two components that dominate the diffuse neutrino flux at different energies with different angular properties Aartsen et al. (2017a, 2015b). This last conclusion is favored by the observation that the TG data sample collects muon neutrinos from the Northern hemisphere only, with energies larger that 200 TeV. The HESE events cover instead the whole sky, so they are sensitive to the galactic centre of the Milky Way; further, they are collected in an energy range starting already at 20 TeV. Investigations in the framework of two different power law sources have led to the conclusion that the data are compatible with the sum of a hard isotropic extragalactic neutrino flux and an additional softer one with a potential galactic origin Chen et al. (2015); Palladino and Vissani (2016); Vincent et al. (2016); Palladino et al. (2016); Anchordoqui et al. (2017); Palladino and Winter (2018); Sui and Bhupal Dev (2018). This two-component hypothesis is at the same time supported by the first combined analysis of IceCube and ANTARES data Chianese et al. (2017a). It turns out that both IceCube and ANTARES telescopes have measured in the same energy range (about 40–200 TeV) a slight excess with respect to an astrophysical power-law flux deduced by TG data ($\gamma\leq 2.2$), after the background subtraction Chianese et al. (2017a, 2016, b, c).

An alternative source for this diffuse UHE neutrino flux is decaying Dark Matter (DM) Chianese et al. (2016, 2017b, 2017c); Anisimov and Di Bari (2009); Feldstein et al. (2013); Esmaili and Serpico (2013); Bai et al. (2016); Ema et al. (2014a); Esmaili et al. (2014); Bhattacharya et al. (2014); Higaki et al. (2014); Rott et al. (2015); Ema et al. (2014b); Murase et al. (2015); Dudas et al. (2015); Fong et al. (2015); El Aisati et al. (2015); Ko and Tang (2015); Dev et al. (2016a); Re Fiorentin et al. (2016); Di Bari et al. (2016); Anchordoqui et al. (2015); Dev et al. (2016b); Chianese and Merle (2017); Borah et al. (2017); Boucenna et al. (2015); Hiroshima et al. (2018); Bhattacharya et al. (2017, 2019); Kopp et al. (2015). While another available source might be identified with annihilating DM Chianese et al. (2016); Feldstein et al. (2013); Zavala (2014); El Aisati et al. (2017); Sui and Bhupal Dev (2018); Kachelriess et al. (2018); Chianese et al. (2018); Bhattacharya et al. (2019), the unitarity limit leads in general to small neutrinos fluxes which are therefore not detectable. Bounds are available in previous studies both for decaying Aartsen et al. (2018c); Abeysekara et al. (2018) and annihilating DM Albert et al. (2017); Gozzini et al. (2019). Analyses mainly devoted to the neutrino energy spectrum are able to distinguish among different DM decay channels. Moreover, different DM models are further constrained through the informations coming from gamma-ray observations. The current data (neutrinos and gamma-rays) show that hadronic final states are strongly disfavored or excluded, while leptonic ones are slightly disfavored or allowed Cohen et al. (2017); Blanco and Hooper (2019). The most favorable case is a heavy DM candidate coupled only to neutrinos.

The foregoing observations have shown that gamma-rays constraints, mainly coming from observations of the Fermi-LAT experiment, have been critical in excluding or disfavoring various DM decay channels. It is therefore all the more interesting to find out whether gamma-rays experiments in different energy ranges would be able to put further constraints on the existence of decaying DM fluxes. While the lower energy range (up to few hundreds GeV) is mainly explored by Fermi LAT Atwood et al. (2009), the higher energy range (from about $100$ TeV) is explored by air shower experiments Kalashev and Kuznetsov (2016) and has already been used in exploring constraints on annihilating or decaying Dark Matter through data by Pierre Auger Aab (2015), CASA-MIA Chantell et al. (1997) and KASCADE Kang et al. (2015); a future experiment along these lines will be LHAASO Cao (2010). We are much interested in the intermediate range. Here the primary sources of information come from Cherenkov telescopes experiments such as the bounds coming from MAGIC Acciari et al. (2018), HESS Rinchiuso (2019), Veritas Zitzer (2018), HAWC Abeysekara et al. (2018) and CARPET-3 Dzhappuev:2018bnl . HAWC in particular has been able to put bounds on decaying Dark Matter, which are not, however, competitive with the Fermi-LAT ones. Data coming from these experiments have mainly been used in point-source investigations of Dark Matter decays and annihilation. A crucial development in this class will be the future CTA Ong (2019), which will overcome the previously mentioned Cherenkov experiments by an order of magnitude in its sensitivity. A number of studies on Dark Matter detection at CTA are already available, both for annihilation Ibarra et al. (2015); Silverwood et al. (2015); Hryczuk et al. (2019) and for decay Pierre et al. (2014). It is further interesting to notice that an excess in the TeV region has also been observed by Fermi LAT Neronov:2018ibl ; Neronov:2019ncc ; Ackermann:2014usa . This excess has been interpreted in terms of a decaying Dark Matter signal as well.

The primary aim of this paper is to provide an updated analysis of the IceCube data at $7.5$ years in the two-component hypothesis of an astrophysical spectrum, parameterized as a power law with a spectral index assumed as a fit parameter in the range between $1.5$ and $5$, and a diffuse spectrum produced from the decaying Dark Matter. We will further explore the possibilities of constraining the decaying Dark Matter hypothesis both through diffuse and through point-like investigations in the region of interest to IceCube.

The analysis we have performed is described in the following sections. In Section II we discuss the computation of the expected number of neutrino events at IceCube coming from decaying Dark Matter. In Section III the same ideas are generalized to the expected gamma fluxes on Earth: this case is more complicated, due to the presence of photon absorption and secondary production. In Section IV we describe the statistical procedure which has been adopted to compare the theoretical expectations with the latest HESE IceCube data, and the results of such analysis. Finally, in Section V we perform an investigation of the detection possibilities for the gamma rays fluxes both from diffuse and point-like searches. In Section VI, we draw our conclusions.

II Neutrinos from decaying Dark Matter

The neutrino flux from decaying Dark Matter unavoidably depends upon the decay mechanism, and, specifically, upon the final products of the decay. A complete analysis of the possible models which can be constructed is beyond the scope of this paper Bertone et al. (2005). A practical way to proceed is to consider the decays of Dark Matter into a pair of Standard Model particle and antiparticle. In the following we will therefore consider separately the following decaying channels:

[TABLE]

where we have denoted by $\chi$ the Dark Matter candidate. We do not consider here the $\gamma\gamma$ decay channel, which is severely constrained by gamma-ray experiments. The $d$ and $s$ quark channels give essentially the same results as the $u$ channel. Neutrino fluxes from Dark Matter decay are generated both inside the Galaxy and from extragalactic Dark Matter sources:

[TABLE]

While for the former we can assume a Navarro-Frank-White (NFW) profile density, for the latter we have to assume an isotropic distribution in accordance with the Cosmological Principle. The relevance of the choice in the Galactic Dark Matter profile density is not as strong as in the analysis of annihilating Dark Matter, due to the dependence in the latter on the square of the Dark Matter density. More detailed analyses of the effects of the profile choice can be found in Cirelli et al. (2011); Pierre et al. (2014). The Galactic component of the neutrino flux can then be represented as:

[TABLE]

where $m_{\rm DM}$ is the Dark Matter mass, $\tau$ is the Dark Matter lifetime, $dN_{\nu+\bar{\nu}}/dE$ is the differential spectrum of neutrinos produced per Dark Matter decay and the integral is taken over the line-of-sight. The angular dependence is parameterized through the density dependence on the Galactic coordinates, dictated by the NFW profile:

[TABLE]

where:

[TABLE]

The numerical values for the parameters appearing in the equations are $\rho_{0}=0.33$ $\rm{GeVcm^{-3}}$, $r_{c}=20$ $\rm{kpc}$, $R_{s}=8.5$ $\rm{kpc}$. The extragalactic component is in the same way obtained from summing the neutrinos produced from different cosmological distances parameterized by the redshift. Here we have to take into account the physical difference between the energy at the moment of production and detection of the neutrino, as well as the difference between the decay rate at production and detection of neutrinos. The result for the flux is:

[TABLE]

where $H(z)=H_{0}\sqrt{\Omega_{\Lambda}+\Omega_{m}(1+z)^{3}}$ is the Hubble parameter, $\Omega_{i}=\rho_{i}/\rho_{cr}$ and $\rho_{cr}$ is the critical density of the Universe. The fluxes have to be computed for all three flavors separately, since, as we will see, they have to be treated in different ways in the analysis of the neutrino detection. We also take into account the neutrino oscillations. We are now going to describe in more detail the computation of the neutrino spectra per decay.

The spectrum of particles produced in a decay or a collision of Dark Matter particles is typically obtained using simulation programs such as HERWIG Bahr et al. (2008) or PYTHIA Sjöstrand et al. (2015). For Dark Matter masses lower than about $100$ $\rm{TeV}$ a package named PPPC is available Cirelli et al. (2011) which already contains the spectra interpolated for most of the decay channels. A critical factor for such high masses of the Dark Matter is the showering which immediately follows the initial decay. While the electromagnetic part of the shower is amenable to analytical treatment and generally included in the softwares for simulation, and the hadronic part is already included both in HERWIG and PYTHIA, the weak showering is harder to describe. The reason is that, while in the first two cases the Sudakov double logarithms appearing in the probabilities for gauge boson emission cancel with the virtual corrections for the process with no gauge boson emission, in the weak case the $W$ boson emission changes the final state, so that the two processes lead to physically distinguishable final states. The appearance of double logarithms renders the perturbation results invalid already at energies of about $100$ $\rm{TeV}$, and in fact the results of Cirelli et al. (2011), obtained through perturbation theory, are only valid up to this energy. The higher energy range requires a resummation to all orders of the process. This has been partly implemented 111With some limitations which are described in detail in Christiansen and Sjöstrand (2014). in PYTHIA. Different approach has been given for instance in Kachelriess et al. (2018); Berezinsky et al. (2002).

Although the resummation of logarithms is not a trivial procedure, a guess which can be made based on pure dimensional analysis is that these differential spectra only depends on the ratio $E/m_{\rm DM}$. In fact, the only other quantity coming from weak physics and entering the result could be the mass $M_{W}$ of $W$ boson, but at energies much higher than this mass we would expect the result to be independent of this parameter. Of course the guess is too naive, since logarithms depend singularly on the ratio $m_{\rm DM}/M_{W}$ and might give rise to anomalous behaviors. We have therefore tested this hypothesis by comparing the spectra obtained using PYTHIA at masses of orders up to the $\rm{PeV}$ and the same spectra obtained by applying the scaling hypothesis to the spectra in Cirelli et al. (2011) (see figure 1). The comparison shows that the scaling hypothesis is justified and we have therefore used it in our calculations. This amounts to using the form

[TABLE]

as energy spectrum for Dark Matter masses $m_{\rm DM}$ larger than a cutoff mass $M_{\Lambda}=200$ TeV, which is the upper limit of the interpolated spectra provided in Cirelli et al. (2011).

III Gamma rays from decaying dark matter

A further source of information from decaying Dark Matter might come from the photons produced during the decay and in the subsequent propagation of particles from their decay place. We here discuss this aspect of our analysis. The production of gamma rays from Dark Matter decays has already been studied in various works (see, e.g., Cohen et al. (2017), Esmaili and Serpico (2015) and, at lower energies, Buch et al. (2015)).

This production can be divided, as in the case of neutrinos, in a Galactic and an Extragalatic component. There are, however, two main differences between photons and neutrinos. On the one hand, while the latter have roughly no interaction whatsoever during their propagation to Earth, the former are able, at high energies, to scatter off photons and electrons, accounting for an effective attenuation which, while for the Galactic component starts giving sensible consequences at electron energies higher than about $100$ $\rm{TeV}$ Esmaili and Serpico (2015), is already important at lower energies for the Extragalactic component. On the other hand, while neutrinos are only produced during Dark Matter decays or in the electroweak and hadronic cascade which immediately follows, photons can also be produced by secondary electrons and positrons by the mechanisms of bremsstrahlung, synchrotron radiation (in the Galactic magnetic field) and Inverse Compton from various sources. A numerical estimate shows that at high energies the relevant process is Inverse Compton. Furthermore, the analysis in Esmaili and Serpico (2015) shows that at such high energies the most important target for Inverse Compton is the Cosmic Microwave Background (CMB). A summary of all the relevant components is given in Table 1.

We can therefore make the simplifying assumption that the four main sources of gamma rays are the so called prompt gamma rays, produced directly during the Dark Matter decay, both in Galaxy and outside Galaxy; and the secondary gamma rays, produced by Inverse Compton of the electrons and positrons produced during the Dark Matter decay.

III.1 Prompt

The prompt Galactic spectrum is obtained in the same way as the prompt neutrino spectrum described previously. Thus, we have:

[TABLE]

For the prompt Extragalactic spectrum the neutrino form needs to be corrected by an attenuation factor, coming from collisions $\gamma-\gamma$ and $\gamma-e$. This attenuation factor ${\rm Att}(E,z)$ as a function of the photon detected energy and the redshift has been obtained in Cirelli et al. (2011). The Extragalactic prompt spectrum can therefore be expressed as:

[TABLE]

III.2 Secondary

Both the Galactic and the Extragalactic secondary gamma ray spectra are produced by electrons and positrons, and therefore depend critically upon the distributions of the latter. Such distribution should in principle be obtained by solving the kinetic equations for the positrons in the Galaxy and outside the Galaxy. For these equations we can require the equilibrium solution, since the decaying Dark Matter can be regarded on cosmic time scales as a steady source of particles. As is known from the theory of cosmic rays, these equations would have the schematic form:

[TABLE]

Here $D$ is a diffusion constant originated from the scattering of positrons on the Alfvèn waves and magnetic field fluctuations, $b(E)$ is the energy loss of the electrons, which at high energies is mainly connected with synchrotron and Inverse Compton losses, $q$ is the source term, in our case described by the injected spectrum of electrons and positrons from decaying Dark Matter and $n(E)$ is the number density of electrons per unit energy interval. While the diffusion constant is known to weakly increase with energy through a power law with exponent smaller than $1$ (see, e.g., Buch et al. (2015)), the energy loss components increase much more rapidly. For example, the synchrotron energy loss grows as the square of energy, and the Inverse Compton losses increase in a slightly lower fashion. We can draw the conclusions that the diffusion in the energy space is much more important than the spatial diffusion, and therefore we can neglect the first term in (9). The solution of the resulting equation is then:

[TABLE]

The energy losses due to Inverse Compton on CMB are:

[TABLE]

where, using the notation of Esmaili and Serpico (2015), $\Gamma_{\epsilon}=4\epsilon\gamma/m_{e}$ and $\gamma$ is the Lorentz factor of the electron. The quantity $\sigma_{T}$ is the Thomson cross section. The synchrotron energy losses are Esmaili and Serpico (2015):

[TABLE]

where for the Galactic magnetic field we have adopted the parameterization

[TABLE]

with $B_{0}=4.78$ $\rm{\mu G}$, $R_{sun}=8.3$ $\rm{kpc}$, $r_{B}=10$ $\rm{kpc}$, $z_{B}=2$ $\rm{kpc}$, and $r$ and $z$ being the radial and vertical distances from the Galactic center, respectively. We have here neglected the turbulent halo magnetic field. This of course has its most important effects at the high energies, where the synchrotron energy losses become dominant over the Inverse Compton losses. However, as shown in Esmaili and Serpico (2015), these effects become relevant at photon energies of order $100$ $\rm{TeV}$.

The Inverse Compton spectrum $d^{2}\phi_{\rm GIC}/dEd\Omega$ will now be the convolution of this electron distribution $n(E)$ given in (10) with the Inverse Compton emissivity on CMB, namely:

[TABLE]

where $\epsilon=E_{\gamma}/E_{e}$, $E^{0}_{\gamma}=m_{e}^{2}E_{\gamma}[4qE_{e}(E_{e}-E_{\gamma})]^{-1}$ and $n(E)$ is the number density of CMB photons per unit energy, which is homogeneous. For the Galactic secondary spectrum, the result is:

[TABLE]

For the Extragalactic secondary component the only modification to take into account is that the integral over the line-of-sight is converted into an integral over the redshift, and the CMB temperature is redshifted, as well as the photon energies. The final result is therefore:

[TABLE]

where we have introduced the explicit dependence of $P_{\rm IC}(E_{\gamma},E_{e},z)$ through the CMB temperature dependence.

IV Fitting $7.5$ years HESE IceCube data

It is here assumed that neutrinos observed at IceCube arise from a double component (together with the atmospheric neutrino background):

[TABLE]

where $\phi_{\chi}$ has been discussed above while the astrophysical spectrum is parameterized as a power law spectrum isotropic and independent of flavors

[TABLE]

where $\gamma$ is the spectral index and $\phi_{0}$ is the normalization of the power law.

Having delineated the neutrino flux arriving at Earth, we now discuss the conversion of this flux into the number of events detected by IceCube. The complete procedure would require the implementation of the stochastic relation between the true energy of each neutrino and the electromagnetic energy deposited in the detector, as described in detail in Palomares-Ruiz et al. (2015). We have instead followed the procedure of connecting these two energies through a deterministic function which gives the expected true energy given the deposited one. This function has been obtained in D’Amico (2018) for the three different cases of charged current shower events, track events and neutral current shower events. We report in Table 2 the numerical values of the expected value for the true neutrino energy corresponding to some benchmark values of the deposited neutrino energy, as given in D’Amico (2018):

The procedure of prediction of the number of events per unit time and energy interval is therefore as follows. For a certain detected energy interval we have obtained the corresponding true energy interval for all three topologies of events mentioned above. For each such topology we have integrated the expected neutrino flux for the three flavors by the effective area $A_{eff}(E,\Omega)$ (function of the true energy and the arrival direction of neutrino) provided by the IceCube Collaboration Aartsen et al. (2013) over the mapped energy interval. In the case of the Galactic flux, which is anisotropic due to the inhomogeneous distribution in the Galaxy plane, we have integrated over the solid angle using the angular dependent effective areas provided by the Collaboration. The weight assigned to each topology for different flavors is obtained by the procedure delineated in Palladino et al. (2015). The probabilities for each flavor $\alpha=e,\mu,\tau$ to produce one of the three topologies, namely shower (S), track (T) and neutral current (NC) showers are given by

[TABLE]

where $\sigma_{CC}$ and $\sigma_{NC}$ are the charged and neutral current cross section Gandhi et al. (1998) and $M_{NC}$, $M_{CC}^{e,\mu,\tau}$ are the effective detector masses Aartsen et al. (2013). So given a generic set of three neutrino fluxes $\phi_{\alpha}$ the number of track and shower events in a certain energy interval is given by

[TABLE]

For the decaying dark matter case, the three neutrino fluxes $\phi_{\alpha}$ are function the specific decaying channel considered.

The fitting procedure of a generic neutrino input spectrum is then to use the methods discussed above to predict the expected number of events in each bin. We can then assume a multi Poisson distribution for each of the energy bin at energies above $60$ TeV, in agreement with the methods used by the IceCube Collaboration Ishihara (2016), to adapt the method of maximum likelihood. Let us denote the total number of expected events from astrophysical and Dark Matter sources in each energy bin by $\mu_{i}$, obtained according to the procedure described above, and by $b_{i}$ the expected number of background events according to the IceCube Collaboration. Furthermore, let $n_{i}$ be the number of measured events in each energy bin. Then the likelihood function adopted is:

[TABLE]

In this way we have a complete procedure for obtaining the number of total events expected from Dark Matter decay in an energy interval. In order to compare with the experimental data, of course, it is necessary to have an estimate of the background. This has been done by using the same binning as the one given by the IceCube Collaboration, which provides in graphical form for each bin the expected number of background events Schneider (2019).

The computation of the likelihood for a given channel and given values of the parameters is most easily accomplished by observing that $\tau$ and $\phi_{0}$ just serve as normalization respectively of the Dark Matter and astrophysical fluxes, whose non trivial dependence is therefore upon $m_{\rm DM}$ and $\gamma$. We have therefore obtained, for a reference normalization flux $\tilde{\phi}_{0}$, the expected number of astrophysical events for each energy bin as a function of the spectral index taken in a range between $1.5$ and $5$, and, in the same way, we have obtained, for a reference Dark Matter decay lifetime $\tilde{\tau}$, the expected number of Dark Matter events for each energy bin. The fluxes for a generic normalization flux $\phi_{0}$ and a generic lifetime $\tau$ are therefore obtained from these by multiplying the factors $\phi_{0}/\tilde{\phi}_{0}$ and $\tilde{\tau}/\tau$ respectively.

IV.1 One-component fit: pure astrophysical power law

A first analysis which can be done with the tools developed in the previous sections is the analysis of the IceCube data in the hypothesis of the pure astrophysical power law spectrum, which is the cumulative spectrum from the astrophysical sources. For this case with no Dark Matter, we find as best fit values $\gamma=3.02\pm 0.13$ and $\phi_{0}=\left(1.9\pm 0.3\right)\times 10^{-15}$ $\rm{TeV^{-1}}$ $\rm{cm^{-2}}$ $\rm{s^{-1}}$ $\rm{sr^{-1}}$. The uncertainties provided here are the $68\%$ confidence levels.

IV.2 Two-component fit: astrophysical power law and decaying dark matter

We now describe the results obtained for the cases of two components. The statistical analysis of this case can be performed with the aid of the maximum likelihood method, with the likelihood specified above. The determination of the confidence regions can be performed once the distribution of the likelihood as a test statistic is known. A theorem which is often adopted in this context is the Wilks theorem Wilks (1938), which states that the distribution of the log-likelihood ratio

[TABLE]

between the likelihood evaluated at its maximum ($\hat{\theta}$) and the maximum of the likelihood in a region which is being tested ($\theta_{0}$) is a chi-squared distribution with a number of degrees of freedom equal to the difference in the number of the free parameters between the null hypothesis and the whole parameter space. The hypothesis under which this theorem is proved is that the likelihood possesses a single maximum and behaves like a Gaussian distribution around this point. As will be seen below, this hypothesis is not realized in the distributions we find. Nevertheless, it can be proven that this hypothesis can be considerably relaxed. In fact, it turns out that the likelihood distribution is a chi-squared one for each point of the parameter space which sufficiently well reproduce, in its mean value, the experimental data. This implies that each set of parameters which produces an expected flux comparable to the data within the statistical uncertainties possesses a chi-squared distribution for the likelihood ratio, and this result does not depend on the geometry of the parameter space. Therefore, we can safely adopt the results of the theorem even if we are beyond the range of application of the theorem as originally stated. Furthermore, the Neyman-Pearson lemma grants us that, at a fixed level of significativity, this is the most powerful test statistic which we can adopt.

For the various channels analyzed, we represent the results in the form of contour plots of the likelihood ratio $\Lambda$ defined above in figures 2, 3 and 4 for the various channels, together with the $68\%$, $95\%$ and $99.7\%$ confidence levels. In these plots, the astrophysical parameters are treated as nuisance parameters, and the likelihood is correspondently maximized with respect to them: in particular, $\gamma$ is chosen in a range between $1.5$ and $5$. As anticipated above, these plots highlight a non trivial structure in the likelihood function, which posses different competing maxima. We expect a strong sensitivity in the choice of the dominant maximum on the details of our analysis, such as the adopted functions for the relation between deposited and real energy shown in Table (2).

Qualitatively, the data can be reproduced either by a softer spectrum, with a spectral index nearer to $2$, and a Dark Matter with a mass of order $\sim 100$ $\rm{TeV}$ which can fit the excess in this region, or by a harder spectrum with a DM mass of order $\sim 1$ $\rm{PeV}$, which can fit the high energy data. A major factor in determining which regions of the $m_{\rm DM}-\tau$ plane will be favored by the likelihood ratio is the choice of the interval over which $\gamma$ varies: in fact, according to our qualitative understanding of the solution, restricting $\gamma$ to be near the expected spectral index, which is around $2$, should favor the region of masses of order $\sim 100$ $\rm{TeV}$.

In figures 2, 3 and 4 we also represent, for the channels for which it is available, the exclusion curves obtained in Cohen et al. (2017) from the Fermi-LAT data with a red line. We also give in Table 3 the best fit values for the astrophysical and Dark Matter parameters. We underline that for the various channels the $2\sigma$ contours are always open, which means that the pure astrophysical spectrum cannot be rejected at more than the $2\sigma$ levels; further, both the hadronic and the neutrinophilic channels are not analyzable with more than the $1\sigma$ confidence level.

A point of special interest is the comparison between the best fit astrophysical parameters and the astrophysical parameters found by the analysis of the through-going muons. As was mentioned in the introduction, the latter predict a spectral index $\gamma$ of $2.28$. The comparison between the spectral indices can be directly performed by virtue of our expectation that the Dark Matter gives a relatively small contribution to the through-going muons, because of the higher energy range and the more limited angular range which are explored. As we mentioned before, the likelihood maxima with $m_{DM}\sim$ PeV have a $\gamma$ larger than $3$, while those with $m_{DM}\sim 100$ TeV have a $\gamma$ smaller than $3$. Therefore the former case leads to a worse agreement with the through-going muons, while in the latter we have an improvement. However, being the two maxima almost degenerate, we cannot unambiguously say that the Dark Matter deteriorate the agreement. The next section will be dedicated to a partial answer to this question.

IV.3 Astrophysical spectrum properties

At the end of the previous section we mentioned the importance of the spectral index in estimating the quality of our fit. In fact, as we stated in the introduction, Dark Matter is a possible candidate especially by virtue of the inconsistency of the spectral index found through the through-going muon neutrinos analysis and the one predicted by the HESE fit.

A rigorous investigation in this respect should be obtained by performing a combined fit on the through-going muon neutrinos data and the HESE data. This is however beyond the scope of the present paper. It is still of interest to take into account the information provided by the through-going muon neutrinos data approximately. This can be done by incorporating, in the likelihood (21), a Gaussian distribution for the astrophysical parameters centered around the values predicted by the IceCube fit with the corresponding uncertainties. In this way, we are introducing a Gaussian approximation to the through-going muon neutrinos likelihood which would be found by fitting them with a pure astrophysical spectrum. The likelihood (21) is therefore modified through the multiplication by a factor:

[TABLE]

where the parameters are extracted by the 10-year through-going muon neutrinos fit to be $\bar{\gamma}=2.28$, $\sigma_{\gamma}=0.08$, $\bar{\phi}=1.44\times 10^{-15}\rm\,TeV^{-1}cm^{-2}s^{-1}sr^{-1}$ and $\sigma_{\phi}=0.25\times 10^{-15}\rm\,TeV^{-1}cm^{-2}s^{-1}sr^{-1}$.

The incorrectness in our method steals from the fact that we are modeling the HESE data with a double component spectrum, while the through-going muon neutrinos are modeled by a pure astrophysical spectrum. However, due to the higher energy range and the different angular range, which does not involve the Galaxy center where the Dark Matter neutrino production is stronger, we expect the latter to be subdominant for the through-going muon neutrinos. This means that we can get at least a qualitative understanding of the changes which would be introduced in our results by taking into account the through-going muon neutrinos data. The contour plots for some significant channels, which best illustrate the qualitative changes, are represented in Figure 5 and Figure 6.

These plots evidences that the introduction of the through-going muons prior information strongly disfavors the PeV region for the DM masses, so that even those channels which previously had their best fit in that region are moved to the $\sim 100$ TeV region. The most dramatic consequence of this is that the hadronic channels are even more disfavored by the gamma ray constraints. The reason for this is that, as was mentioned above, the PeV region requires spectral indices even larger than $3$ to obtain a consistent fit of the data. An efficient comparison can be made by looking at the best fit parameters obtained in this new analysis: for the decay channels which have been analyzed these results are given in Table 4.

Another important consequence is that for all channels the $1\sigma$ contours are now closed. This can also be understood on the grounds that the pure astrophysical hypothesis requires, as we know, a spectral index near $3$, which is therefore relatively far from the through-going muon neutrinos value. Therefore, the region of large $\tau$ has a lower likelihood than before and is more strongly disfavored. Even with this information, however, the pure astrophysical spectrum can never be rejected at the $2\sigma$ level.

The question of the competing maxima, one for a Dark Matter mass of the order of PeV and one at the order of hundreds of TeV, can be made clearer by comparing the relative contribution of the astrophysical and Dark Matter spectrum in either of the two cases above mentioned, both with the introduction of the prior information and without it. We have done so in Figure 7 for two representative cases: in the case $\chi\rightarrow\nu_{e}\nu_{e}$ the best fit value of the Dark Matter mass is of order $\sim 1$ PeV, while in the case $\chi\rightarrow WW$ it is of order $\sim 100$ TeV. For both decay channels we represent the histogram for the best fit values in the case with prior information and without it. It can now be seen that in the two channels, for the case without prior information, the Dark Matter excess is relevant in two different regions: in the former only the higher energies are affected, while in the latter the lower energies of order $100$ TeV are influenced. The introduction of the prior distribution moves the Dark Matter contribution in both cases at lower energies, while favoring a slightly harder astrophysical spectrum, as is clear from the presence of more events at higher energies.

V High energy gamma from decaying dark matter: background comparison

As has been detailed in previous sections, a possible way of dealing with the investigation of Dark Matter decay is the analysis of the gamma ray fluxes. There are in fact a large number of experiments detecting high energy gamma rays which might be able to put further constraints on the models which we are analyzing. Some of these experiments generally work through point-like sources investigations: these include MAGIC Acciari et al. (2018), HESS Rinchiuso (2019), Veritas Zitzer (2018), HAWC Abeysekara et al. (2018), CARPET-3 Dzhappuev:2018bnl and the future CTA Pierre et al. (2014). A second class of experiments is instead able to detect the diffuse flux: these include Fermi-LAT Atwood et al. (2009), Pierre Auger Aab (2015), CASA-MIA Chantell et al. (1997) and KASCADE Kang et al. (2015). In the following we will qualitatively investigate their potentialities in detecting the Dark Matter produced photon fluxes for both kinds of experiments. In order for this strategy to be useful, a clear necessity is that the photon fluxes produced by the decay should be distinguishable from the background. We will therefore test whether this requirement is met in the cases of interest for IceCube.

Since most of the high energy experiments mentioned above typically investigate point-like sources, we will make the comparison for the gamma ray fluxes coming from a region of $2^{\circ}\times 2^{\circ}$ around the Galactic Center. We follow the analysis performed in Ibarra et al. (2015), studying a region with an aperture which is comparable to the field of view of the CTA experiment. The estimation of the background is performed similarly to Ibarra et al. (2015), taking into account different background components. The detailed procedure is left for the appendix.

In figure 8 we compare the background spectrum with the photon spectra from decaying Dark Matter for all the analyzed channels: the mass and lifetime of the Dark Matter is fixed to the best fit value found in the previous section for each different channel. The DM spectra are lower than the background by at least two orders of magnitude: we therefore reach the conclusions that CTA would not be sensitive to such signals and could therefore not provide further constraints on the region of the parameter space which is of special interest to IceCube. Accordingly, the same can be said for the other high energy experiments such as HESS, VERITAS and MAGIC having a lower sensitivity.

A more detailed analysis can be performed along the guidelines of Ibarra et al. (2015) to determine the significance with which CTA would be able to discriminate the presence of the predicted fluxes. In order to do so, we need to predict, for a given incident flux per unit energy, the number of events expected at CTA. This requires a knowledge of the CTA effective area $A_{eff}(E)$ and reconstruction function $R(E,E^{\prime})$. The latter is assumed to be Gaussian, with a $68\%$ containment interval provided by the CTA Collaboration, together with the effective area Ong (2019). In terms of these elements the expected number of events in an energy bin $\Delta E$ is:

[TABLE]

where $\Delta t$ is the $50$ hours exposition time.

We divide the energy interval from $25$ GeV to $100$ TeV into $200$ bins per decade and predict the number of background events $b_{i}$ and of Dark Matter decaying events $a_{i}$ for each energy bin. This analysis has been separately done for each decay channel, with the DM parameters fixed to their IceCube best fit value.

To foresee the sensitivity of CTA we need to assess the distribution of the likelihood ratio between the case of pure background and the case of background plus signal. This distribution has to be computed in the case of events distributed according to the background plus signal case222If the events were distributed according to the background only case, we would assess the exclusion p-value.. The likelihood is assumed to take the multi-Poissonian form (21). If we call $\Lambda$ the likelihood ratio, its distribution in the vicinity of the maximum can be well approximated by a Gaussian with a mean value of:

[TABLE]

and a standard deviation of:

[TABLE]

To assess the number of standard deviations with which it would be possible to exclude the pure background hypothesis, we notice that in such hypothesis the expected value of the likelihood ratio would be:

[TABLE]

It follows that the number of standard deviations for rejection of the pure background hypothesis would be the ratio:

[TABLE]

Performing the preceding analysis for each channel, we find the results in Table 5, where for each channel we provide the predicted number of standard deviations with which, in the case of signal plus background, it would be possible to reject the hypothesis of pure background. Confirming the intuitive conclusions coming from the comparison with the background, we see that for none of the analyzed channels it would be possible to discriminate the signal. This result is in agreement with Pierre et al. (2014). We also notice that in (28) we can expand in the small ratio $\frac{a_{i}}{b_{i}}$ to obtain the form:

[TABLE]

which shows that the number of standard deviations scales linearly with the signal amplitude for small signals. This means that, for example, the $t\bar{t}$ channel would reach the $5\sigma$ detection limit for a signal of order $\sim 10$ times larger, which corresponds to a lifetime smaller by a factor of $10$. By comparison with figure 4 we see that this value is lower than the Fermi-LAT exclusion limits, which confirms the expectation that analyses of point-like sources are generally not competitive in giving constraints on the Dark Matter fluxes.

A further possibility of investigation of the Dark Matter fluxes, as explained above, is the diffuse flux. We provide for reference the diffuse fluxes expected for all channels with a DM mass and lifetime fixed at their best fit values, together with the data collected by KASCADE Chantell et al. (1997) and CASA-MIA Schatz et al. (2003). In Figure LABEL:fig4 we show the diffuse fluxes which, as foreseen, are more constrained by the experimental data than the corresponding point-like investigations.

VI Conclusions

The multimessenger analysis has recently become a guiding instrument in the investigation of cosmological Dark Matter. Following these guidelines, in this work we have focused on the neutrino and gamma ray production from decaying Dark Matter. The unavoidable model dependence of such study has been reduced to the decay into particle-antiparticle pair from the Standard Model for various decay channels. In particular, we have analyzed the latest IceCube HESE data under the assumption of a dark component superimposed to an astrophysical cumulative spectrum. This constitutes an update of previous analyses, and therefore allows the determination of more stringent bounds in the decaying Dark Matter parameter space. The results for the best fit parameters show that all of the channels favor either a Dark Matter mass around $100$ TeV or a Dark Matter mass in the PeV region: in the former case the spectral index is expected to become nearer to the values of $2.3$, while in the latter it is expected to increase to higher values. The statistical test evidences that for none of the channels the pure astrophysical spectrum can be rejected with more than the $2\sigma$ level, and in some of the channels it cannot be rejected even at the $1\sigma$ level. This implies that the Dark Matter introduction does not substantially improve the data explanation.

Moreover, we have used, for each decay channel, the best fit scenario to the IceCube data to predict the expected gamma ray spectra. The purpose of this examination was to provide an estimation of the gamma rays experiment possibilities in further constraining the model. We have done so both for point-like and diffuse searches. In the former case, we have predicted the expected spectra from Dark Matter for all channels in the region around the Galactic Centre, where these are expected to be higher, with the modeled background from the same region. This analysis has been extended, for the specific case of the Cherenkov Telescope Array, to a statistical investigation of the p-value with which the future experiment should be able to assess the presence of a Dark Matter originated component. The results show that it would not be possible through a point-like analysis to significantly constrain the Dark Matter parameter space beyond the already existing constraints, in that the IceCube best fits are not testable even at the $1\sigma$ level. In the case of the diffuse searches, we have compared the predicted diffuse gamma ray spectra for all channels with the data coming from KASCADE, CASA-MIA and Fermi LAT. While the former two are able to constrain only Dark Matter with masses above the $10$ PeV, the latter is more informative in determining the admissible region of the parameter space, and has in fact already been adopted in the previous analysis of the IceCube data.

Acknowledgments

The authors thank Kohta Murase for his useful comments during the 36th International Cosmic-Ray conference in Madison. This work was partially supported by the research grant number 2017W4HA7S "NAT-NET: Neutrino and Astroparticle Theory Network” under the program PRIN 2017 funded by the Italian Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR). D.F.G.F., G.M., S.M. and O.P. acknowledge partial support by the research project TAsP (Theoretical Astroparticle Physics) funded by the Instituto Nazionale di Fisica Nucleare (INFN).

Appendix A Gamma rays background estimation for the Galactic Center

There are three main components to the gamma rays background at the Galactic Center: one coming from the showers induced by cosmic protons, one coming from the electron and positron induced showers and a gamma ray Galactic background. The component of background due to the showers produced by the cosmic protons is parameterized as Hoerandel (2003):

[TABLE]

Due to the different yield of the proton showers, this spectrum is actually calculated for an energy $E^{\prime}_{p}=\rho E_{p}$, where $\rho\simeq 3$. Another background component is the spectrum of the cascades induced by the cosmic electrons. For energies below $\sim 100\rm{GeV}$ this spectrum is measured by AMS-02 Aguilar et al. (2014) as:

[TABLE]

while above $\sim 30\rm{GeV}$ is measured by HESS as Aharonian et al. (2006a):

[TABLE]

To ensure a smooth transition between these two spectra, we take

[TABLE]

Finally, the gamma ray background is made up of mainly two components. On the one hand there is the Galactic ridge emission, which is here parameterized as Aharonian et al. (2006a):

[TABLE]

On the other hand, the second component can be identified as a point-like source. Its spectrum is in the range $5-100$ $\rm{GeV}$ Chernyakova et al. (2011):

[TABLE]

and in the range $160$ $\rm{GeV}-30$ $\rm{TeV}$ Aharonian et al. (2006b):

[TABLE]

The two spectra are interpolated by the form

[TABLE]

In the overall background flux the different efficiency in the detection of cosmic protons is taken into account by multiplying the corresponding flux by a factor $\epsilon_{p}=0.001+0.002\left(\frac{E_{p}}{20\rm{TeV}}\right)^{1.4}$.