Heavy decaying dark matter and large-scale anisotropy of high-energy cosmic rays

TL;DR

This paper investigates how large-scale anisotropy in high-energy cosmic rays can be used to detect or constrain heavy decaying dark matter particles, comparing anisotropy constraints with gamma-ray limits.

Contribution

It provides new constraints on dark matter particle lifetime using cosmic ray anisotropy data and discusses the measurement precision needed for future detection.

Findings

Anisotropy constraints are weaker than gamma-ray limits.

Estimated precision needed for anisotropy measurements.

Discussed prospects for dark matter detection with current facilities.

Abstract

We examine the role of the large--scale anisotropy of the high--energy cosmic ray distribution in a search for the heavy decaying dark matter (DM) signal. Using recent anisotropy measurements from the extensive air shower (EAS) observatories we constrain the lifetime of the DM particles with masses $10^{7}~\leq~M_X~\leq~10^{16}$ GeV. These constraints appear to be weaker than that obtained with the high energy gamma-ray limits. We also estimate the desired precision level for the anisotropy measurements to discern the decaying DM signal marginally allowed by the gamma-ray limits and discuss the prospects of the DM search with the modern EAS facilities.

Equations

$$J_{\exp }(E)=J_1\{\begin{array}{l}{E}_{\text{-1}}^{{\gamma }_{-1}-{\gamma }_0}{E}_0^{{\gamma }_0-{\gamma }_1}{E}^{-{\gamma }_{-1}}E<E_{-1}\\ E_0^{{\gamma }_0-{\gamma }_1}{E}^{-{\gamma }_0}{E}_{\text{-1}}\le E<E_0\\ E^{-{\gamma }_1}{E}_{\text{0}}\le E<E_{\text{1}}\\ E_{\text{1}}^{{\gamma }_2-{\gamma }_1}{E}^{-{\gamma }_2}{E}_{\text{1}}\le E<E_2\\ E_{\text{1}}^{{\gamma }_2-{\gamma }_1}{E}_2^{{\gamma }_3-{\gamma }_2}{E}^{-{\gamma }_3}E>E_2\end{array},$$

$$J_{\mathrm{DM}}(\delta ,\alpha ,E)=J_p^{\mathrm{G}}(\delta ,\alpha ,E)+J_{\gamma }^{\mathrm{G}}(\delta ,\alpha ,E)+J_p^{\mathrm{EG}}(E)+J_{\gamma }^{\mathrm{EG}}(E)$$

$$J_{\mathrm{tot}}(\delta ,\alpha ,E)=J_{\exp }(E)+J_p^{\mathrm{G}}(\delta ,\alpha ,E)+J_{\gamma }^{\mathrm{G}}(\delta ,\alpha ,E)$$

$$R_g=\frac{E}{qe{B}_{\perp }}\simeq 1.1\times \frac{1}{q}\left(\frac{E}{10^{18}\mathrm{eV}}\right){\left(\frac{B_{\perp }}{\mu \mathrm{G}}\right)}^{-1}\mathrm{kpc},$$

$$J(\alpha ,E)=A_0(E)+\sum _n\left[A_n(E)\sin (n\alpha )+B_n(E)\cos (n\alpha )\right],$$

$$A_0(E)=\frac{1}{2\pi }\underset{-\pi }{\overset{\pi }{\int }}J(\alpha ,E)d\alpha ,$$

$$A_n(E)=\frac{1}{\pi }\underset{-\pi }{\overset{\pi }{\int }}J(\alpha ,E)\cos (n\alpha )d\alpha ,$$

$$B_n(E)=\frac{1}{\pi }\underset{-\pi }{\overset{\pi }{\int }}J(\alpha ,E)\sin (n\alpha )d\alpha .$$

$$r_1=\sqrt{a_1^2+b_1^2}$$

$$\omega (a,\delta ,{\theta }_{\text{max}})\sim (\cos a\cos \delta \sin {\alpha }_m+{\alpha }_m\sin a\sin \delta ),$$

$${\alpha }_m=\{\begin{array}{ll}0& ;\xi >1,\\ \pi & ;\xi <-1,\\ \arccos \xi & ;-1<\xi <1;\end{array}$$

$$\xi =\frac{(\cos {\theta }_{\text{max}}-\sin a\sin \delta )}{\cos a\cos \delta }.$$

$$a_1(E)=\frac{2\underset{\Omega }{\int }{J}_{\mathrm{DM}}(\delta ,\alpha ,E)\omega (a,\delta ,{\theta }_{\text{max}})\cos (\alpha )d\Omega }{J_{\exp }(E)\underset{\Omega }{\int }\omega (a,\delta ,{\theta }_{\text{max}})d\Omega +\underset{\Omega }{\int }{J}_{\mathrm{DM}}(\delta ,\alpha ,E)\omega (a,\delta ,{\theta }_{\text{max}})d\Omega }$$

$$b_1(E)=\frac{2\underset{\Omega }{\int }{J}_{\mathrm{DM}}(\delta ,\alpha ,E)\omega (a,\delta ,{\theta }_{\text{max}})\sin (\alpha )d\Omega }{J_{\exp }(E)\underset{\Omega }{\int }\omega (a,\delta ,{\theta }_{\text{max}})d\Omega +\underset{\Omega }{\int }{J}_{\mathrm{DM}}(\delta ,\alpha ,E)\omega (a,\delta ,{\theta }_{\text{max}})d\Omega }.$$

$$J(\delta ,\alpha ,E)=\sum _{l\ge 0}\sum _{m=-l}^l{a}_{lm}(E){\mathrm{Y}}_{lm}(\delta ,\alpha )$$

$$a_{lm}(E)=\underset{\Omega }{\int }J(\delta ,\alpha ,E){\mathrm{Y}}_{lm}(\delta ,\alpha )d\Omega ,$$

$$C_l=\frac{1}{2l+1}\sum _m|a_{lm}{|}^2.$$

Full text

Keywords: dark matter, cosmic-ray anisotropy.

1 Introduction

The hypothesis of dark matter (DM) consisting of heavy long–living particles has attracted significant attention in the context of inflationary cosmology [1, 2]. There are several scenarios of effective DM particles production on various stages of early Universe evolution [3, 4, 1, 2, 5, 6, 7, 8, 9, 10, 11]. Although, heavy DM was discussed in other contexts as well [12, 13, 14]. From the experimental point of view the most appropriate method to search for such DM particles is to look for the secondary high–energy cosmic–ray flux from the particle decay. Historically the first indication on super–heavy DM existence came from the observation of super–GZK cosmic rays in AGASA [15]. However later on the GZK cut–off existence was confirmed by the next generation cosmic ray experiments [16, 17]. Several DM decay based interpretations [18, 19, 20, 21, 22, 23] have been proposed for the detection of PeV neutrinos in IceCube [24, 25]. While most of these interpretations are disfavored by the recent studies in which the respective gamma-ray signal [26, 27] is analyzed, some DM models with suppressed photon production are still viable [28, 21, 29].

Technically the heavy DM candidate $X$ has two main parameters: mass $M_{X}$ and lifetime $\tau$. Absolutely stable $X$–particles are not so interesting from the experimental point of view since its annihilation cross–section is bounded by unitarity: $\sigma_{X}^{\rm ann.}\sim 1/M_{X}^{2}$, which makes its indirect detection impossible for the todays experiments [30]. There are several sources of constraints for the heavy DM parameters. The mass is subjected to cosmological constraints [31, 8, 32, 10, 33, 34], and the lifetime of the DM particles can be effectively constrained with the observed fluxes of various high–energy particles or with the upper limit on these fluxes. For example, in Ref. [35] the constraints have been put using the shape of charged cosmic–ray spectra. However, with the modern cosmic ray data this method is not as effective in constraining $\tau$ as gamma–ray [36, 27, 37, 38, 39] and neutrino [27, 26] data.

Another observable sensitive to the heavy DM decay is the cosmic–ray anisotropy. Apart from the DM searches, the anisotropy is a powerful tool for the elucidation of the cosmic–ray origin and propagation. In particular it is useful in study of the galactic magnetic field structure imprinted in cosmic rays [40, 41, 42, 43] or in search for the extended cosmic–ray sources such as large scale structure [44, 45, 46, 47]. While for TeV–PeV energies the existence of large–scale anisotropy has been confirmed by several experiments [48, 49], for higher energies the deviations from isotropy are either not observed or have low significance [50, 51, 52, 53, 54, 55]. In the present work we use the upper–limits on the cosmic ray flux anisotropy mentioned above to obtain the conservative constrains on the lifetime of the heavy DM with masses $10^{7}\leq\leavevmode\nobreak\ M_{X}\leq\leavevmode\nobreak\ 10^{16}$ GeV. We use the parameters of DM allowed by the gamma–ray and neutrino limits to reveal the possible DM contribution to the anisotropy observables. This study complements our previous works [39, 26], where constraints on the heavy decaying DM lifetime were obtained using the high–energy gamma–ray and neutrino flux upper limits.

2 Cosmic ray flux from the dark matter decay

We consider DM consisting of scalar particles $X$ decaying through the primary channel $X\rightarrow q\bar{q}.$ This implies hadronisation and subsequent decay of unstable hadrons. The final products of the decay cascade are photons, protons, neutrino etc. In this study we are interested in the decay products that can contribute to the cosmic–ray flux anisotropy observable at Earth — that is photons and protons. We follow the method of Refs. [56, 57] in calculation of the decay spectra. The details were reviewed in our previous works [39, 26] along with the justification of the approximations used. Here we just list the key points. We assume that all quark flavors are coupled to $X$ similarly. The photon production is assumed to be saturated by the pion decays while the $\sim 10\%$ contribution from kaons and other mesons is neglected. We also neglect the electro–weak corrections to the decay spectrum and assume that final state particles momenta are distributed isotropically on average.

The spectra of photons and protons from the decay of particle of mass $M_{X}$ can be obtained by the DGLAP evolution [58, 59, 60, 61] of the low energy scale fragmentation functions. For this calculation we use the code provided by the authors of Ref. [56] that solves DGLAP equations numerically in the leading order of $\alpha(s)$. It assumes that all quark flavours are coupled to gluon similarly and implies the mixing of gluon fragmentation function with the quark singlet fragmentation function. The initial fragmentation functions parametrized on scale of $1$ GeV are taken from Ref. [62] and extrapolated to the interval $10^{-5}\leavevmode\nobreak\ \leq\leavevmode\nobreak\ x\leavevmode\nobreak\ \leq\leavevmode\nobreak\ 1$. The examples of photon and proton spectra from the decay of $X$ particles with various masses are shown in Fig. 1.

The large–scale anisotropy model predictions should be calculated for the total flux of the high–energy cosmic rays. It is known that the flux is dominated by the isotropic contribution of charged particles which we denote as $J_{\rm exp}(E)$, while the possible decay of the DM gives only a small anisotropic admixture denoted by $J_{\rm DM}(\delta,\alpha,E)$, where $\{\delta,\alpha\}$ is equatorial coordinates. The “experimental” flux is taken from Telescope Array [63] and Tibet [64] spectrum measurements. We choose these two experiments, since their fluxes are close to the average of all the experimental spectra reported in Ref. [65]. The uncertainty in anisotropy predictions due to discrepancy between spectra measured in different experiments is estimated in the next section. We parametrize the “experimental” flux in the following way:

[TABLE]

where $E_{-1}=4.0\cdot 10^{15}$ eV, $E_{0}=1.0\cdot 10^{17}$ eV, $E_{1}=5.2\cdot 10^{18}$ eV and $E_{2}=6.3\cdot 10^{19}$ eV corresponds to the energies of “knee”, “second knee”, “ankle” and GZK cut-off respectively and the normalization factor $J_{1}$ is taken at $10^{18}$ eV. The values of spectral indexes are: $\gamma_{-1}=2.72$, $\gamma_{0}=3.12$, $\gamma_{1}=3.23$, $\gamma_{2}=2.66$, $\gamma_{3}=4.65$.

In turn, the DM part of the flux consists of proton and photon contributions 2)2)2)We neglect the comparable neutrino flux because the sensitivity of EAS experiments to neutrino are at least two orders of magnitude smaller than that to photons and protons. .

