Gamma-ray bounds from EAS detectors and heavy decaying dark matter constraints

The very high energy Galactic $\gamma$-ray sky is partially opaque in the ($0.1-10$) PeV energy range. In the light of the recently detected high energy neutrino flux by IceCube, a comparable very high energy $\gamma$-ray flux is expected in any scenario with a sizable Galactic contribution to the neutrino flux. Here we elaborate on the peculiar energy and anisotropy features imposed upon these very high energy $\gamma$-rays by the absorption on the cosmic microwave background photons and Galactic interstellar light. As a notable application of our considerations, we study the prospects of probing the PeV-scale decaying DM scenario, proposed as a possible source of IceCube neutrinos, by extensive air shower (EAS) cosmic ray experiments. In particular, we show that anisotropy measurements at EAS experiments are already sensitive to $\tau_{\rm DM}\sim {\cal O}(10^{27})$~s and future measurements, using better gamma/hadron separation, can improve the limit significantly.

the associated gamma-ray flux should still be detectable at VHE, albeit with a significantly suppressed spectrum. This has been noted soon after the discovery, see for instance [12] or [13], and current γ-ray constraints are one argument disfavoring close-by discrete Galactic sources for the majority of the events.
The calculation of both the expected signal and the observational constraints is however more involved in the case of truly diffuse local sources associated with the IceCube data, such as in an astrophysical origin in the Galactic halo [14], or a decaying Dark Matter (DDM) origin, see [15][16][17][18][19] (see also [20][21][22][23] for some DM related interpretations of the IceCube data). On the one hand, the expected signal is usually anisotropic, at very least due to our off-set position with respect to the center of the Galactic Halo, so its detailed calculation requires at least a 2D modeling of the problem, and multiple numerical integrations, a complication that typically is not fully taken into account even in the most recent calculations [24]. On the other hand the observational bounds (the most constraining ones being the one published by CASA-MIA [25] and those reported in a proceeding by KASCADE [26]) are derived: a) based on observations of a limited portion of the sky; b) typically assuming a isotropic γ-ray flux.
The hypothesis b) is clearly incorrect, since even an incoming isotropic flux would acquire an anisotropy due to the anisotropic absorption. The hypothesis a) means that a proper use of these constraints would require to re-run ad hoc analyses by the collaborations, based on specific energy and angular-dependent templates, to be convoluted with the detector characteristics. This is a pity, since as we have argued in [16], for an important class of scenarios like the DDM ones, this kind of data are what currently comes the closest to an independent test of the hypothesis 1 .
Although the important role of EAS probes of this scenario has been discussed in the past (see e.g. [12,16,24]), here we revisit the calculation of the expected γ-ray flux, with a triple goal: i) To estimate more precisely the spectral and angular shape of a DDM signal, with state of the art treatment for the primary γ-ray absorption and the inverse Compton component. ii) To point out that due to the generically anisotropic nature of the VHE γ-ray component, even detectors with little or without gamma-hadron rejection capability should be able to put constraints on these contributions based merely on the expected anisotropy. iii) To motivate experimental collaborations to specifically constrain some angular-energy templates, to optimize their constraining power for specific models. For the DDM case, for instance, an intermediate step in this direction would be to derive (energy-dependent) bounds in coronas around the GC, in Galactic cylindrical coordinates. We also discuss the greatly improved potential of detectors with significant hadron vs. γ-ray rejection capabilities.
As a case study we consider throughout this paper the peculiarities of the expected γ-ray flux from DDM. Yet, similar considerations would apply to any other Galactic diffuse flux model (a few examples have been listed e.g. in [27]). We will consider both the prompt γ-ray flux and the flux from inverse Compton (IC) scattering of e ± off the ambient photon bath. Both contributions would be present also in other diffuse flux models: in the commonly considered astrophysical hadronic production of neutrinos, they are associated to prompt γ- 1 Note that a superficial look at Fermi-LAT isotropic gamma-ray background (IGRB) results might suggest that they are already very constraining, notably thanks to the few high-energy points. We stress here that they are unfortunately not robust with respect to the IGRB extraction procedure from the Extragalactic gammaray background (EGB). In more technical terms, the determination of the IGRB by the Fermi-LAT team does not take into account the uncertainty in the subtraction of the point sources contribution (very uncertain at high energies, relying on an extrapolation), which would constitute the dominant source of uncertainty in the last IGRB points. If the EGB is conservatively used, instead, the bounds are degraded, as can be seen in [24]. rays from π 0 decay and e ± from π ± decay, respectively. The specific features studied in this paper, mainly the inherent anisotropy due to absorption at Galactic scale and the peculiar profile of IC flux due to diffusion/losses of PeV e ± , would similarly provide powerful diagnostic tools in probing alternative diffuse flux models. The only differences would be in the initial spectra and the geometric distribution of the source term in the Galaxy/Galactic halo.
This article is structured as follows: in section 2 we describe the peculiar energy-angular dependence of the γ-ray flux absorption, and our computation of the γ-ray opacity. In section 3 we compare the expected γ-ray flux from DDM with current constraints from EAS experiments, as well as the diagnostic power of forthcoming experiments. The two components of the γ-ray flux from DDM, prompt and IC flux, are discussed in sections 3.1 and 3.2, respectively. Section 4 is devoted to the discussion of the expected anisotropy of the total cosmic ray flux. Finally, in section 5 we conclude.

