First implications of Tibet AS$_\gamma$ data for heavy dark matter

Extensive air shower detectors of gamma rays in the sub-PeV energy region provide a new and relatively unexplored window for dark matter searches. Here we derive some implications of the recently published Tibet AS$_\gamma$ data for decaying dark matter candidates. The available spectral information is already useful in obtaining competitive constraints, surpassing existing limits above 10 PeV mass for hadronic or massive boson final states. This is particularly true if accounting for a benchmark astrophysical background of Galactic cosmic rays in the (0.1-1) PeV range. By relying on the arrival distribution of the photons, we show that significantly better sensitivity can be attained, comparable or better than IceCube also for most leptonic final states. Full data exploitation requires however further information disclosure.


I. INTRODUCTION
The last decade has seen the opening of two windows in astrophysics: The high-energy neutrino one [1,2] and the gravitational wave one [3]. As it is usually the case, these seminal discoveries have defied some of the expectations and have stimulated new fascinating questions. For instance, the origin of the bulk of IceCube events is still a puzzle, with mounting indications that they could involve new classes of objects, either astrophysical or exotic ones (for recent studies in that sense, see [4,5]). Decaying dark matter (DM) has long been recognized as one of these possibilities [6,7], capable at the same time to account for the peculiar energy and angular distributions of the events [8] 1 . In [11], we made the point that extensive air shower (EAS) detectors sensitive to sub-PeV gamma-rays could provide amazing sensitivity to this class of DM models, either supporting such an interpretation of IceCube data or significantly improving the constraints. We forecasted that the forthcoming generation of observatories would attain the needed sensitivity for crucial studies. With the first detection of Galactic gamma-rays up to PeV energies by Tibet AS γ [25] and the commissioning of the LHAASO observatory [26], we are witnessing the dawn of this new astrophysical window, and should prepare to harness its full potential, in particular in a multiwavelength and multimessenger context (see [27][28][29] for initial examples of this approach). In this article, we provide an assessment of the importance of this first detection for DM searches, a task whose relevance has already been recognized in [30]. This energy range has its own peculiarities, notably the absence of any extragalactic contribution (since the gamma-rays are fully absorbed on the CMB and IR background, cascading down to sub-TeV energies) and the need to take into account partial and anisotropic absorption of the gamma-rays even within the Galaxy, as detailed in [11]. In Sec. II, we briefly introduce the dataset and the model used, deriving spectral constraints. Sec. III illustrates the power of the angular analysis, whose full exploitation will have to wait for the disclosure of further experimental information. In Sec. IV, we report our conclusions. For the first time in this context, we also include a model of Galactic diffuse gamma-rays according to [31] to assess the impact of astrophysical backgrounds on DM searches in this window. Due to the very rudimentary understanding of the CR sources and propagation properties in this regime, this has to be seen as a motivation for a better assessment of these predictions and of their uncertainties.