[TABLE]

where G stands for “Galactic” and denotes contribution from the DM decay in Milky Way while EG stands for “Extra–Galactic” and denotes the contribution from the rest of Universe. Since here we discuss DM decay and not interesting in annihilation the anisotropic patterns related to the matter clustering in the extragalactic contribution are washed out [66] and this contribution can be considered as isotropic. For the purposes of anisotropy study it is convenient to consider the extragalactic DM contributions as a part of $J_{\rm exp}(E)$. Indeed, at low energies where $J_{\rm exp}$ dominates over $J_{\rm DM}$ the predicted anisotropy is small and the accounting of $J_{\rm DM}^{\rm EG}$ provides only few percent correction to it. Above the GZK threshold energy the $J_{\rm DM}$ starts to dominate, but extragalactic photon and proton DM fluxes are suppressed by attenuation effect in the same way as $J_{\rm exp}$ and only galactic contributions are relevant. These assumptions were justified by direct calculations of anisotropy with the actual parameters of the DM. Therefore, we take the total flux in the following form

[TABLE]

For the galactic flux calculation we use the Navarro-Frenk-White DM distribution [67, 68] with the parametrization for Milky Way from Ref. [66]3)3)3)For comparison we also calculate the resulting anisotropy using Burkert DM profile [69].. For galactic gamma–ray flux we take into account only prompt photon spectra of DM decay and neglect the smaller amount of photons from inverse Compton scattering (ICS) of prompt $e^{\pm}$ on the interstellar background photons. This assumption was discussed in our previous paper [39]. However, we allow for the modification of photon spectra due to interactions with CMB photons. This correction becomes important for the $E_{\gamma}\lesssim 10^{19}$ eV i.e. for $M_{X}\lesssim 10^{14}$ GeV. We use numerical code [70] which simulates development of electron-photon cascades on CMB driven by the chain of $e^{\pm}$ pair production and inverse Compton scattering. Although the code allows to calculate the flux of the cascade photons it doesn’t take into account deflections of $e^{\pm}$ by the halo magnetic field. Since electrons in the code propagate rectilinearly they produce less cascade photons. Therefore the calculated flux of photons should be considered as conservative lower bound. The code also includes attenuation of photons on extragalactic background light (EBL), though the effect of EBL is negligible on distances which we consider. The corrections related to the production of electro–magnetic cascades are important for the energies of photon lower than $\sim 10^{19}$ eV — above $10^{19}$ eV the correction is less then $1\%$. The comparison of cascading spectrum of photons with the injected spectrum both calculated for the decay of DM in Milky Way and propagated to Earth is presented in Fig. 2. In turn, the galactic proton contribution is affected by the galactic magnetic field, which deflects the protons and therefore washes out the anisotropy pattern.