Absorption of γ-rays at Galactic scale
The γ-ray flux in the approximate range 10 −2 ÷ 10 2 PeV will suffer attenuation in the Galaxy due to the pair production γγ → e − e + process onto photon baths: at the lower energies, starlight (SL) and infrared (IR) photons constitute important targets (mostly however for directions towards the inner Galaxy), while at ∼ PeV energies and above the homogeneous cosmic microwave background (CMB) is dominant. In the following we calculate the optical depth τ γγ for both CMB and SL+IR, for different incoming directions and energies.
For the technically simpler case of pair production on CMB photons, the optical depth for photons of energy E γ coming from a source at distance L can be calculated as (here and in the following, we use natural units with c = k B = 1) where σ γγ is the pair production cross section given by where α is the fine-structure constant, m e the electron mass and θ is the angle between the momenta of photons. The n CMB (ε) is the differential number density of CMB photons given by where T CMB = 2.348 × 10 −4 eV. By changing variable ε → ε c , where ε c = εE γ (1 − cos θ)/2 is the photon center of momentum energy, and performing the integral on θ, the expression for τ CMB γγ can be reduced to a single integral to be performed numerically, namely (2.5) Figure 1a shows the τ CMB γγ as function of E γ for three different values of L = 4 kpc, 8.3 kpc and 20 kpc. As can be seen, for a source of γ-ray at Galactic center (GC), at about L = 8.3 kpc, the absorption is ∼ 70% at E γ ∼ 2 PeV. Figure 1b shows the contour plot of exp[−τ CMB γγ ], as function of photon energy E γ and source distance L. The optical depth due to pair production on the SL+IR photon bath can be calculated similarly to eq. (2.1), with the extra complication that the integral along the line of sight is non-trivial, since the photon bath number density n SL+IR also depends on position x, and the optical depth also depends on the Galactic coordinates (b, l). In the approximation that the photon field is inhomogeneous but isotropic one can write where the line-of-sight parameter s is related e.g. to the cylindrical coordinates (r, z), with the origin at the GC, by r = R 2 + s 2 cos 2 b − 2sR cos b cos l and z = s sin b , where R 8.3 kpc is the distance of the Sun to the GC. The number densities of SL and IR photons have been extracted from the GALPROP code [28] and their energy densities for some representative positions are plotted in figure 2. Obviously, the CMB radiation field is homogenous and thus pervades the whole Galaxy uniformly, while the SL and IR components of radiation field are clearly position dependent: larger at GC and in the Galactic disk, decreasing rapidly by moving perpendicularly from Galactic disk, along the z direction.
The optical depth due to SL+IR photon bath for two different distances and various directions are shown in figure 3. It is clear that the absorption effect is relevant around energies of O(100) TeV, but only for directions towards the inner Galaxy (b l 0). The calculated optical depths in this section are consistent with the results reported in [3]. The effect of the total opacity of Galactic medium (i.e., τ γγ = τ CMB γγ + τ SL+IR γγ ) will be discussed in the following sections.