II. SPECTRAL ANALYSIS
The Tibet air shower and muon detector array has detected for the first time a sub-PeV gamma-ray flux from the Galaxy [25]. The collaboration reports energy spectra in two angular windows in Galactic coordinates: Region I ( |b| < 5 • and 25 • < l < 100 • ) overlaps with the region from which the ARGO-YBJ collaboration reported a diffuse detection in the TeV range [32]; Region II (|b| < 5 • and 50 • < l < 200 • ) matches the one where CASA-MIA previously only reported upper limits [33]. DM bounds derived from both regions are comparable, with Region II marginally better (at the 10% level) in * Electronic address: arman@puc-rio.br † Electronic address: serpico@lapth.cnrs.fr 1 For further studies of heavy DM signatures in neutrino/gamma-ray experiments see [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].  . The blue and black solid curves show the tentative bounds from angular and spectral analyses (without background), respectively (see text for methodological explanations); the red curves show the 95% C.L. bounds from the spectral analysis including the Galactic contribution from [31]. The brown dashed and green dot-dashed curves show the 95% C.L. bounds from IceCube [34] and Fermi-LAT [35], respectively. some regions of parameter space. Bounds from ARGO-YBJ data in Region I kick in and dominate at sub-PeV masses. Below, for each DM decay channel, the best limit from Regions I and II and of either experiment is reported at each value of DM mass. For the DM signal template, we consider the DM decay models discussed in [34], which we address the reader to for further details. For the spectra of gamma-rays and e ± from DM decay we use the results of [37]. Anisotropic absorption of the gamma-rays in Galaxy via e ± production on IR/optical/UV photons (plus the isotropic absorption onto CMB) is accounted for as in [11]. Also, the inverse-Compton scattering of e ± from DM decay off the ambient photons (mostly CMB photons) has been taken into account (for technical aspects, see [11]), although this contribution remains < ∼ a few percent not only in the Galactic plane region, but also for |b| < 40 • , which is considered in Sec. III. For the DM halo we use the Navarro-Frenk-White density profile [38] with the critical radius 24 kpc and Sun-GC distance of 8.3 kpc, yielding the DM density at Solar System 0.39 GeV cm −3 (for further details see section 3 of [11]).
The bounds derived from the spectral analysis of Tibet AS γ data, reported in three energy bins 100 < E γ /TeV < 158, 158 < E γ /TeV < 398 and 398 < E γ /TeV < 1000, are presented in Fig. 1 for twelve decay channels. The black solid curves in Figure 1 show the bounds derived simply by requiring the γ-ray flux from DM to remain below the 1σ upper limits reported by Tibet AS γ in all the three energy bins. The upturns in sub-PeV DM masses are the consequence of ARGO-YBJ data at ∼ 1 TeV. Inclusion of the "space-independent" astrophysical model of diffuse gamma-rays from ref. [31], leads to stronger exclusion. The red solid curves in Fig. 1 show the 95% C.L. bounds obtained from χ 2 analysis that includes this Galactic contribution. For DM masses m DM < ∼ 1 PeV, the bounds from spectral analysis are weaker than the limits obtained from diffuse gamma-ray data of Fermi-LAT [35] for almost all the decay channels. For some channels, like the u and b-quark, the spectral analysis bounds improve over the IceCube's limits [34] at m DM < ∼ 1 PeV. Note that these bounds are based on spectra extracted from Pythia's results (for details, see [36]), while ours rely on [37]: We checked that differences for masses below O(10) PeV are typically at the few percent level, hence irrelevant for our purposes.
The newly derived spectral limits surpass Fermi-LAT background limits for masses above ∼ 10 PeV for all hadronic final state channels, as well as massive gauge-bosons and Higgs. While for light leptons they do not improve over Fermi, they do for third generation leptons already at the PeV scale, although they are not yet competitive with IceCube. It is worth noting that, for light or b-quark final states, Tibet AS γ 's spectral bounds surpass both Fermi and IceCube bounds above ∼ 10 PeV, in particular if the astrophysical background is accounted for.