The Milky Way magnetic field can be decomposed to regular (large scale) and random components [71]. The large-scale magnetic field obtained from Faraday rotation of pulsars and extragalactic sources is typically around 1.5-2$\mu$G. The total magnetic field in the Solar neighborhood is about 6 $\mu G$, which suggests presence of comparable random component. Towards the Galactic center, the magnetic field strength increases, reaching values 7.6-11.2 $\mu$G at a Galactocentric radius of 4 kpc. The magnetic field strength in the gaseous halo, or thick disk, is comparable to that in the disk, with an uncertainty of a factor 2-3. The Larmor radius of a particle with energy $E$ and electric charge $qe$ in a regular magnetic field is

[TABLE]

where $B_{\perp}$ is the field component perpendicular to the particle’s motion. The critical energy $E_{c}$ for protons in the Milky Way magnetic field i.e. the energy where the Larmor radius equals to the coherence length of the turbulent component, is estimated as $E_{c}\simeq 0.3$ EeV. The flux of protons with energies $E<E_{c}$ is completely isotropic due to randomization of their momenta directions by the turbulent magnetic field component, while protons with higher energies spiral around the regular component of the field. Below we conservatively assume that only protons with energies above $10^{19}$ eV contribute to the flux anisotropy. We assume rectilinear propagation of these protons and neglect possible contribution of lower energy protons to the anisotropy. We justify this approximation in the next section by the comparison of proton and photon contributions to the anisotropy.