Prompt component
The prompt component of the Galactic γ-ray flux from DM decay in the direction (b, l) is given by where m DM and τ DM are respectively the DM mass and lifetime, and τ γγ is the total optical depth. ρ h is the density profile of DM particles in our Galaxy as a function of radial distance (in spherical coordinates) from the Galactic center, . For our fiducial model we adopt a Navarro-Frenk-White density profile [29] where r c 24 kpc is the critical radius and ρ h = 0.18 GeV cm −3 , which yields a DM density at the Solar System ρ = 0.39 GeV cm −3 [30]. The line-of-sight integration parameter s is related to radial distance via The dN γ /dE γ is the energy spectrum of photons produced in the decay of a DM particle (here obtained from the PYTHIA 8.2 [31], including the weak gauge boson radiation corrections as from [32]). To illustrate the typical spectra from DM decay, in figure 4 we plot the E γ dN/dE γ for various decay channels of a DM particle with m DM = 4 PeV. In a specific model of the DM (i.e, specific decay channels with branching ratios determined by the model) the spectrum of γ-ray can be obtained by the appropriate weighting of the spectra in figure 4. In this paper we adopt the scenario introduced in [33] where the heavy DM particle is a sterile neutrino with mass ∼ 4 PeV and lifetime ∼ 10 28 s, with the branching ratios of the decay channels given by where U i are the the elements of the neutrino mixing matrix (for details see [33]). It is shown in [16] that this scenario provides a reasonable fit to the energy distribution of IceCube neutrino data. However, let us emphasize that these choices of branching ratios and decay channels are not extremely constrained. In fact, as discussed in [15], any model with a sizable branching ratio into hard (leptonic) channels, with the remaining (even dominant) branching ratio into soft (hadronic and gauge bosons) channels, would provide decent fit to the IceCube data.

