Gamma-ray bounds from EAS detectors and heavy decaying dark matter constraints
The very high energy Galactic -ray sky is partially opaque in the () PeV energy range. In the light of the recently detected high energy neutrino flux by IceCube, a comparable very high energy -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 -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 s and future measurements, using better gamma/hadron separation, can improve the limit significantly.
a]Arman Esmaili b]and Pasquale Dario Serpico Prepared for submission to JCAP \trfontLAPTH-026/15
Gamma-ray bounds from EAS detectors and heavy decaying dark matter constraints
INFN, Laboratori Nazionali del Gran Sasso, Assergi (AQ), Italy
LAPTh, Univ. de Savoie Mont Blanc, CNRS, B.P.110, Annecy-le-Vieux F-74941, France
The very high energy (VHE) -ray band, loosely defined as the one above (100)GeV, is peculiar in many respects. For instance, on the detection side, it is at the upper edge of the energy range which can be probed from space via direct techniques, and can be effectively studied only by ground-based telescopes. An important theoretical feature in the study of VHE -rays is that absorption is relevant, i.e. the extragalactic sky is not transparent to VHE -rays due to the pair-production absorption onto the Extragalactic Background Light (EBL) and the cosmic microwave background (CMB). Actually, for photon energies of TeV even the Galactic diffuse radiation field starts playing a role, while at PeV energies the mean-free path of -rays due to the absorption on CMB photon bath approaches kpc, comparable to the distance to Galactic Center (GC). Hence, for VHE -rays even the Galactic sky becomes partially opaque. While the opacity of the extragalactic sky is routinely taken into account in the analysis of cosmologically distant active galactic nuclei, both by the Fermi-LAT space telescope  and by the ground based instruments, see e.g. , the absorption at Galactic scale is usually neglected, although the theoretical importance of the phenomenon has been acknowledged in the past, see e.g. .
The most obvious reason for the relative lack of interest in this phenomenon is that the “PeVatron” window in -ray astrophysics has yet to be opened, with current telescopes running out of statistics below (100) TeV (even for the most powerful sources) to be sensitive to spectral distortions. Another reason may simply be that PeVatrons (in gamma-rays) in our Galaxy are expected to be rare, if present at all, requiring extreme acceleration parameters. The situation has however changed in the last couple of years, following the discovery by the IceCube collaboration of a diffuse neutrino flux in the 30 TeV-2 PeV range [4, 5, 6]. Virtually in any conceivable model for its origin, one expects an associated -ray flux with similar energy budget in the PeV energy range. The main unknown concerns the distance at which the sources of these events are located: if the neutrino flux is due to a population at cosmological distances, the absorption should be so severe that the initial VHE flux would cascade down below (100)GeV, contributing to the residual isotropic background precisely measured by Fermi . However, if it is due to Galactic or nearby extragalactic sources (see e.g. [8, 9, 10, 11]), 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  or , 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 , 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 . On the other hand the observational bounds (the most constraining ones being the one published by CASA-MIA  and those reported in a proceeding by KASCADE ) 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 , 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 111Note 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 gamma-ray 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 ..
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 ). We will consider both the prompt -ray flux and the flux from inverse Compton (IC) scattering of 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 -rays from decay and 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 , 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.
2 Absorption of -rays at Galactic scale
The -ray flux in the approximate range PeV will suffer attenuation in the Galaxy due to the pair production 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 coming from a source at distance can be calculated as (here and in the following, we use natural units with )
where is the pair production cross section given by
is the fine-structure constant, the electron mass and is the angle between the momenta of photons. The is the differential number density of CMB photons given by
where eV. By changing variable , where is the photon center of momentum energy, and performing the integral on , the expression for can be reduced to a single integral to be performed numerically, namely
Figure (a)a shows the as function of for three different values of kpc, 8.3 kpc and 20 kpc. As can be seen, for a source of -ray at Galactic center (GC), at about kpc, the absorption is at PeV. Figure (b)b shows the contour plot of , as function of photon energy and source distance .
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 also depends on position x, and the optical depth also depends on the Galactic coordinates . In the approximation that the photon field is inhomogeneous but isotropic one can write
where the line-of-sight parameter is related e.g. to the cylindrical coordinates , with the origin at the GC, by
where is the distance of the Sun to the GC. The number densities of SL and IR photons have been extracted from the GALPROP code  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 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 TeV, but only for directions towards the inner Galaxy (). The calculated optical depths in this section are consistent with the results reported in . The effect of the total opacity of Galactic medium (i.e., ) will be discussed in the following sections.
3 Expected fluxes from DDM and constraints from EAS experiments
3.1 Prompt component
The prompt component of the Galactic -ray flux from DM decay in the direction is given by
where and are respectively the DM mass and lifetime, and is the total optical depth. 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 
where is the critical radius and , which yields a DM density at the Solar System . The line-of-sight integration parameter is related to radial distance via
The is the energy spectrum of photons produced in the decay of a DM particle (here obtained from the PYTHIA 8.2 , including the weak gauge boson radiation corrections as from ). To illustrate the typical spectra from DM decay, in figure 4 we plot the for various decay channels of a DM particle with 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  where the heavy DM particle is a sterile neutrino with mass PeV and lifetime s, with the branching ratios of the decay channels given by
where are the the elements of the neutrino mixing matrix (for details see ). It is shown in  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 , 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.
3.2 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
where is given in eq. (3.3), is the IC power due to up-scattering on different photon backgrounds and is the differential number density of plus at steady state. Although the IC flux reported in this article is calculated by taking into account the spatial-dependent 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 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, , 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 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
is the energy spectrum from DM decay, obtained via PYTHIA 8.2 . can be calculated straightforwardly as reported in appendix A. Yet, the energy loss coefficient still depends on the poorly known value of the magnetic field, , 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 PeV and s (chosen to be close to best-fit parameters; the flux scales inversely with and ) 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., . The cyan, black and green dashed curves show the IC flux from the Galactic pole within the assumptions , 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  and KASCADE  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 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 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 leads to a larger energy loss coefficient , 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-) 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 -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 
where is the maximum zenith angle acceptance of the detector. For the CASA-MIA experiment , and for KASCADE . The blue shaded regions in figures (a)a and (b)b show the estimated acceptance of CASA-MIA and KASCADE experiments, respectively, in Galactic coordinates. Darker blue regions correspond to higher effective exposure to a -ray flux. The superimposed red curves on the plots are the contours of at PeV, with values down to , respectively from the inner to outer circles (in units TeV cm s sr), from the decay of DM with PeV, s and branching ratios given in eq. (3.4). From figure 5, it follows that the upper limits on at PeV, are and respectively from CASA-MIA and KASCADE experiments. However, these upper limits apply to the dark blue regions of 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 TeV cm s sr, which is almost one order of magnitude below the KASCADE upper limit.
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 is the total cosmic ray flux, taken from . 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 does not exceed the observed total anisotropy in cosmic rays, . In practice, by requiring that in no energy bin exceeds by more than two sigma the measured value of we can obtain a conservative bound on the DM lifetime as s. The power of this observable is due to the fact that the intrinsic anisotropy in charged cosmic rays is at the level of , 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 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 (without) taking into account the absorption, for the fiducial choice of lifetime discussed previously; while the red dot-dashed curve corresponds to the limiting value when exceeds the measured at . 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 , decreasing with energy, while the observed anisotropy 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 s.
The observable, on the other hand, has a higher sensitivity to the inner Galaxy DM profile. For instance, the previously quoted bound of s for the fiducial NFW profile would degrade to s for a cored isothermal profile  with where and .
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 at low energies to better than 0.4 at high-energies . This is almost irrelevant for a cored DM profile, while a resolution can degrade the bound for a NFW profile down to 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 dex in energy . 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 . Fortunately, despite these uncertainties, figure 7 shows that the expected signal peak is very broad, changing by no more than 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 . 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 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 . 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-of-magnitude 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 of the contaminating hadronic background. The ratio 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  with corresponding -factors never much larger than . 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 at TeV ), and even better angular resolution, about at TeV . But the major improvement is in the rejection capability of the hadron background: at high-energies, 99.9% of the background is routinely rejected, but stringent cuts with and above 10 TeV, i.e. , have already been reported , with even better performances that could be attained . With an effective area, , approaching m at high energy and a field of view of 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 . A high-energy photon-enriched sample in year of HAWC data would consist of about
against a number of background CR events . This would certainly ease the measurement of the gamma-ray spectrum, as already noted in the past . However, it is perhaps even more remarkable that the expected anisotropies of in the gamma-ray sample correspond to variations in gamma-ray counts of , as opposed to anisotropies in the CR background which are expected to be of 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  has two orders of magnitude larger effective area than HAWC at high-energy, but and thus is not ideal for this kind of measurement, although it may still be useful for complementary studies . However LHAASO , thanks to its optimized hadron rejection capability, would provide the ultimate sensitivity for this type of analysis: According to , thanks to the KM2A array (for detecting hadronic induced muons) surrounding the Cherenkov detector, at TeV LHAASO would reach and , 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 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 sub-leading DM induced anisotropy. The optimal energy window for DM signal to CR noise appears to be at 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.
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 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 [53, 52], 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 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 s, only one order of magnitude away from the value needed to fit IceCube events ( 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 for IC40 configuration , with the possibility of increasing to for the whole IC86 configuration. At the South Pole, the GC is located at zenith angle . 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 ), 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  located in Mexico (with latitude and zenith angle cut ) thanks to quality factors could provide first crucial tests of the DDM scenario, not only via spectral studies (see ) but particularly when adding angular information, as discussed in our article. A future experiment such as LHAASO , 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 TeV ) 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.
We thank M. Cirelli, P. Panci, M. Taoso and F. Vissani for useful discussions, and P. Blasi for discussions and comments on a preliminary version of this manuscript. A. E. thanks Sebastian Wild for several helps on PYTHIA. At LAPTh, this activity was developed coherently with the research axes supported by the Labex grant ENIGMASS. For A. E. this research was partially supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. The author would like to express a special thanks to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support.
Appendix A Details of the IC flux calculation
In this appendix we provide the details of the IC flux calculation, which have been used to compute the curves in figure 5. As we discussed, the IC flux is given by eq. (3.5) and here we discuss calculation of the ingredients and .
The steady state spectrum of at position in our Galaxy (we will in fact assume azimuthal symmetry in cylindrical coordinates) can be calculated by
where is the energy loss coefficient and is the function taking into account the diffusion of . Let us elaborate on each of these ingredients:
Energy loss function : The lose energy by two main processes: synchrotron radiation in the magnetic field of Galaxy and energy loss due to IC scattering on the photon bath (SL+IR and CMB). Since both the magnetic field and SL+IR field are position dependent, the energy loss function also is position dependent. So,
where is the Thomson cross section, and ; and is the SL+IR+CMB photons differential number density at position . The energy loss due to synchrotron radiation can be calculated by
where is the magnitude of the total magnetic field in our Galaxy, consisting of regular and turbulent components 222The interstellar magnetic field of our Galaxy off the Galactic plane is determined by the Faraday-rotation measurement of the polarization of extragalactic radio sources. However, these measurements depend on the assumed free-electron distribution at high latitudes which is poorly known. On the other hand, even in the Galactic plane, the estimated turbulent magnetic field is a factor of few larger than the regular magnetic field.. For the regular magnetic field we adopt the following profile 
where kpc, kpc, kpc and G. For the halo (possibly turbulent) magnetic field we assume a uniform constant strength magnetic field. Figure 8 shows the energy loss function , as function of , at GC (black solid), Sun position (blue solid) and at vertical distance kpc from Galactic plane (at GC) for three different assumptions for and 2 G, respectively by solid, dashed and dot-dashed red curves. As can be seen by increasing the value of from 0 to 2 G, the energy loss coefficient at 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 . 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 : The diffusion halo function can be calculated by solving the diffusion-loss equation of in the Galaxy. To avoid repetition we skip reporting the details of calculation, which is done according to the prescription reported in . However, it is worth mentioning that at high energies which we are interested in this paper , 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 (see eq. (3.5)) which can be decomposed into the IC power from each component of the photon bath; i.e., , where is the IC power from the photon bath with CMB and SL+IR. The can be written as
where , and
Figures (a)a and (b)b, respectively, show the IC powers and as function of for TeV at GC. As can be seen, at high energies the IC power sharply peaks at . Namely, in the IC scattering almost all the energy transfers to the photon. Also, by comparing the corresponding curves in figures (a)a and (b)b, it can be seen that at high energies the main contribution to the total IC power comes from . The IC power in figure (a)a 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.
-  M. Ackermann et al. [Fermi-LAT Collaboration], “The Imprint of The Extragalactic Background Light in the Gamma-Ray Spectra of Blazars,” Science 338, 1190 (2012) [arXiv:1211.1671 [astro-ph.CO]].
-  A. Abramowski et al. [HESS Collaboration], “Measurement of the extragalactic background light imprint on the spectra of the brightest blazars observed with H.E.S.S,” Astron. Astrophys. 550, 4 (2013) [arXiv:1212.3409 [astro-ph.HE]].
-  I. V. Moskalenko, T. A. Porter and A. W. Strong, “Attenuation of vhe gamma rays by the milky way interstellar radiation field,” Astrophys. J. 640, L155 (2006) [astro-ph/0511149].
-  M. G. Aartsen et al. [IceCube Collaboration], “First observation of PeV-energy neutrinos with IceCube,” Phys. Rev. Lett. 111, 021103 (2013) [arXiv:1304.5356 [astro-ph.HE]].
-  M. G. Aartsen et al. [IceCube Collaboration], “Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,” Science 342, 1242856 (2013) [arXiv:1311.5238 [astro-ph.HE]].
-  M. G. Aartsen et al. [IceCube Collaboration], “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett. 113, 101101 (2014) [arXiv:1405.5303 [astro-ph.HE]].
-  M. Ackermann et al. [Fermi-LAT Collaboration], “The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV,” Astrophys. J. 799, no. 1, 86 (2015) [arXiv:1410.3696 [astro-ph.HE]].
-  S. Razzaque, “The Galactic Center Origin of a Subset of IceCube Neutrino Events,” Phys. Rev. D 88, 081302 (2013) [arXiv:1309.2756 [astro-ph.HE]].
-  C. Lunardini, S. Razzaque, K. T. Theodoseau and L. Yang, “Neutrino Events at IceCube and the Fermi Bubbles,” Phys. Rev. D 90, no. 2, 023016 (2014) [arXiv:1311.7188 [astro-ph.HE]].
-  C. Lunardini, S. Razzaque and L. Yang, “Multimessenger study of the Fermi bubbles: Very high energy gamma rays and neutrinos,” Phys. Rev. D 92, no. 2, 021301 (2015) [arXiv:1504.07033 [astro-ph.HE]].
-  M. Ahlers, Y. Bai, V. Barger and R. Lu, “Galactic TeV-PeV Neutrinos,” arXiv:1505.03156 [hep-ph].
-  M. Ahlers and K. Murase, “Probing the Galactic Origin of the IceCube Excess with Gamma-Rays,” Phys. Rev. D 90, no. 2, 023010 (2014) [arXiv:1309.4077 [astro-ph.HE]].
-  A. D. Supanitsky, “Gamma rays and neutrinos from a cosmic ray source in the Galactic Center region,” Phys. Rev. D 89, no. 2, 023501 (2014) [arXiv:1312.7304 [astro-ph.HE]].
-  A. M. Taylor, S. Gabici and F. Aharonian, “Galactic halo origin of the neutrinos detected by IceCube,” Phys. Rev. D 89, no. 10, 103003 (2014) [arXiv:1403.3206 [astro-ph.HE]].
-  A. Esmaili and P. D. Serpico, “Are IceCube neutrinos unveiling PeV-scale decaying dark matter?,” JCAP 1311, 054 (2013) [arXiv:1308.1105 [hep-ph]].
-  A. Esmaili, S. K. Kang and P. D. Serpico, “IceCube events and decaying dark matter: hints and constraints,” JCAP 1412, no. 12, 054 (2014) [arXiv:1410.5979 [hep-ph]].
-  A. Bhattacharya, M. H. Reno and I. Sarcevic, “Reconciling neutrino flux from heavy dark matter decay and recent events at IceCube,” JHEP 1406, 110 (2014) [arXiv:1403.1862 [hep-ph]].
-  C. Rott, K. Kohri and S. C. Park, “Superheavy dark matter and IceCube neutrino signals:bounds on decaying dark matter,” arXiv:1408.4575 [hep-ph].
-  C. S. Fong, H. Minakata, B. Panes and R. Z. Funchal, “Possible Interpretations of IceCube High-Energy Neutrino Events,” JHEP 1502, 189 (2015) [arXiv:1411.5318 [hep-ph]].
-  J. Zavala, “Galactic PeV neutrinos from dark matter annihilation,” Phys. Rev. D 89, 123516 (2014) [arXiv:1404.2932 [astro-ph.HE]].
-  A. Bhattacharya, R. Gandhi and A. Gupta, “The Direct Detection of Boosted Dark Matter at High Energies and PeV events at IceCube,” JCAP 1503, no. 03, 027 (2015) [arXiv:1407.3280 [hep-ph]].
-  J. Kopp, J. Liu and X. P. Wang, “Boosted Dark Matter in IceCube and at the Galactic Center,” JHEP 1504, 105 (2015) [arXiv:1503.02669 [hep-ph]].
-  Y. Daikoku and H. Okada, “PeV Scale Right Handed Neutrino Dark Matter in Flavor Symmetric extra U(1) model,” Phys. Rev. D 91, no. 7, 075009 (2015) [arXiv:1502.07032 [hep-ph]].
-  K. Murase, R. Laha, S. Ando and M. Ahlers, “Testing the Dark Matter Scenario for PeV Neutrinos Observed in IceCube,” Phys. Rev. Lett. 115, 071301 (2015) [arXiv:1503.04663 [hep-ph]].
-  M. C. Chantell et al. [CASA-MIA Collaboration], “Limits on the isotropic diffuse flux of ultrahigh-energy gamma radiation,” Phys. Rev. Lett. 79, 1805 (1997) [astro-ph/9705246].
-  G. Schatz et al. [KASCADE Collaboration], “Search for Extremely High Energy Gamma Rays with the KASCADE Experiment” Proceedings of the 28th ICRC, Tsukuba, Japan (2003).
-  O. E. Kalashev and S. V. Troitsky, “IceCube astrophysical neutrinos without a spectral cutoff and 10Ð10 eV cosmic gamma radiation,” JETP Lett. 100, no. 12, 761 (2015) [arXiv:1410.2600 [astro-ph.HE]].
-  http://galprop.stanford.edu/
-  J. F. Navarro, C. S. Frenk and S. D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490 (1997) 493 [astro-ph/9611107].
-  R. Catena and P. Ullio, “A novel determination of the local dark matter density,” JCAP 1008 (2010) 004 [arXiv:0907.0018 [astro-ph.CO]].
-  T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna and S. Prestel et al., “An Introduction to PYTHIA 8.2,” Comput. Phys. Commun. 191, 159 (2015) [arXiv:1410.3012 [hep-ph]].
-  J. R. Christiansen and T. Sjöstrand, “Weak Gauge Boson Radiation in Parton Showers,” JHEP 1404, 115 (2014) [arXiv:1401.5238 [hep-ph]].
-  T. Higaki, R. Kitano and R. Sato, “Neutrinoful Universe,” JHEP 1407, 044 (2014) [arXiv:1405.0013 [hep-ph]].
-  P. Sommers, “Cosmic ray anisotropy analysis with a full-sky observatory,” Astropart. Phys. 14, 271 (2001) [astro-ph/0004016].
-  T. Antoni et al. [KASCADE Collaboration], “KASCADE measurements of energy spectra for elemental groups of cosmic rays: Results and open problems,” Astropart. Phys. 24, 1 (2005) [astro-ph/0505413].
-  M. Aglietta et al. [EAS-TOP Collaboration], “A Measurement of the solar and sidereal cosmic ray anisotropy at E(0) approximates 10**14-eV,” Astrophys. J. 470, 501 (1996).
-  M. Aglietta et al., Proc. 28th Int. Cosmic Ray Conf., Tsukuba, 2003, 1, 183
-  M. Aglietta et al. [EAS-TOP Collaboration], “Evolution of the cosmic ray anisotropy above 10 eV,” Astrophys. J. 692, L130 (2009) [arXiv:0901.2740 [astro-ph.HE]].
-  T. Kifune, T. Hara, Y. Hatano, N. Hayashida, M. Honda, K. Kamata, M. Nagano and K. Nishijima et al., “Anisotropy of Arrival Direction of Extensive Air Showers Observed at Akeno,” J. Phys. G 12, 129 (1986).
-  M. G. Aartsen et al. [IceCube Collaboration], “Observation of Cosmic Ray Anisotropy with the IceTop Air Shower Array,” Astrophys. J. 765, 55 (2013) [arXiv:1210.5278 [astro-ph.HE]].
-  R. Abbasi et al. [IceCube Collaboration], “Observation of an Anisotropy in the Galactic Cosmic Ray arrival direction at 400 TeV with IceCube,” Astrophys. J. 746, 33 (2012) [arXiv:1109.1017 [hep-ex]].
-  K. G. Begeman, A. H. Broeils and R. H. Sanders, “Extended rotation curves of spiral galaxies: Dark haloes and modified dynamics,” Mon. Not. Roy. Astron. Soc. 249, 523 (1991).
-  A. Borione et al. [CASA-MIA collab.] “CASA-MIA Observation of Sun and Moon Shadows” Proc. 24th International Cosmic Ray Conference in Rome, Vol. 2, (1995) p.788.
-  A. Borione et al. [CASA-MIA collab.] “A Procedure to Obtain the Cosmic Ray Energy Spectrum from CASA-MIA Data” Proc. 23th International Cosmic Ray Conference in Calgary, Vol. 2, (1993) p.112.
-  R. Engel [Auger Collab.] “Test of hadronic interaction models with data from the Pierre Auger Observatory” Proc. 30th International Cosmic Ray Conference in Merida, Vol. 4, (2007) p.385.
-  http://www.hawc-observatory.org/observatory/sensi.php
-  John Pretz [HAWC Collaboration] “Limit on an Isotropic Diffuse Gamma-Ray Population with HAWC,” Proc. 34th International Cosmic Ray Conference in The Hague, (2015), PoS(ICRC2015)820, available at http://pos.sissa.it/archive/conferences/236/820/ICRC2015_820.pdf
-  Andrew J. Smith, “HAWC: Design, Operation, Reconstruction and Analysis” Proc. 34th International Cosmic Ray Conference in The Hague, (2015), PoS(ICRC2015)966, available at http://pos.sissa.it/archive/conferences/236/966/ICRC2015_966.pdf
-  http://wwwiexp.desy.de/groups/astroparticle/score/en/
-  Z. Cao [LHAASO Collaboration], “A future project at Tibet: The large high altitude air shower observatory (LHAASO),” Chin. Phys. C 34, 249 (2010). See also: http://english.ihep.cas.cn/ic/ip/LHAASO/
-  S. Vernetto and P. Lipari, “Exploring the gamma ray sky above 30 TeV with LHAASO,” Proc. 34th International Cosmic Ray Conference in The Hague, (2015), PoS(ICRC2015)740, http://pos.sissa.it/archive/conferences/236/740/ICRC2015_740.pdf
-  A. U. Abeysekara et al. [HAWC Collaboration], “Sensitivity of HAWC to high-mass dark matter annihilations,” Phys. Rev. D 90, no. 12, 122002 (2014) [arXiv:1405.1730 [astro-ph.HE]].
-  A. Cuoco, S. Hannestad, T. Haugbolle, G. Miele, P. D. Serpico and H. Tu, “The Signature of Large Scale Structures on the Very High Energy Gamma-Ray Sky,” JCAP 0704, 013 (2007) [astro-ph/0612559].
-  M. G. Aartsen et al. [IceCube Collaboration], “Search for Galactic PeV Gamma Rays with the IceCube Neutrino Observatory,” Phys. Rev. D 87, no. 6, 062002 (2013) [arXiv:1210.7992 [astro-ph.HE]].
-  M. G. Aartsen et al. [IceCube Collaboration], “Atmospheric and astrophysical neutrinos above 1 TeV interacting in IceCube,” Phys. Rev. D 91, no. 2, 022001 (2015) [arXiv:1410.1749 [astro-ph.HE]].
-  A. W. Strong, I. V. Moskalenko and O. Reimer, “Diffuse continuum gamma-rays from the galaxy,” Astrophys. J. 537, 763 (2000) [Astrophys. J. 541, 1109 (2000)] [astro-ph/9811296].
-  M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal and F. Sala et al., “PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection,” JCAP 1103, 051 (2011) [JCAP 1210, E01 (2012)] [arXiv:1012.4515 [hep-ph], arXiv:1012.4515 [hep-ph]].