3 Anisotropy analysis and results

Several observables are used in the literature to study the large–scale angular distribution of cosmic–rays. The most commonly used method is based on the cosmic ray flux harmonic analysis. Unfortunately, none of the currently running EAS experiments is observing the full sky. However the experiments with the full duty cycle cover some band on celestial sphere due to the Earth rotation. In this situation the appropriate solution is the Fourier analysis in right ascension, where the flux assumed to be the average in declination over the particular experiments field of view. The flux can be presented in the form:

[TABLE]

where

[TABLE]

Below we use the normalized coefficients, $a_{n}\equiv 2A_{n}/A_{0}$, $b_{n}\equiv 2B_{n}/A_{0}$. The observable commonly reported by the experiments is the normalized amplitude of the first harmonic:

[TABLE]

To obtain the theoretical prediction for this quantity one needs to take into account the effective exposure of the particular experiment which is given by [72, 52]:

[TABLE]

where $a$ is geographical latitude of the observatory, $\theta_{\text{max}}$ is the maximal zenith angle accessible for fully efficient observation in the experiment and $\alpha_{m}$ is given by

[TABLE]

[TABLE]

After inclusion of the exposure into the analysis we have for $a_{1}$:

[TABLE]

and for $b_{1}$:

[TABLE]

From expressions (13)-(14) it is easy to estimate the uncertainty of theoretically predicted values of $a_{1}$ and $b_{1}$ due to variation of the experimental flux $J_{\rm exp}(E)$. Generally, at lower energies $J_{\rm exp}(E)$ dominates over $J_{\rm DM}$, therefore the flux error linearly maps to the uncertainty of $a_{1}(E)$ and $b_{1}(E)$, while at the higher energies $J_{\rm DM}(E)$ starts to supersede the ”experimental” flux and the impact of its uncertainty on the resulting anisotropy decreases. However, at $E\simeq 10^{20}$ eV, where $J_{\rm DM}$ is still subdominant the experimental flux uncertainty is almost two orders of magnitude. 4)4)4)We estimate the uncertainty by comparison of the flux measurements in Telescope Array and Pierre Auger experiments. Therefore, the predictions for anisotropy above $E\simeq 10^{20}$ eV should be interpreted with these reservations.