Including the inverse Compton component
An additional component of the γ-ray flux comes from the inverse Compton (IC) scattering of the electrons and positrons from DM decay, up-scattering mostly the CMB photons, which writes Figure 4. Spectrum of the γ-ray yield in various decay channels of a DM particle with m DM = 4 PeV. The red curve for DM → quarks, shows the average spectra for the DM decay to all the quark flavors. The spectra are obtained via PYTHIA 8.2 [31].
where is given in eq. (3.3), P IC is the IC power due to up-scattering on different photon backgrounds and dn e /dE e is the differential number density of e + plus e − at steady state. Although the IC flux reported in this article is calculated by taking into account the spatialdependent nature of energy losses and the effect of spatial diffusion (see appendix A for the details of the calculation), in order to grasp the main features of IC flux, in the following we pursue a simplified version of the calculation. At the energies of interest here the transport of electrons and positrons in our Galaxy is determined almost exclusively by the energy losses. Also, one realizes that for directions close to the Galactic plane and for realistic values of the Galactic magnetic field (i.e., ∼ few µG, with a profile that increases towards the inner Galaxy) synchrotron emission is the dominant energy loss mechanism, simply because the synchrotron emission is always quadratic in the electron energy and does not suffer the Klein-Nishina suppression of IC on SL and IR photon baths. Also, at high energies the IC power P IC is almost exclusively due to up-scattering of the CMB photons, and thus independent of the position. The position dependence of the energy loss coefficient, b = −dE e /dt, is more involved and traces the Galactic magnetic field profile. However, in the approximation in which the thin gaseous disk of the Galaxy is embedded in a thick diffusive halo permeated by a constant magnetic field, the loss coefficient b is independent of the position. In this approximation (which we checked to be accurate whenever the IC signal is non-negligible, see appendix A for details), one can write and the total γ-ray flux (i.e., sum of the prompt and IC components) can hence be written as dN e /dE e is the e ± energy spectrum from DM decay, obtained via PYTHIA 8.2 [31]. P IC can be calculated straightforwardly as reported in appendix A. Yet, the energy loss coefficient b still depends on the poorly known value of the magnetic field, B halo , permeating the thick halo which extends to several kpc away from the disk, and for the lack of better information we approximated it as constant.
In figure 5 we show the γ-ray flux from DM decay, assuming m DM = 4 PeV and τ DM = 10 28 s (chosen to be close to best-fit parameters; the flux scales inversely with m DM and τ DM ) and for the decay pattern with the branching ratios of the decay channels given by eq. (3.4), for different directions in Galactic coordinates. The solid curves depict the prompt flux, eq. (3.1), from GC (red, top), anti-GC (blue, bottom) and Galactic Pole (orange, intermediate). In each of these curves the dot-dashed curve deviating from the solid curve at higher energies shows the flux neglecting the absorption of γ-rays discussed in section 2. When comparing the expected γ-ray flux from DDM with the experimental bounds, the importance of accounting properly for the absorption of γ-ray on CMB photons is manifest, particularly at high energy. The dashed curves in figure 5 show the IC flux: the red (blue) dashed curves are the IC flux form GC (anti-GC) direction. The orange dashed curve shows the IC flux from the Galactic pole direction, with the assumption that the Galactic magnetic field only consists of the (thin disk) regular field, given in eq. (A.5); i.e., B halo = 0. The cyan, black and green dashed curves show the IC flux from the Galactic pole within the assumptions B halo = 0.5, 1 and 2 µG, respectively. Finally, the green and brown bar lines with arrows show, respectively, the upper limits on the γ-ray flux inferred by CASA-MIA [25] and KASCADE [26] experiments.
Let us elaborate on the various IC fluxes shown in figure 5: the IC component is clearly sub-leading with respect to the prompt emission for directions along the Galactic plane (note the red and blue dashed curves). However, this is not necessarily the case for the IC flux from the Galactic poles. The enhancement of the IC flux from the Galactic pole direction originates from the fact that for vertical directions the b coefficient drops fasters than the DM density along the line of sight integration of eq. (3.7). This enhancement is sizable if one assumes that the magnetic field exponentially decreases for vertical directions (with the profile of eq. (A.5))-dashed orange curve in figure 5-so that the IC flux can become comparable to the prompt flux towards the Galactic pole. However, as we have mentioned earlier, it is realistic to assume that a non-zero magnetic field permeates a thick halo to large distances, as consistent with the assumption that a charged cosmic ray population still propagate diffusively in a region several kpc away from the disk. The constant B halo is a toy-model representation of this field, and its effect on the IC flux can be seen by the cyan, black and green dashed curves. In all cases, the emission is suppressed with respect to the "unmagnetized halo" situation. The reason is that a growing B halo leads to a larger energy loss coefficient b, and thus more suppressed IC flux, since a growing fraction of energy is channeled into synchrotron. In conclusion, the IC flux from directions close to the Galactic plane (low-b) is quite robustly predicted to be small. The exact value of IC flux towards the Galactic poles is hard to predict due to the uncertain thickness and B-field strength of the magnetized halo, with the orange dashed curves in figure 5 providing a reasonable upper limit to this uncertain component.
It is worth noting that the CASA-MIA and KASCADE experiments would have already probed interesting parameter space for DM models, if they had accumulated significant exposure towards inner Galaxy, e.g. if they had been located in the Southern hemisphere. Unfortunately, their acceptance mostly peaks in regions far away from the GC and hence they would have been exposed to more modest fluxes, comparable to the orange curve in figure 5, insufficient to test the model even for optimistic IC expectations. To illustrate this point, in the following we briefly describe some notions on the geometrical acceptance of EAS experiments. An EAS is often classified as γ-like event, as opposed to a hadronic-like event, based on a significantly poorer muon content of the former shower with respect to the latter (at a fixed primary energy). Only for events which are not too inclined with respect to the vertical this separation can be done meaningfully, thus imposing a cut on maximum zenith angle of the shower. Assuming that the detector is continuously operational (i.e., the acceptance is uniform with respect to azimuth, or right ascension in equatorial coordinate), the geometrical acceptance efficiency ω of an EAS experiment located at the latitude λ as function of declination δ, can be written as [34] ω(δ) ∝ cos λ cos δ sin α m + α m sin λ sin δ , (3.9) where  figure 6 and, as can be seen, the limits relax by moving to the regions where the γ-ray flux from DM increases. In fact in the regions where these experiments are mostly sensitive, the expected flux from decaying DM is ∼ 3 × 10 −13 TeV cm −2 s −1 sr −1 , which is almost one order of magnitude below the KASCADE upper limit.