III. ANGULAR DISTRIBUTION ANALYSIS
Interestingly, the measured diffuse gamma-ray flux by Tibet AS γ and reported in [25] is compatible with an excess with respect to the existing evaluation of Galactic diffuse gamma-ray flux of [31], in particular in Region II which is farther away from the inner Galaxy. This is qualitatively in agreement with expectations from decaying DM models, whose flux declines very gently at directions away from the Galactic Center. Merely based on spectral information, an excess could admit a DM decay interpretation, as illustrated in the left panel of Fig. 2: The black dashed curve shows the gamma-ray flux from DM → bb with the mass and lifetime (m DM , τ DM ) = (30 PeV, 6 × 10 27 s) which is at the edge of exclusion (see the red solid curve in the middle panel at the top row of Fig. 1), and is compatible with the limits from IceCube and Fermi-LAT. The blue solid curve, showing the sum of DM and Galactic contributions, provides a better fit to the Tibet AS γ data than solely the Galactic contribution. The DM model in the left panel of Fig. 2 is only one example among various spectrally viable DM models, particularly if one allows DM decay into multiple channels. Angular information of Tibet AS γ 's data provides an exquisite handle to probe a wider range of DM models and test the viability of those that may even appear favoured by spectral considerations only and is exemplified in the left panel of Fig. 2. Unfortunately, a complete angular analysis is not possible at the moment: Although the collaboration reports in [25] the Galactic latitude profiles of data in the three energy bins, the use of these data is made difficult by the lack of precise information on the actual experimental aperture (or exact exposure) and γ-efficiency, especially its energy dependence. Also, the authors of [25] renormalize (separately for each energy bin) the Galactic model of [31] to the observed number of events within |b| < 5 • , which prevents inferring the detector's characteristics. In addition, the background estimation at high Galactic latitudes needs clarification and likely improvement. The method currently employed in [25] for the first two energy bins (100 < E γ /TeV < 158 and 158 < E γ /TeV < 398) yields 1σ upper limits on the number of events (or excess of events, subtracting the estimated background from observed events) which are negative already at |b| ∼ 10 • . Taking these upper limits at face value, no model (either conventional Galactic models or DM decay ones) would survive. As a consequence of these limitations, here we just illustrate how powerful such an analysis could be, by using the Galactic latitude distribution of observed events in the third energy bin 398 < E γ /TeV < 1000. The right panel of Fig. 2 shows the latitude profile of events in the highest energy bin. The black curve is the expectation of the DM model considered in the left panel and the red points show the Tibet AS γ data with 1σ error bars. In the calculation of the latitude profile for DM we assume that the acceptance efficiency of the detector depends on declination (and is uniform with respect to right ascension) according to the geometrical arguments of [39]. The blue shaded histogram shows the expected number of events from the Galactic model of [31] rescaled back to its true normalization by a factor equal to the ratio of the averaged Galactic flux to the measured flux. Clearly, the angular distribution of events at |b| > 5 • in the highest energy bin, without considering the Galactic model, excludes the DM model of the left panel in Fig. 2.
Following the crude application of not-overshooting the 1σ upper limits on angular profile of observed events in the highest energy bin (without taking into account the Galactic contribution), the blue curves in Fig. 1 represent our tentative bounds from angular considerations. For hadronic final states, as well as massive gauge and Higgs bosons, Tibet AS γ data constitute the most sensitive probe of decaying DM above a few PeV. Even for leptonic final states, they become comparable or better than IceCube at high DM masses, and noticeably better than Fermi-LAT limits. Once the angular distribution of events in the lower energy bins could be properly analyzed, we anticipate a significant improvement in sensitivity also at lower DM masses. Proper availability of these data may be also useful to probe other large-scale astrophysical Galactic features [40].

IV. CONCLUSIONS
In this article we have presented a first examination of the impact of the recent detection of Galactic gamma-rays up to PeV energies by Tibet AS γ [25] for heavy decaying DM models. Our main conclusions are that the spectral data alone improve over previous EAS bounds, making them competitive (i.e. sometimes better than best existing ones) above 10 PeV or so, particularly for hadronic or massive boson final states. This is particularly true if the astrophysical Galactic background is included. The angular distribution of the data holds an even greater potential. We illustrated how the angular analysis would make Tibet AS γ competitive or better than IceCube even for leptonic final states which generally produce more neutrinos, and undoubtedly better than Fermi-LAT bounds.
Our two main recommendations are: i) To the Tibet AS γ collaboration, to disclose more information in order for the sensitivity power of the angular analysis to be fully exploited. In particular, detailed information of the detector exposure and its latitude dependence (in case of difference with respect to the geometrical corrections) and the energy-dependence of γ-efficiency are essential. Also, an improved estimation of background events at high latitudes is crucial for the viability of angular analysis in lower energy bins. ii) To the cosmic ray theory community, to refine our understanding of the astrophysical background at PeV energies and of its uncertainties, which may prove instrumental also to tighten bounds by typically 30% or so, according to the current benchmark.