The way to improve the sensitivity of the experiments to large–scale anisotropy is to consider summarized data of two experiments with fields of view jointly covering the whole celestial sphere. In that case one can expand the flux in spherical harmonics. This method has been adopted to study the anisotropy of the ultra-high energy cosmic rays by Pierre Auger and Telescope Array experiments in Ref. [52] 5)5)5)In recent studies of Pierre Auger [54] and IceCube[49] the advanced analysis techniques was used to extract the full–sky angular power spectrum from the partial sky data of these experiments. However, the restrictions on the anisotropy patterns that could be extracted by these techniques makes questionable its applicability to the DM decay anisotropy search.. The expansion of the flux into spherical harmonics has the form:

[TABLE]

with the coefficients defined as

[TABLE]

where integration goes over the whole celestial sphere and the exposure effects assumed to be eliminated from the resulting experimental quantities. The quantity analogous to $r_{n}$ is the angular power spectrum which is defined as

[TABLE]

Since here we are interested anisotropy only, below we redefine $a_{lm}\rightarrow\sqrt{4\pi}a_{lm}/a_{00}$, i.e. normalize coefficient to the monopole one.

The alternative approach for the DM signal search in the full–sky cosmic–ray map would be the fitting of the map with the signal plus background model, using the profile likelihood method. This technique was effectively employed by Fermi-LAT collaboration for the search of MeV–GeV gamma-ray flux of DM decay/annihilation origin [73]. This method would be in general more sensitive to the DM signal than the harmonic analysis. However, as the appropriate analysis tools are yet undeveloped, we leave this issue for the future studies.

For any mass $M_{X}$ the lifetime $\tau$ can be constrained using the upper–limits on the amplitudes (9, 17) or on the coefficients (16). We use the data from EAS-TOP [48], IceCube [49], KASCADE [50], KASCADE-Grande [51], Yakutsk [53] and Pierre Auger [55] experiments. All data is interpreted in terms of $r_{1}$ amplitude $C.L.=95\%$ upper-limits, except the KASCADE-Grande data which is presented as $C.L.=99\%$ upper-limit. We also employ the result of joint Telescope Array and Pierre Auger full–sky anisotropy study [52] presented in the form of separate upper–limits on $a_{lm}$ coefficients. We conservatively assume that all the anisotropy is given by the DM decay. To obtain the constraints we vary the DM lifetime $\tau$ until the amplitude (9) touches one of the experimental constraints from below. The procedure for $a_{lm}$ limits is similar but we use all the coefficients $a_{lm}$ for given energy bin. To obtain 95 % CL limit taking into account the respective statistical penalty we use $(0.95)^{\frac{1}{2l+1}}$ CL limits derived from all $(2l+1)$ coefficients $a_{lm}$. Since the values of $a_{lm}$ lie in the sign-changing band in each case we chose the edge of the band which has the sign similar to the respective theoretically predicted coefficient as a limit value.

The results are shown in Fig. 3. As one can see all the anisotropy constraints lie in the parameter area which is already excluded by the high–energy gamma–ray and neutrino based limits. Since the anisotropy and gamma–ray limits are set by the same experiments the above result indicates that EAS experiments are more sensitive to the DM decay photons than to the respective anisotropy. This fact has some interesting consequences that we discuss below.