Anisotropy
Despite the fact that current EAS bounds are not yet constraining enough for the DDM explanations of IceCube events, the interesting parameter space appears within reach. Even relatively modest optimizations of current sensitivities might thus prove crucial. In fact, the main reason for the degradation of the bounds discussed in the previous section relates to the incorrect assumption that the gamma-ray flux is isotropic. In this section, we discuss to what extent one may turn that weakness into an opportunity, suggesting that anisotropy studies alone, even without shower property discrimination capabilities, might contribute to the constraints. EAS experiments in fact routinely measure cosmic ray anisotropy, albeit often only in terms of some "partial estimator" like the dipolar anisotropy (averaged with respect to right ascension). Let us define a characteristic "gamma-ray induced anisotropy" as where dΦ CR /dE is the total cosmic ray flux, taken from [35]. The anisotropy variable as defined in eq. (4.1) mainly arises from the prompt flux and the contribution of IC flux is negligible, not only because the IC is sub-leading but also since it is expected to be relatively more isotropic. An immediate constraint on DM lifetime can be obtained by requiring that a γ does not exceed the observed total anisotropy in cosmic rays, a. In practice, by requiring that in no energy bin a γ exceeds by more than two sigma the measured value of a we can obtain a conservative bound on the DM lifetime as τ DM > 2.5 × 10 27 s. The power of this observable is due to the fact that the intrinsic anisotropy in charged cosmic rays is at the level of 10 −4 ÷ 10 −3 , while a much larger (by two to three orders of magnitude!) relative anisotropy in gamma-rays is expected, at very least due to the off-center position of the Sun in the Galaxy. This means that, despite the fact that gamma-rays only constitute a small fraction of the overall CR flux at 0.1 − 1 PeV energies, in the anisotropy observable one can benefit from a larger signal to noise ratio. Accounting for absorption, however, suppresses the gamma-ray anisotropy, since pair-production is more severe in the GC direction than the anti-GC direction. In figure 7 the blue solid (dashed) curve shows the expected anisotropy a γ (without) taking into account the absorption, for the fiducial choice of lifetime discussed The blue solid (dashed) curve depict the anisotropy by taking into account (neglecting) the absorption of γ-rays, for τ DM = 10 28 s. The red dot-dashed curve shows the anisotropy for τ DM = 2.5 × 10 27 s which is the lower limit on lifetime at 2σ from anisotropy data. The data points show the measured anisotropies by EAS-TOP [36][37][38], Akeno [39], IceTop [40] and IceCube [41] experiments.
previously; while the red dot-dashed curve corresponds to the limiting value when a γ exceeds the measured a at 2σ. For comparison, we also report the amplitudes of dipolar anisotropies measured by different experiments. A few remarks are in order: • The suppression of anisotropy due to the absorption can be clearly seen. It also contributes to the peculiar energy dependence of a γ , decreasing with energy, while the observed anisotropy a moderately increases with energy.
• perhaps surprisingly, the bounds following from anisotropy are at least comparable in strength with the previously obtained bounds coming from comparisons with the (prompt) flux limits from EAS detectors and Fermi-LAT diffuse isotropic data, at the level of 10 27 s.
• The a γ observable, on the other hand, has a higher sensitivity to the inner Galaxy DM profile. For instance, the previously quoted bound of 2.5 × 10 27 s for the fiducial NFW profile would degrade to 1.9 × 10 26 s for a cored isothermal profile [42] with where r c = 4.38 kpc and ρ h = 1.387 GeV cm −3 .
One may wonder to what extent the above considerations are spoiled by a realistic account of experimental angular and energy resolutions. For instance, CASA-MIA had an angular resolution going from 2 • at low energies to better than 0.4 • at high-energies [43]. This is almost irrelevant for a cored DM profile, while a O(1 • ) resolution can degrade the bound for a NFW profile down to 6 × 10 26 s. Nonetheless, this constitutes a sub-leading uncertainty with respect to the one coming from the unknown shape of the DM profile in the inner Galaxy. Concerning energy resolution, it was demonstrated that CASA-MIA could detect spectral features not larger than 0.2 − 0.3 dex in energy [44]. Recent cross-calibrations between the energy scale retrieved via surface detectors and the one inferred via fluorescence light suggest that there might be an over-estimate on the absolute energy scale of the former experiments of the order of 30%, see for instance [45]. Fortunately, despite these uncertainties, figure 7 shows that the expected signal peak is very broad, changing by no more than ∼ 20% between 150 TeV and 350 TeV. We do not expect thus that accounting for energy resolution and energy scale uncertainty would degrade our conclusions by more than O(10%). It is also interesting that at few hundreds TeV the expected anisotropy from DDM matching IceCube observations may be only ∼ one order of magnitude below the measured overall dipolar anisotropy, while at higher energies (∼ PeV) the suppressed anisotropy is smaller by a factor of few, and its ratio to the charged cosmic ray signal is significantly less favorable. This suggests a first potential strategy to improve the constraints by using the energy-dependence and phase information of the anisotropy: although no deterministic prediction of the expected anisotropy due to charged cosmic rays is possible, one could calibrate a model for the stochastic phase fluctuation on the high-energy bins (say, 700 to 2000 TeV) where the charged cosmic ray contribution to a is definitely expected to dominate, and use the low energy band (around 200-300 TeV), where the fractional contribution of DM is expected to be maximal, to put a constraint on the amplitude of a sub-leading contribution to the total anisotropy due to a γ . The latter is characterized by a fixed direction (constant phase) and a specific energy dependence, very different from the competing charged cosmic ray anisotropy.
Despite some room for improvement, with past data one cannot really expect order-ofmagnitude gains in sensitivity. However, a yet more powerful way to improve over the present analysis would be to rely on the gamma/hadron discrimination possibility attained by the present generation of EAS detectors. In general, cuts based on the morphology of the shower (sometimes dubbed "compactness" criteria) allow one to select a photon-rich sample, keeping γ fraction of the initial photons, while retaining only h γ of the contaminating hadronic background. The ratio Q ≡ γ / √ h allows to quantify the gain in sensitivity to a photon signal when this cut is applied. While some rejection capability was already present in old experiments, even gamma-ray astrophysics oriented EAS detectors of the past generation, such as MILAGRO, were limited by h 0.05 − 0.1 [46] with corresponding Q-factors never much larger than ∼ 2. On the other hand, the situation is significantly different already at currently operating water Cherenkov EAS gamma observatories, such as HAWC. Such an experiment has similar energy resolution performances as the above-mentioned ones (about 40% at E γ > 10 TeV [46]), and even better angular resolution, about 0.1 • at E γ > 10 TeV [46]. But the major improvement is in the rejection capability of the hadron background: at highenergies, 99.9% of the background is routinely rejected, but stringent cuts with h 10 −4 and γ 25% above 10 TeV, i.e. Q 30, have already been reported [47], with even better performances that could be attained [48]. With an effective area, A eff , approaching ∼ 10 5 m 2 at high energy and a field of view of ∆Ω ∼ 2 sr, HAWC is expected to reach a sensitivity below the level of the IceCube diffuse neutrino flux, thus providing a unique constraint on the electromagnetic counterpart of the neutrino signal [47]. A high-energy photon-enriched sample in T = 1 year of HAWC data would consist of about against a number of background CR events h T ∆Ω A eff Φ CR (E CR > 10 TeV) 3 × 10 6 . This would certainly ease the measurement of the gamma-ray spectrum, as already noted in the past [12]. However, it is perhaps even more remarkable that the expected anisotropies of O(10%) in the gamma-ray sample correspond to variations in gamma-ray counts of ∼ 10 4 , as opposed to anisotropies in the CR background which are expected to be of ∼ 10 3 events. Put otherwise, in the gamma enriched sample the anisotropy should be fully dominated by the gamma-ray contribution and, with such statistics, HAWC may provide a crucial test of the DM hypothesis through anisotropy studies, besides spectral ones. The situation may be even more favorable with future detectors. HiSCORE [49] has two orders of magnitude larger effective area than HAWC at high-energy, but 1 Q 2 and thus is not ideal for this kind of measurement, although it may still be useful for complementary studies [12]. However LHAASO [50], thanks to its optimized hadron rejection capability, would provide the ultimate sensitivity for this type of analysis: According to [51], thanks to the KM2A array (for detecting hadronic induced muons) surrounding the 10 5 m 2 Cherenkov detector, at E γ > 80 TeV LHAASO would reach γ 1 and h 10 −7 , which ensures that observations would be essentially CR background-free.
In summary, anisotropy data offer a complementary tool to constrain DDM contribution to the gamma-ray flux at 0.1−2 PeV energy. A simple analysis shows that current constraints from the normalization of the anisotropy are competitive with the other methods, and we sketched two possible approaches to improve upon them: with a reanalysis of current data, the addition of phase information could already allow one to achieve sensitivity to a subleading DM induced anisotropy. The optimal energy window for DM signal to CR noise appears to be at ∼ 200 − 300 TeV, while higher energies should be dominated by the CR component, and could be used to calibrate the dominating background anisotropy, whose phase is fluctuating with energy. A major improvement relies however on the current (HAWC) and forthcoming (LHAASO) generation of EAS gamma-ray detectors: thanks to their greatly enhanced photon/hadron rejection capabilities, already in HAWC the anisotropy signal should be dominated by the gamma-ray one, allowing for a first stringent test of any "Galactic based" model for the origin of a significant fraction of IceCube events. In LHAASO, we expect the sample to be essentially background-free, allowing for the ultimate test (or detailed studies, in case of positive detection) of these models. Needless to say, these experiments have a great potential for heavy dark matter constraints, which has only recently started to be studied.

Conclusions
The IceCube discovery of a PeV flux of astrophysical neutrinos has several implications for high energy astrophysics and astroparticle physics, as proven by the broad community whose attention has already triggered. If the sources of (a fraction of) these events are Galactic, this discovery paves the way to "Galactic γ-ray PeVatrons". In this regard it is timely to investigate in more detail the peculiarities of VHE γ-ray propagation in our Galaxy. In fact, the yet undisclosed 0.1 − 10 PeV γ-ray window in Galactic astrophysics would be affected by the pair production absorption. In this article, we discussed some effects that this phenomenon has onto expected signals in extensive air shower (EAS) detectors. We selected the benchmark case of a continuous emission from the Galactic halo in a decaying dark matter (DDM) scenario, although most of our results apply mutatis mutandis to other source distributions. Our choice was also motivated by the observation that, while the potential of these EAS instruments for DM detection has sometimes been considered in the past, see for instance [52,53], their potential for serendipitous DM discoveries in more unconventional scenarios has passed mostly unnoticed.
In this article, we have calculated the expected γ-ray flux from DDM, with values for the mass and lifetime motivated by the IceCube flux, by considering both the prompt and IC scattering components. While the angular dependence of the prompt γ-ray flux is robust (with the detailed features of its energy spectrum depending however on poorly known branching ratios, and only computable within the theoretical uncertainties unavoidably affecting PeV scale physics), the IC scattering component can potentially exhibit unusual angular features. Indeed, the IC scattering profile depends on the properties of the magnetized halo of our Galaxy, which is typically poorly known at large vertical distances from Galactic disk. By simple modeling of the Galactic magnetized halo, we have calculated the expected IC flux from directions near Galactic poles, arguing that in the most optimistic case the IC flux can be enhanced, becoming comparable to the prompt flux from this direction. However, a relatively large halo magnetic field at high latitudes will suppress the IC flux significantly.
Typical EAS bounds on the γ-ray fraction in cosmic ray flux have been derived under the hypothesis of a isotropic flux. We argued that this approximation is untenable in the energy range of interest, even in the limit of an isotropically emitted flux, due to the direction (and energy) dependent optical depth of the Galactic sky. We quantified this observation showing that the effect is so relevant that the current constraints from KASCADE or CASA-MIA experiments-which naively would appear to constrain some DDM scenarios-are still at least several times too weak. Ideally, an exposure only marginally better than the above mentioned ones, but at an observatory located in the Southern hemisphere, would be much more promising in this respect. We also argued that anisotropy may offer an independent handle to constrain DDM (as well as other similar scenarios): the expected γ-ray flux induces an anisotropy in the overall cosmic ray flux only a few times smaller at ∼ O(100) TeV (and about one order of magnitude at ∼ PeV) than the current measurements of the dipolar anisotropy routinely performed in EAS experiments. Turning the argument around, existing data are already sensitive to DM lifetimes of O(10 27 ) s, only one order of magnitude away from the value needed to fit IceCube events (∼ 10 28 s), showing the power of anisotropy analyses and motivating an attempt to improve over the current situation.
Some progress is expected from the experimental point of view. For instance, the IceTop facility at the top of the IceCube detector can look at the Southern sky. Unfortunately, while located at the South Pole, the GC is not in standard analyses involving the IceTop array: the IceCube detector plays the role of penetrating muon detector for IceTop facility, which requires that the axis of air shower should pass through the volume of IceCube. This requirement leads to a cut on the zenith angle of shower ∼ 30 • for IC40 configuration [54], with the possibility of increasing to ∼ 45 • for the whole IC86 configuration. At the South Pole, the GC is located at zenith angle ∼ 61 • . Although the expected sensitivity of IceTop, after 5 years of data taking, is close to the CASA-MIA limit at ∼ PeV (see figure 14 at [54]), due to the closeness of the field of view to GC, IceTop can still moderately improve the limits. Both for past and forthcoming data, we argued that a dedicated analysis might be sensitive to sub-leading anisotropies expected from DDM, notably if the shape of the anisotropy and its energy dependence, whose expectations are relatively well-known, are imposed a priori in the analysis.
Eventually, however, we have argued that greatly improved photon/hadron discrimination capabilities are needed for a decisive jump in the sensitivity in both existing and forthcoming experiments: The recently inaugurated HAWC observatory [52] located in Mexico (with latitude λ ∼ 19 • and zenith angle cut θ m ∼ 45 • ) thanks to quality factors Q = γ / √ h ∼ 30 could provide first crucial tests of the DDM scenario, not only via spectral studies (see [12]) but particularly when adding angular information, as discussed in our article. A future experiment such as LHAASO [50], benefiting from the KM2A array, is expected to provide a gamma-enriched data set which is almost background (CR)-free, paving the way to exquisite constraints or, in case of detection, to detailed studies of the spectrum and morphology of the signal.
Definitely, the opening of the PeV astrophysical window may offer new opportunities for interesting multi-messenger studies, probably shedding light on intriguing astroparticle questions. As already noticed in the past, and as we further argued here, searches in this new window significantly benefit from EAS experiments. Note that, while the low energy part of the neutrino flux observed by IceCube (recently extended down to ∼ 10 TeV [55]) can naturally receive contribution from Galactic sources, perhaps easing their identification, pinpointing the origin of the high energy part of the flux is more challenging. In that respect, EAS experiments appear a unique and powerful probe. For such a task, as illustrated in this article, rather than considering the Galactic-scale horizon imposed by the finite optical depth as a limitation, we should perhaps reconsider it as an original opportunity to exploit the specific capabilities of EAS detectors.
where σ T is the Thomson cross section, Γ ε = 4εγ/m e and γ = E e /m e ; and n(ε, x) is the SL+IR+CMB photons differential number density at position x. The energy loss due to synchrotron radiation can be calculated by where B is the magnitude of the total magnetic field in our Galaxy, consisting of regular and turbulent components 2 . For the regular magnetic field we adopt the following profile [56] where R = 8.3 kpc, r B = 10 kpc, z B = 2 kpc and B 0 = 4.78 µG. For the halo (possibly turbulent) magnetic field we assume a uniform constant strength magnetic field. Figure 8 shows the energy loss function b, as function of E e , at GC (black solid), Sun position (blue solid) and at vertical distance z = 5 kpc from Galactic plane (at GC) for three different assumptions for B halo = 0, 1 and 2 µG, respectively by solid, dashed and dot-dashed red curves. As can be seen by increasing the value of B halo from 0 to 2 µG, the energy loss coefficient at z = 5 kpc increases by about one order of magnitude, which justify the suppression of IC flux (black and green dashed curves in figure 5) for larger B halo . Obviously, the effect of halo field is smaller at lower energies, since the synchrotron emission is the main mechanism of energy loss at higher energies. Diffusion halo function I diff (E e , E e , x): The diffusion halo function can be calculated by solving the diffusion-loss equation of e ± in the Galaxy. To avoid repetition we skip reporting the details of calculation, which is done according to the prescription reported in [57]. However, it is worth mentioning that at high energies which we are interested in this paper I diff 1, and so the results of of IC flux reported in figure 5 are only marginally dependent on the diffusion halo function. Put otherwise, the approximation described in the main text is actually excellent.
The other ingredient in the calculation of IC flux is the IC power P IC (see eq. (3.5)) which can be decomposed into the IC power from each component of the photon bath; i.e., P IC = i P i IC , where P i IC is the IC power from the photon bath n i with i = CMB and SL+IR. The P i IC can be written as Figures 9a and 9b, respectively, show the IC powers P CMB IC and P SL+IR IC as function of E γ for E e = 1, 10, 10 2 , 10 3 , 10 4 TeV at GC. As can be seen, at high energies the IC power sharply peaks at E e . Namely, in the IC scattering almost all the e ± energy transfers to the photon. Also, by comparing the corresponding curves in figures 9a and 9b, it can be seen that at high energies the main contribution to the total IC power comes from P CMB IC . The IC power in figure 9a is independent of the position in Galaxy (due to the uniform CMB photon bath), while the IC power due to SL+IR strongly decreases by distancing from the GC, especially in the vertical direction with respect to Galactic plane.    (E γ , E e , x = 0), for E e = 1, 10, 10 2 , 10 3 , 10 4 TeV, from the left to the right, respectively.