Another important feature is that the anisotropy constraints are stronger for the higher energies. This fact can be understood if we notice that the background isotropic cosmic–ray flux grows faster with the decrease of energy than the precision of anisotropy measurements. As it was anticipated the full–sky constraints surpass that of the limited sky coverage. The surprising fact is that second harmonic of full–sky analysis bounds $\tau$ stronger than the first one, although the theoretically predicted amplitude of the former is generally smaller (see below). This is the effect of the incidentally small value of the upper limit for one of the coefficients which is not compensated by the statistical penalty. The constrains by the third harmonic which are not shown in the figure are weaker than those by the first and second harmonic. The small bump around $M_{X}=10^{14}$ GeV on each curve is due to the accounting of galactic proton flux along with the photons, its impact on result is expectedly small.

For further experimental analysis development it is worth to discuss the maximal expected large–scale anisotropy from DM decay allowed by the most recent gamma–ray constraints [39]. We obtain it for the particular experiments and for the full–sky measurements by fixing the value of $\tau$ which is marginally allowed by the gamma–ray limits. The results for the range of masses $M_{X}$ are shown in Figs. 4–5 together with the recent limits. For the individual experiments we calculate the desired anisotropy sensitivity in terms of $r_{1}$. Variation of $r_{1}$ with the choice of an experiment reflects the fact that the anisotropy observed in particular experiment depends on its field of view. We use Pierre Auger and Telescope Array for medium and high energies and IceCube for lower energies. Increase in sensitivity and expanding of the energy range is expected in Telescope Array with the planned deployment of the low energy extension TALE [74], which will allow to collect events with the energies down to $10^{16}$ eV. The high area extension TA$\times$4 [75, 74] should give significant increase of the statistics at higher energies. One should note that the Southern hemisphere based experiments — Pierre Auger and IceCube have an advantage in galactic anisotropy study over the Northern hemisphere based Telescope Array because of more convenient position relative to the Galactic Center. We see that significant increase in experimental sensitivity to the large–scale anisotropy is need to detect the maximal signal expected from the DM decay. The IceCube sensitivity should be increased by roughly two orders of magnitude, while for Auger the respective values are from $\sim 10^{4}$ to $\sim 20$ times at 1 EeV energy.

For the full–sky analysis we present the predictions in therms of the first two angular power spectrum coefficients $C_{1}$ and $C_{2}$. While the full–sky constraints shown in Fig. 3 were imposed using the sets of coefficients $a_{1m}$ and $a_{2m}$, the power spectrum $C_{l}$ reveals the overall sensitivity of certain harmonic to the respective theoretical anisotropy pattern. In this sense the predictions shown are conservative. From the Figs. 4–5 one can learn that the large energies area is most opportune for the DM decay anisotropy search, while at the energies of $\sim 10$ EeV the sensitivity need to be increased by at least several hundred times comparing to the recent searches.

4 Discussion

The obtained results indicate that current EAS experiments are more sensitive to the photons from DM decay than to the respective anisotropy. This occurs due to relatively the good hadron–photon primaries separation in EAS analysis and due to insufficient sensitivity of the ground based EAS experiments to the large scale anisotropy. A natural obstacle here is the necessity to combine the results of two experiments for the full–sky analysis. We should also note the connection between the anisotropy and gamma–ray signal. The large–scale anisotropy if observed at a particular energy not accompanied by the gamma–rays should be attributed to physics other than the DM decay, e.g. to the imprint of Large Scale Structure of the Universe [44, 45, 46, 47] or the anisotropic particle acceleration in the local cosmos [40, 41, 80]. In other words, until the gamma–rays of the respective energies are detected the DM signal should not interfere with the study of astrophysical large–scale anisotropy.

Some of the future experiments may be more sensitive to anisotropy than to gamma–ray DM signal. For instance the EUSO experiment is planned to have high sensitivity to anisotropy[81] while its ability of photon–hadron primaries separation is expected to be lower that in current experiments [82]. In Ref [79] the anisotropy detection prospects from the DM decay signal allowed by current photon limits were found favourable for EUSO experiment at ultra–high energies. At the same time planned photon sensitivity improvements in the currently running experiments - Pierre Auger and Telescope Array (see Ref.[83] for details) would make them even more effective in search for the signal of heavy decaying DM.

Acknowledgements

We would like to thank V. Rubakov, S. Troitsky and G. Rubtsov for helpful discussions. We are especially indebted to R. Aloisio, V. Berezinsky and M. Kachelriess for providing the numerical code solving the DGLAP equations. This work has been supported by the Russian Science Foundation grant 14-22-00161.