TeV Dark Matter Searches in the Extragalactic Gamma-Ray Sky

High-energetic gamma rays from astrophysical targets constitute a unique probe for annihilation or decay of heavy particle dark matter (DM). After several decades, diverse null detections have resulted in strong constraints for DM particle masses up to the TeV scale. While the gamma-ray signature is expected to be universal from various targets, uncertainties of astrophysical origin strongly affect and weaken the limits. At the same time, spurious signals may originate from non-DM related processes. The many gamma-ray targets in the extragalactic sky being searched for DM play a crucial role to keep these uncertainties under control and to ultimately achieve an unambiguous DM detection. Lately, a large progress has been made in combined analyses of TeV DM candidates towards different targets by using data from various instruments and over a wide range of gamma-ray energies. These approaches not only resulted in an optimal exploitation of existing data and an improved sensitivity, but also helped to level out target- and instrument-related uncertainties. This review gathers all searches in the extragalactic sky performed so far with the space-borne Fermi-Large Area Telescope, the ground-based imaging atmospheric Cherenkov telescopes, and the High-Altitude Water Cherenkov Gamma-Ray Observatory (HAWC). We discuss the different target classes and provide a complete list of all analyses so far.


Introduction
For about a hundred years, powerful evidence has been accumulating that about 80% of the cosmic mass budget is yet of unknown nature, coined the dark matter (DM) [1]. Among the many proposed physical origins of DM, a promising hypothesis is the existence of weakly interacting massive particles (WIMPs) on the scale of GeV to TeV within some theory beyond the current Standard Model of particle physics [2]. However, despite the tremendous experimental efforts to test the WIMP hypothesis and various other DM candidates, the nature of DM is still elusive [3].
Searches in astrophysical gamma rays play a crucial role in these efforts. The use of astrophysical gamma rays to probe DM candidates was already proposed in the late 1970s in the GeV regime: possible signatures from annihilation of massive particles were first discussed for the extragalactic diffuse gamma-ray background [4] and the center of the Milky Way [5] to be searched in the data of early satellite gamma-ray detectors. Searches were followed up with the EGRET satellite [6] and lately perfected with unprecedented sensitivity by the Fermi-Large Area satellite Telescope (Fermi-LAT) [7].
Within supersymmetric extensions of the Standard Model, new particles on the TeV scale are favored as DM candidates [8], and gamma-ray searches close to the spectral endpoint at TeV energies can naturally provide an important puzzle piece for the identification of DM. At these energies, DM signals might be better distinguishable from decreased astrophysical backgrounds and may show prominent spectral features. However, such energies can hardly be probed by satellite detectors. Searching for DM particles in the TeV regime with ground-based Atmospheric Cherenkov telescopes was first discussed by [9] in 1992, the same year as the first extragalactic TeV gamma-ray source was detected by the Whipple telescope [10]. First    The background shows for illustration the Fermi-LAT counts map after 705 weeks (13.5 years) above 1 GeV. Note that data on galaxies and galaxy clusters were taken for multi-science purposes, while dSph targets were observed exclusively in hunt for DM. Targets written in italics are still unconfirmed dSph candidates.
Whipple 10-meter telescope and the HEGRA array were conducted a decade later towards nearby galaxies [11,12] and galaxy clusters [13]. DM searches with the current generation of imaging atmospheric Cherenkov telescopes (IACTs) were already detailed in the late 1990s for probing the Galactic center (GC) region [14] and external galaxies [15]. While the GC region promises a bright gamma-ray signal from DM annihilation, it suffers from background emission up to TeV energies [16], and major uncertainties about the expected DM signal were discussed [17]. These caveats resulted in an increased interest in alternative targets. To date, around 2500 h have been invested by IACTs in the search for DM on the TeV scale, with more than half of this time distributed on dozens of targets in the extragalactic sky. 1 In Figure 1, we illustrate the observations on all these targets outside the Galactic plane with the various IACTs, along with the respective exposures. Besides IACTs, also water Cherenkov detectors like the High-Altitude Water Cherenkov Observatory (HAWC), similar to Fermi-LAT constantly surveying a large fraction of the sky, are excellent instruments for TeV DM searches. In Figure 2, we show more than 100 DM targets in the extragalactic sky analyzed with HAWC and in the public data of Fermi-LAT.
This contribution reviews the current status of all these DM searches in the extragalactic sky, what we have learnt from them, and their particular challenges. In Section 2, we introduce the basic concepts of indirect DM searches. In Sections 3-6, we present the results of TeV DM searches obtained so far on the various targets. In Section 7, we draw a focus on multi-target, multi-instrument, and multi-messenger searches. Additionally, not only WIMP DM may show gamma-ray signatures-also axion-like particles (ALPs) and primordial black holes (PBHs) are viable DM candidates that can be searched for with gamma rays. While ALPs are covered in the contribution [19] to this Special Issue, we will discuss searches for PBHs in Section 8. Finally, we give a short outlook and conclude in Section 9.
Here, it is assumed that the DM is made of Majorana particles, so that DM particles are identical to their antiparticles. The signal from decaying DM reads as: In these equations, m DM denotes the mass of the WIMP, E γ the observed gamma-ray energies, dN/dE the number of photons per energy resulting from a single annihilation or decay process, and l.o.s the integration along the line of sight, the co-moving radial coordinate l. ∆Ω is the viewing cone. σv and t DM denote the velocity-averaged annihilation cross section and WIMP lifetime, respectively. For far-away gamma-ray sources, the flux is attenuated by absorption of photons on the extragalactic background light with the gamma-ray optical depth τ. While this attenuation is negligible for sources in the Local volume (z 10 −3 ) and sub-TeV energies, it may play a significant role for sources at higher redshift and/or if the WIMP mass is located at higher energies. In particular, for DM models on the PeV scale beyond the unitary bound [28], attenuation on the cosmic microwave background (CMB) can significantly suppress the signal on subgalactic scales [29].
The integrals over the density distribution in Equations (1) and (2) are called the astrophysical factor, usually denoted in short J-factor (or sometimes D-factor in case of decay). In case of annihilation, the signal scales with the number of DM particles times the DM density, as two particles from the same density distribution are involved in the annihilation process. For decay in turn, the signal scales merely with the number of particles in the viewing cone ∆Ω. Further, for annihilation, the flux scales with m −2 DM (see Figure 3, bottom panel), while only with m −1 DM for decaying DM. We will later see that this different scaling for decay allows for strong gamma-ray constraints on t DM for m DM larger than several TeV. The J-factor usually constitutes the largest uncertainty in indirect DM searches. This holds especially in case of annihilation, for which it depends on the density subclustering of a DM halo [30]. For the same reason, the factor (1 + z) 3 in front of the density integral in Equation (1) expresses that a higher density in proper coordinates at the time when the annihilations took place boosts the annihilation rate. 3 For most sources, their angular extension dJ/dΩ is expected to be strongly peaked towards the DM halo center and appears almost point-like to current-generation gamma-ray telescopes. The extension is larger in the case of decay, and being able to identify the extended emission profile may help to pinpoint the DM origin of a signal. In Figure 4, we show the total J-factors, J tot = J(∆Ω → 4π), 4 for both annihilation and decay, and for different source classes. It can be seen the larger uncertainties of J ann compared to J dec , especially when targets are expected to host a significant number of DM substructure (error bar caps marked by arrows). Note that comparing the merit of different targets is more complex, as some targets are affected by background gamma-ray emission or are very extended; the latter case poses a challenge in particular for instruments with a small field of view and residual isotropic backgrounds.   The vertical axis is scaled such that the number of photons in an energy interval can be read by eye when integrating over one order of magnitude (e.g., for annihilations into bb quarks, approximately 35 photons are created per annihilation between 1 and 10 GeV). Bottom panel: Different DM particle masses for a given annihilation channel. The vertical axis is scaled by m −2 DM , indicating that about 100 times less photons reach the observer when the DM mass increases by a factor 10. The sensitivity range of Fermi-LAT [44] (P8R3_SOURCE_V3 all events), MAGIC [45] (0-30 • after cuts), and HAWC [46] (with γ/h sep.) is illustrated by their effective areas (i.e., solid-angle-averaged acceptance for Fermi-LAT).
Additionally, the expected spectral energy distribution of gamma rays dN/dE is crucial for the DM interpretation of a signal. At the astrophysical source, it depends on the composition of particles produced after DM annihilation or decay. This branching composition, dN/dE = ∑ i f i (dN/dE) i into different channels i, is largely unknown and depends on the specific DM particle model. However, given a specific channel of annihilations or decays into Standard Model particles, the spectra can be robustly computed and are shown in Figure 3 (top). For interpreting null measurements, indirect DM searches usually scan individual annihilation or decay channels and derive limits under the assumption f i = 1 for each channel. The actual limits may indeed be roughly bracketed between the limit predicted by a soft gamma-ray spectrum (i.e., few gamma rays arriving in the instrument's energy range) and that by a hard channel. Assuming the DM to be cold, i.e., at rest at annihilation, all gamma-ray spectra share a sharp cut-off at the DM particle mass; momentum conservation guarantees that no gamma ray can be produced being more energetic than the rest mass of a single DM particle. In case of decay, this cut-off occurs at half the DM particle mass. After DM annihilation or decay into Standard-model matter, most of the high-energetic gamma-ray emission is produced by π 0 decays in hadronic cascades. Therefore, below the cut-off, most channels result in a rather featureless continuum spectrum (see Figure 3, top panel). Distinguished from these continuum spectra are cases where one or several gamma rays are directly produced in the annihilation or decay and carry away a large fraction of the DM particles' energy. This includes direct annihilation or decay into a gamma-ray pair, a gamma ray and another vector boson, or three-body processes including gamma rays like virtual internal bremsstrahlung. Such emission would not only provide compelling evidence for some exotic new physics, but also make it easier to detect for background-dominated detectors. However, such processes are for most DM models suppressed by orders of magnitude compared to annihilation or decay into pairs of heavy leptons or gauge bosons [41][42][43].   [31]. dIrr galaxies: values from [32] (light orange, 2021) and [33] (dark orange, 2018), distance values additionally from [34,35]. M31 (magenta): [36]. Galaxy clusters (green): from [37], with a 12% distance uncertainty according to [38]. Where error bar caps are shown as arrows, the uncertainty is assumed to be dominated by the uncertain signal boost from DM substructure. Bottom panel: same for decay according to Equation (2). dSph galaxies (purple): values from [31]. dIrr galaxies (orange): values from [39], with a 40% error according to [33]. M31 (magenta): [36]. Galaxy clusters (green): from [40].

Dwarf Spheroidal Galaxies
Dwarf spheroidal galaxies (dSphs) are low-luminosity satellite galaxies of the Milky Way or Andromeda (at least those detected to this date) and form the majority of the Local group galaxies [47]. Their name is due to their relatively small size, containing typically not more than 10 3 -10 7 stars [48], and to the fact that these stars follow more or less a spherical distribution. They are formed in the DM subhalos surrounding galaxies with enough baryonic matter to start and sustain stellar formation. With DM-dominated masses of about 10 6 -10 9 M [31,49], they are the DM targets with the largest mass-to-light (M/L) ratio detected to this date, although they could be surpassed by the predicted but yet still undetected dark DM subhalos (see Section 4). dSphs are old objects that only had a few generations of stars. As little to no star formation is occurring in these systems today, dSphs are considered as very clean targets for indirect DM searches and are expected to host no background gamma-ray emitters. dSphs typically orbit the Milky Way or the Andromeda galaxy at a few dozen to a few hundred kiloparsecs. In the Milky Way, this closely located to Earth, it allows to resolve the spectra of their stars individually and to estimate their DM profiles with a relative good precision in most of the cases, allowing relatively robust predictions of their J-factors. A relatively low amount of DM substructure is expected in dSphs, translating in an negligible boost of gamma-ray emission from DM annihilation [30].
The observed number of dSphs and their DM profile is not without posing questions though. The same simulations that predict the existence of such DM-dominated objects orbiting the Milky Way also predict cuspy profiles of the DM density profile near the centers [50,51]. In turn, spectroscopic measurements of dSphs show a variety of DM profiles from cores to cusps, a discrepancy known as the "core-cusp" problem [52]. Furthermore, simulations suggest the existence of more and heavier satellites than what has been observed so far, a problem known as the "missing satellites" problem [53,54]. While the problem is claimed to be largely understood [55][56][57], one explanation is that some dSphs lack baryonic mass (or lost some via unaccounted tidal effects) which could explain why their luminosity is too low to have yet observed them. However, this is in contradiction to what is expected for such massive objects which should host visible stars. This contradiction is known as the "too big to fail" problem [58]: these objects should be too big to fail to host star formation, at least in our current understanding of galaxy evolution and star formation.
dSphs were extensively used for the search of DM and adopted as dedicated targets for IACTs. Instruments with a large field of view such as Fermi-LAT and HAWC observed them as part of their continuous sky monitoring. The IACTs used their sensitivity to accumulate deep exposures on the promising targets shown in Figure [67], Triangulum II in 2020 by MAGIC [68], Willman 1, Draco, Ursa Minor, and Boötes I by VERITAS in 2010 [69], and a deep exposure with VERITAS on Segue 1 in 2012 [70]. At first, IACTs concentrated on these single targets and on deriving limits for single promising objects. At the same time, the lower sensitivity of Fermi-LAT and HAWC, with however steadily growing exposure and without the need of choosing a target to point at, let them focus on multi-target analyses from the beginning: the Fermi-LAT collaboration has published constraints on DM annihilation based on 11 months of data on 14 dSphs in 2010 [71], 2 years of 10 dSphs in 2011 [72], 4 years of 25 dSphs in 2014 [73], 6 years of 25 dSphs in 2015 [74], and 6 years of 45 objects in 2017 [75]. The latest limits by [76] using 11 years of the public Fermi-LAT data are shown in Figure 5 for annihilations into bb quarks. It can be seen that constraints are done with Fermi-LAT up to DM masses of 10 TeV. In 2018, HAWC published an analysis of 15 dSphs with 507 days of data [77]. The stacking approach pursued by Fermi-LAT and HAWC has the advantage of maximizing the sensitivity of searching a universal DM signal and was gradually also adopted by IACTs: H.E.S.S. published analyses of ∼140 h on 5 dSphs (Carina, Coma Berenices, Fornax, Sagittarius I, Sculptor) in 2014 [78] (continuum) and 2018 [79] (line emission) and of ∼82 h on Reticulum II, Tucana II, Tucana III, Tucana IV and Grus II in 2020 [80]. VERITAS published ∼230 h on 5 dSphs (Boötes I, Draco, Segue 1, Ursa Minor, Willman 1) in 2017 [81], and MAGIC ∼354 h on 4 dSphs (Coma Berenices, Draco, Segue 1, Ursa Major II) in 2022 [82]. In 2016, Fermi-LAT and MAGIC published a joint analysis of 15 dSphs observed by Fermi-LAT plus the ∼160 h accumulated by MAGIC on Segue 1 [83]. That analysis inaugurated multi-instrument gamma-ray searches for DM annihilation in dSphs (see Section 7). In Figure 5, we include all the latest combined searches by IACTs, compared against the Fermi-LAT and HAWC results. It can be seen that, while Fermi-LAT provides the strongest constraints in the GeV regime, competitive limits are obtained from IACTs and HAWC for DM masses above tens of TeV. In particular, the multi-instrument combination [83] still provides one of the strongest constraints for DM masses on the TeV scale. 5 While dSphs were so far mostly studied to search for DM annihilation, [84] have analyzed 6 years of Fermi-LAT towards 19 dSphs, providing strong constraints on the lifetime of GeV DM ( Figure 6). This illustrates the potential of dSphs also for searching for decaying DM.  [80]. Gray arrows indicate the projected sensitivity of future instruments at their best sensitivity in the m DM parameter space for 500 h observations of Segue I dSph with CTA [86] and 1 year combined dSph observations with LHAASO [87]. Analyses assumed different J-factors for the same targets. Thermal WIMP cross section prediction: analytic calculation from [88].
A complete list of the analyzed dSphs searched so far for DM annihilation or decay with Fermi-LAT and HAWC is given in Table A1, as completion of a similar table gathering all IACT observations in [18]. Furthermore, another recent review describing in depth the statistical analysis used in the search of DM in the direction of dSphs with Fermi-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS, can be found in [89]. Note that while a worsening of most of the dSph constraints by a factor of 2 to 7 had been recently claimed [90], a most recent re-evaluation suggests the validity of earlier J-factor estimations and constraints based on those [91]. Moreover, it was recently claimed for Sagittarius I that a gamma-ray emitting population of millisecond pulsars may be present in dSphs [92]. This would render the identification of a DM signal in these targets more challenging than long thought.

Dark Subhalos
Not all of the cosmic DM structures may have been massive enough to attract sufficient baryonic matter to form stellar systems. Numerical simulations predict a huge population of DM halos on scales lighter than those of dwarf galaxies, down to the order of Earth mass (10 −6 M ) and lighter [27]. Those halos are called "dark", as no baryonic interactions occur inside them which could reveal the positions of DM halos by electromagnetic radiation, except the emission from DM interactions itself. While the absolute rate of relic annihilations or decays in these light objects is lower, they may be found close to the Milky-Way disk and therefore even outshine dSph galaxies in their gamma-ray flux [93,94]. However, the existence of such low-mass DM halo population is hypothetical: attempts are done to identify those objects gravitationally on statistical and individual basis by lensing signatures [95] or stellar tidal streams [96][97][98][99][100][101]. 6 Moreover, their mass and spatial distribution in the Galaxy is uncertain, largely affecting the gamma-ray detection prospects at Earth [93,94,103].
Further information about the nature of unidentified gamma-ray sources can be obtained by follow-up observations with IACTs, with higher angular resolution and extending the energy range to search for TeV DM candidates. With respect to Fermi-LAT gamma-ray sources unassociated at that time, so far eight promising DM subhalo candidates have been subsequently observed by IACTs with the purpose of clarifying their possible DM origin: in 2010, the sources 1FGL J2347.3+0710 and 1FGL J0338.8+1313 have been observed by MAGIC [118], and between 2013 and 2015, the sources 2FGL J0545.6+6018 and 2FGL J1115.0− 0701 by VERITAS [119]. In 2018 and 2019, H.E.S.S. observed the four sources 3FHL J0929.2−4110, 3FHL J1915.2−1323, 3FHL J2030.2−5037, and 3FHL J2104.5+2117 [120]. None of the observations showed a signal in very high energy (E > 100 GeV) gamma rays. Due to no association with some gravitationally presumed DM budget, no individual constraints on the DM particle physics can be derived from null detections in such searches. However, reversely, under the assumption of some DM particle physics, constraints on the J-factor along the corresponding line of sight can be obtained [120]. Since their discovery, several of the GeV unassociated gamma-ray sources now have been classified to be probably not of DM origin, as summarized in Table 1. Notably, the nature of three of these sources is still elusive. In addition, follow-up observations of unassociated sources in Table 1. Current identification status of dark subhalo candidates so far observed by IACTs among Fermi-LAT previously unidentified gamma-ray sources. The likely non-DM association of 3FHL J2104.5+2117 was already noted by authors of [120].  [126,[129][130][131] the HAWC 2HWC catalog have been done by MAGIC [121] and VERITAS [122], though without significant detection and DM interpretation. A second approach is to head for serendipitous discoveries of dark subhalos in the field of view of IACTs. While no DM subhalo candidate has been detected with IACTs so far, projections have been made for discovery with the HAWC observatory [117] and the future Cherenkov Telescope Array (CTA) during large-scale surveys [93,132,133] or serendipitously in the field of view of other targets [133].

Source
Third and last, current non-detection or remaining dark subhalo candidates in gammaray observations can be translated into constraints on the DM annihilation cross section or lifetime, through assumptions on the global population of DM subhalos in the Galactic halo and their J-factors. Such limits on the annihilation cross section have been derived multiple times for the Fermi-LAT detector [104,110,111,113,115,116,[134][135][136][137][138] and also HAWC [117,139]. Analogous constraints from non-detection of DM subhalos in the gammaray data of IACTs are yet to be done.

Irregular and Spiral Local Galaxies
Indirect DM searches in individual galaxies other than Galactic dSphs have recently gained increased interest. In particular, this concerns local dwarf Irregular (dIrr) galaxies: these are dwarf disc galaxies rich in gas without a defined disk shape, and showing a relatively low rate of star formation that occurs outside of well-defined spiral arms [140]. In fact, this type of dwarf galaxies dominates the Local volume [141]: while most galaxies in the Local group are dSphs, around twice as many dIrr galaxies than dSphs are expected within a distance of ∼10 Mpc [33], and around 500 dIrrs have already been identified within that sphere [47,140]. The analysis of 36 spectroscopically well covered dIrrs of distances between 520 kpc and 16.6 Mpc revealed that their DM density profiles can be determined more robustly than for dSphs, and shows average masses of ∼10 10 M , larger than those of most studied dSphs [33,140]. Due to their large distances and a cored DM profile preferred by kinematic data, the authors of [33,140] have found annihilation J-factors for dIrrs generally 10 to 100 times fainter than the brightest dwarfs (see Figure 4, top panel). However, thanks to the large masses of these dIrrs, an annihilation boost from substructure may significantly increase the signal up to more than one order of magnitude in some cases [32] (Figure 4, top panel, light orange sample). Finally, the cored inner DM profile of dIrrs is still a matter of debate [32]. Notably, in contrast to dSphs, dIrrs are in principle expected to show an intrinsic gamma-ray background from star-forming regions. However, Ref. [33] has shown that this astrophysical background is negligible for most objects. 7 At GeV energies, a combined analysis for DM annihilation signals from the direction of 7 dIrr objects in the data from Fermi-LAT has been pursued [32]. Moreover, the Large Magellanic Cloud (LMC) [145], Small Magellanic Cloud (SMC) [146], and Smith High-Velocity Cloud [147] were studied. No DM signal from any of these targets was discovered, with constraints about a factor 10 to several hundreds weaker than those from stacked dSph galaxy data.
These searches have been complemented by recent searches in the TeV regime: H.E.S.S. has observed the Wolf-Lundmark-Melotte (WLM) dIrr galaxy [148] in the Southern hemisphere (see Figure 1) for 18 h, constraining an annihilation cross section in the order of magnitude σv 10 −21 cm 3 s −1 in the range of DM masses of a few TeV. A preliminary stacked analysis by HAWC [39] on both DM annihilation and decay in sources in the HAWC field of view (most of them clustered in the upper left of Figure 2) out of the sample from [140] observed for 1017 days provides annihilation constraints similar to those by H.E.S.S. on WLM, in the same energy range.
In addition, remote spiral galaxies like our own Galaxy, providing classical evidence for the existence of DM [149][150][151], serve as indirect DM targets. Spiral galaxies consist of a flat rotating disk with a central bulge and contain usually 10 10 to 10 12 stars, with DM halos of masses between 10 11 to 10 12 M , about a factor 5 to 10 heavier than the galaxies' stellar masses [152]. However, only the closest spiral galaxies are of interest for DM searches with gamma rays, among them most prominently the M31 (Andromeda) and M33 (Triangulum) galaxies. These two galaxies have been studied using Fermi-LAT by [153][154][155] for M31 and by [154] for M33. M33 (together with M32) was also searched for DM annihilation before by Whipple [11], and, as one of the very first DM searches with IACTs, the HEGRA data on M31 was analyzed for line signals from DM annihilation [12]. Lately, M31 was searched for DM at multi-TeV energies using HAWC [36]. Despite the very large J-factor associated with M31 (see Figure 4), DM searches in M31 are intricate. Extended gamma-ray emission from the inner 5 kpc (corresponding to 0.4 • in radius) of M31 [156,157] has already been detected by Fermi-LAT. Interestingly, similar to our own Galaxy, some excess above the expected astrophysical emission has been claimed [157], leading to the discussion of DM interpretations [158][159][160][161]. Moreover, like in the Milky Way, gamma radiation can occur in star-forming regions of the M31 disk or Fermi-bubble-like structures, visible under distances up to ∼1 • from the center of M31. In fact, a weak gamma-ray signal in the Fermi-LAT data has been claimed up to tens of degrees from the center of M31 [162], leading as well to DM interpretations [155,163].
All these caveats lead to gamma-ray constraints on DM annihilation on the TeV scale in M31, and similarly in M33, a factor 2 to 10 weaker than the best limits from dSph observations (also for HAWC, which does not see background gamma-ray emission). However, this picture changes when searching for DM decay: here, as shown in Figure 6, the best current limit on decaying DM with masses above 10 TeV is obtained by HAWC data on M31. This can be understood by the large M31 J dec -factor and the, compared to annihilation, less decreased flux from TeV DM decay, as described in Section 2.
Besides M31 and M33, many more spirals are present in the Local volume, although at much larger distances than M31 and M33. An analysis on four low-brightness galaxies (three of them classified as spirals, one as dIrr) resulted in limits about 4 orders of magnitude weaker than from stacked dSphs [164]. However, recently, Ref. [165] found that stacking a data set of the order of 10 5 low-brightness galaxies-which could be readily available with LSST-could provide constraints competitive with other targets. Despite the very early searches with the Whipple telescope, no spiral galaxies have been searched for DM with the current generation of IACTs.

Galaxy Clusters
Galaxy clusters form the largest overdensity structures in the Universe, with masses up to more than 10 15 M , out of which ∼80% is composed of DM [166]. The average DM density within the central hundreds of kiloparsecs of galaxy clusters is still tens of thousands times higher than the mean Universe's DM density. 8 This makes nearby cluster structures well suited DM targets.
The search for DM in galaxy clusters with gamma rays is generally hampered twofold: first, being located at cosmological distances, absorption at the extragalactic background light comes significantly into play (see Equations (1) and (2)). This attenuates the detectable radiation from these far-away objects in addition to the inverse square-distance law. Second, a variety of other physical processes generating high-energy gamma rays in clusters forms a non-negligible background for DM searches. While, so far, no diffuse gamma rays from cosmic-ray interactions in the inter-cluster medium of galaxy clusters have been clearly detected either by Fermi-LAT in the GeV [169][170][171] or by ground-based instruments in the TeV regime [172][173][174], many clusters host gamma-ray sources such as active galactic nuclei. As blazars, these may form a bright gamma-ray background, with even multiple blazars being present within a single cluster. For DM annihilation, additionally, the detection prospects in clusters strongly depend on the signal boost from the hierarchical DM substructure [30,175]. The J ann -factors (Equation (1)) of the average DM mass profiles of galaxy clusters are generally somewhat smaller than those of nearby dwarf galaxies. However, multiple levels of substructure are present in galaxy clusters, from the DM halos of their member galaxies, galaxy subhalos, and dark halos down to subhalo masses of 10 −6 M [176]. As illustrated in Figure 4 (top), this may boost the gamma-ray emission by  [180], stacked dwarf spheroidal galaxies [84], and stacked galaxy clusters [40] (from left to right). For DM masses larger than about 10 TeV, however, constraints are so far only obtained by ground-based instruments from observations of the Perseus cluster by MAGIC [178] or M31 by HAWC [36]. Gray arrows indicate the projected sensitivity reach of future instruments at their most sensitive point in the m DM parameter space for 300 h observations of the Perseus galaxy cluster with CTA [86] and 1 year of combined dSph observations with LHAASO [181].
factors of several tens, resulting in J ann -factors competitive with those of Galactic dSphs [37].
In Figure 5, we show the limits from H.E.S.S. observations of the Fornax cluster [85], where a boost of the emission by factor 100 over the considered integration area was assumed compared to the emission from a cluster halo without substructure. It can clearly be seen that such boost results in constraints similar to those from dSph observations. In the case of decay, a signal boost from substructure is not present. However, due to the large cluster masses seen under a small integration angle, nearby clusters can outplay Galactic dSph galaxies for DM decay searches [177], as illustrated in Figure 4 (bottom) by the J dec -factors.
Within an observation campaign covering multiple science purposes, the MAGIC telescopes have observed the Perseus galaxy cluster for almost 400 h [178]. As shown in Figure 6, using more than 200 h of good-quality data from this data set allowed to present one of the most stringent lower limits as of today on the lifetime of heavy particle DM with masses above 20 TeV. Based on less data, MAGIC [179] (also from Perseus), as well as Whipple [13] (from Perseus and Abell 2029), VERITAS 9 [173] (from the Coma cluster), and, as mentioned already above, H.E.S.S. [85] (from the Fornax cluster) derived limits also on DM annihilation from cluster observations. Many analyses have been performed so far with Fermi-LAT [40,169,[182][183][184][185][186][187][188][189][190][191][192]. Facilitated by Fermi-LAT's all-sky coverage, stacked analyses of many objects were performed, for continuum annihilation [40,186,191,192] and decay [40,183] spectra, and for line emission from annihilation [185,[188][189][190]192]. So far, HAWC has published a preliminary analysis for DM annihilation in the Virgo cluster [193]. While most stacked Fermi-LAT analyses focused on a selections of dozens or less most promising targets, others included hundreds [191] or thousands [194] of objects from galaxy catalogs. 10 Stacking of large target samples is also done in searches for cross-correlation between gamma rays and galaxy catalogs [195][196][197][198], gravitational tracers [199][200][201] or constrained simulations of the large-scale structure [202]. Additionally, without correlating to spatial emission structures, the diffuse gamma-ray background was searched for annihilation [203] or decay [40,204,205], as well as the diffuse background's angular power spectrum [206,207]. In Figure 2, we have included the 18 individual galaxy cluster targets studied multiple times for DM signals so far. Particularly for searches for DM decay, deep exposures on galaxy clusters by TeV instruments have turned out to be very powerful, as shown in Figure 6 (yellow curve) for the result of MAGIC towards the Perseus cluster of galaxies [178]. Such competitive limits compared to measurements at lower energies are possible because, as mentioned previously, the gamma-ray flux after DM decays does not drop as fast with the particle mass as for annihilation. Moreover, the angular emission profile from decay is much broader than the one from annihilation, and can be resolved, as in the case of Perseus, by the high angular resolution of IACTs. This allows to separate regions where only emission from DM is expected from conflicting astrophysical emitters, as e.g., the active galactic nucleus NCG 1275 in the Perseus cluster.

Multi-Wavelength, Multi-Messenger, and Multi-Instrument Combined Analyses
In recent years, there have been increasingly growing efforts in expanding the reach of DM searches with gamma rays. These expansions are following three main avenues: multiinstrument, multi-wavelength, and multi-messenger analyses. The first two already started with the first joint analysis by Fermi-LAT and MAGIC published in 2016 [83] and cover a wide range of masses from 10 GeV to 100 TeV. Since then, efforts have been expanded to a joint analysis by five gamma-ray detectors, namely, the gamma-ray satellite Fermi-LAT, the water Cherenkov observatory HAWC, and the IACTs H.E.S.S., MAGIC, and VERITAS. These efforts largely resulted in a standardization of the analyses developed in each collaboration and a consistent analysis framework. A preliminary combined analysis on DM masses from 5 GeV to 100 TeV could already exclude a thermal relic cross section of DM masses up to ∼200 GeV in the most optimistic scenario [208,209]. In the future, such joint searches can easily be extended to other gamma-ray detectors such as LHAASO which shows excellent capabilities to DM searches in the 1 TeV to 1 PeV range [87,181].
Beyond these experiments covering the GeV to PeV range, DM searches are also done using photons at lower energies, while using the same or similar targets, see e.g., [210][211][212][213][214][215][216][217][218][219]. For example, very strong constraints have been recently obtained on GeV to TeV DM annihilation from radio data on the LMC [220]. There is nothing preventing the effort of extending a combined and harmonized analysis to such wavelengths. Finally, the neutrino detectors ANTARES and IceCube recently reported on a combined analysis of their data searching for DM in the direction of the GC [221], and also a combination of DM searches between gamma rays and neutrinos would be promising. For further details on multi-messenger TeV DM searches, we refer the reader to the review by [222].

Primordial Black Hole Evaporation
At the end of Section 4, we have remarked that a blind search for gamma-ray signals from DM subhalos in the accumulated data sets of IACTs is yet to be done. However, archival "blind" searches in large IACT data sets have already been pursued for a DM candidate unrelated to WIMPs, namely, PBHs. This shows that gamma-ray data sets can be used to also constrain a completely different class of DM candidates. PBHs are thought to have formed in the radiation-dominated early Universe, and thus do not belong to the baryonic mass budget inferred by CMB measurements [223][224][225]. Therefore, they would be viable candidates to contribute to the cosmic DM budget. They can be searched with many different methods among which gravitational lensing, measurements of dynamical gravitational distortions like in stellar streams [226] or emission of gravitational waves [227].
In particular, lighter PBHs can be probed with the gamma rays from PBHs evaporating via Hawking radiation: any black holes lose mass, i.e., evaporate, via Hawking radiation. From the Hawking radiation into leptons, quarks, and gauge bosons, secondary gamma rays are produced traveling away from the evaporating black hole [228]. The closer a PBH comes to the end of its life, the stronger and energetic the gamma-ray emission becomes [228,229]. Specifically, the PBH lifetime is given by: where M is the mass of the PBH. Consequently, only PBHs created in the early Universe with masses heavier than ∼5 × 10 14 g = 2.5 × 10 −19 M have survived until today. This allowed for strong constraints to PBH masses slightly above masses of ∼5 × 10 14 g from measurements of the isotropic gamma-ray background [227]. At the very end of the PBHs life, the gamma-ray flux light curve of the evaporation is expected to increase following a power law, resulting in a sudden burst of energy in the last seconds. The gamma-ray spectrum of this final-stage emission extends up to TeV energies, and can be detected by gamma-ray instruments as sudden TeV gamma-ray bursts. However, there are two crucial differences to known gamma-ray bursts of cosmological origin: 11 First, unlike gamma-ray bursts from collapse or merger of compact extragalactic objects, detectable gamma-ray flashes from evaporating PBHs originate only from nearby PBHs of distance within several parsecs and therefore, are not attenuated by absorption from the extragalactic background light. Secondly, the duration of the peak PBH burst emission is very short, happening on timescales of seconds [228]. Unlike in cosmological gamma-ray bursts, no afterglow emission is expected from PBHs, and follow-up observations with pointed telescopes of PBH evaporation candidates seen by gamma-ray burst monitors are therefore virtually impossible. As a consequence, only the archival search for serendipitously recorded burst events in the data of IACTs and air shower arrays is possible. In addition, while measurements of the isotropic gamma-ray background have limited PBHs of masses 5 × 10 14 g to only contribute less than a fraction of 10 −6 to the cosmic DM budget, PBHs of this mass scale could have accumulated in galactic halos, resulting in a significantly higher burst rate detectable in their galactic neighborhood. While current limits correspond to a cosmic average PBH burst rate of 10 −6 bursts pc −3 yr −1 , in the case of galactic PBH clusters, these limits correspond to about one burst pc −3 yr −1 in the Galactic neighborhood according to [230]. Such a rate is in reach to be probed by gamma-ray detection of individual bursts. However, the search for individual burst events would be independent from isotropic measurements and from the assumption of a given clustering scale [228]. IACTs have been used to set constraints on the evaporation rate since Whipple [231] reached an upper limit of 1.08 × 10 6 pc −3 yr −1 (99% CL) on the evaporation rate in the Solar neighborhood (several parsecs) in the Galaxy, similar to the constraints found by the CYGNUS [232] and TIBET [233] air shower arrays. Since then, H.E.S.S. [234] (2013), [235] (2019) and [236] (2021) has improved these searches up to constraints of 10 3 pc −3 yr −1 for time intervals of the final evaporation between 10 and 100 seconds. VERITAS [237] (2012), [238] (2017) derived similar preliminary limits to a burst rate of ≤2.2 × 10 4 pc −3 yr −1 for a search window of 30 s using 747 h of data [239], with ongoing analyses to improve that result [239,240]. Limits have also been derived by Milagro [241], Fermi-LAT [242], and recently, HAWC [243]. A recent study for the Southern Wide-field Gamma-ray Observatory (SWGO) showed how the sensitivity to detect PBH evaporation bursts in TeV gamma rays can be pushed down to rates below 100 bursts pc −3 yr −1 [244] in the Solar neighborhood. These constraints show how individual PBH evaporation burst events can be competitively searched for with TeV gamma-ray instruments, and that small field of view telescopes like IACTs show capabilities comparable to large field of view detectors like Fermi-LAT and air shower arrays. Moreover, competitive constraints can be achieved with the MAGIC telescopes and the future CTA [245].

Future Searches and Conclusions
In this review, we have demonstrated the merits of the many targets in the extragalactic sky to search for TeV WIMP DM and PBHs with gamma rays. In Figure 5, the current best constraints on the WIMP annihilation cross section obtained from various targets and (combined) instruments are compiled for a smooth final-state spectrum, e.g., for annihilations into bb quarks). While for WIMP masses below ∼1 TeV, these constraints are mainly dominated by Fermi-LAT, for higher masses, results from IACTs and water Cherenkov detectors like HAWC become strongly competitive. Moreover, constraints from different targets, or after combining different targets, give comparable results: this illustrates the robustness of the constraints obtained with multiple targets. In the case of searches for decaying WIMPs, observations of extragalactic targets even compete with constraints from the Galactic halo, and for masses above several TeV, measurements with ground-based instruments compete with space-borne searches as illustrated in Figure  6. In addition, high-energy gamma-ray telescopes can contribute to searches for PBHs as potential alternative DM candidate by trawling their (archival) data for short-term gamma-ray bursts from evaporating PBHs in the Solar neighborhood.
In the recent years, the combination of the data from several instruments allowed to further improve the sensitivity for WIMP DM searches with gamma rays and to reduce the uncertainties from possible instrumental systematics. In the latest ongoing work by [208,209], a joint analysis over the data on 20 dSphs, obtained with five instruments using three different detection principles is performed. Such analyses may play an even increasingly important role in the future. Besides a variety of targets and different detection techniques of gamma rays, also different messengers like neutrinos can be included in a joint analysis. Such multi-instrument analysis approaches could be the key answer to the quest for heavy particle DM. This is already illustrated by strong constraints from radio signals [246], like e.g., towards the LMC [220] or Galactic dSph galaxies [247]. In any case, a diverse pool of observations on various targets and with various instruments will be crucial. This also includes dedicated searches for exotic spectral features like lines, "boxes", or sharp cut-offs, which may provide a distinctive smoking gun evidence for DM.
In parallel, next-generation instruments will provide a further boosted sensitivity to probe WIMP annihilations on the mass scale beyond several hundreds of GeV to multiples of TeV: big hopes lie here on CTA, which is currently being built. Consisting of a total of about a 100 IACTs divided into two arrays in the Northern and Southern hemisphere respectively, CTA is predicted to achieve a sensitivity to gamma-ray signals from TeV DM a factor 5 to 10 higher than current IACTs in the same observation time. This sensitivity is facilitated by (i) a decreased energy threshold, allowing to detect more of the low energy photons after WIMP annihilation or decay (see Figure 3), (ii) a factor 10 increased effective detection area, and (iii) even better angular resolution than current IACTs. In Figures 5 and 6, we show with the gray arrows how CTA can correspondingly reach a regime beyond the current limits towards dSphs. In Figure 5, one can see how the best sensitivity for CTA also shifts to lower energies compared to current IACTs, thanks to its improved low-energy threshold. Besides this, also the air shower array LHAASO [87,181] and the planned water Cherenkov detector SWGO [248] promise improved sensitivity to WIMP masses at beyond-TeV energies; in Figures 5 and 6, it is shown how LHAASO may provide unprecedented constraints from dSph targets beyond 10 TeV.
In conclusion, indirect searches for various DM candidates with high-energy gamma rays constitute a central building block in ongoing endeavor to unveil the nature of DM. In particular, indirect detection of TeV WIMP DM in gamma rays from multiple targets in the extragalactic sky will provide a crucial piece of information for a clear DM identification: the signal will be seen from the same objects that show gravitational evidence for DM, and with a universal signature from different targets, many of them free of confounding backgrounds. With this, we advocate these searches to be continued in the future, with novel detectors and improved analysis techniques. To this end, we hope that this review, and in particular Figures 1 and 2 and Table A1, may provide a clear guidance to which targets and searches are yet to be complemented with improved methods and observations. Acknowledgments: The authors would like to thank the Special Issue Editors John Quinn, Deirdre Horan, and Elisa Pueschel for the invitation to compile this review. Moreover, the authors would like to thank Javier Rico and Francesco G. Saturni for providing useful comments that helped to substantially improve the quality of the manuscript. This article makes use of matplotlib [249] and https://carto.com/carto-colors color schemes provided by https://jiffyclub.github.io/palettable (websites accessed on July 30, 2022).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: ALP Appendix A Table A1. Names and references to the DM analyses in Fermi-LAT and HAWC data of the sources displayed in Figure 2. For galaxy clusters, only sources are listed that are individually analyzed for DM in at least two publications in our literature search. Further clusters are, e.g., studied in [186,190,191] and implicitly in cross-correlation or power spectrum studies (see Section 6 for further references.) Figure 2 Source Name Reference  Figure 2 Source Name Reference  Centaurus (A 3526) [182,183,185,186,188,189,191,192] 103 Coma cluster [40,[182][183][184][185][186][188][189][190]192 Ophiuchus [188-190] 111 Perseus [188-190,192] 112 Virgo cluster [169,184,185,[187][188][189][190][191][192][193] Notes 1 In this review, we consider all targets outside the Galactic plane. However, we omit DM searches in high-latitude globular clusters. A recent review including DM searches in the inner Galaxy, Galactic globular clusters, and intermediate black holes, as well as their particular challenges, can be found in [18]. 2 Gamma rays can additionally be produced in the astrophysical surrounding by secondary processes like inverse Compton scattering of high-energetic e + /e − produced after DM annihilation or decay [25,26]. 3 This factor is usually neglected due to the large uncertainties in the density distribution entering the J-factor. 4 The size of a DM halo is usually limited by a defined virial or characteristic overdensity radius. Accordingly, the integral of Equations (1) and (2) vanishes at angular distances from the halo center larger than this radius. 5 Note, however, that these compared analyses do not use consistent J-factor values for the same targets. 6 Additionally, a reverse approach was undertaken to search for faint dwarf galaxies with optical telescopes towards the sky directions of Fermi-LAT unidentified gamma-ray sources [102]. 7 Exceptions are the Large and Small Magellanic Clouds, from which gamma-ray emission has already been detected [142,143]. Note that the classification of Magellanic spirals like the Large Magellanic Cloud as dIrr galaxies is ambiguous [144]. 8 This corresponds to the densities within the so-called characteristic scale radius, r s . Usually, the total extension of a DM halo is given by R 200 , defined as the radius within which the mean density is 200 times higher than the critical density, or about 600 times higher than the mean DM density, ρ m . For common clusters, the value of R 200 is several megaparsecs. For comparison, the average DM density in the inner r s ≈ 20 kpc of the Milky Way is about (1 ± 0.2) × 10 5 ρ m [167,168]. 9 While VERITAS [173] published a joint analysis of Coma together with Fermi-LAT data, their DM limits are only given for the VERITAS observations alone. 20 of 31 10 Note that in these cases of very large samples, sky regions of cluster structures and other closely located objects overlap, as e.g., individual local dIrr galaxies in [191]. 11 See [228] for a complete list of differences between cosmological and PBH gamma-ray bursts.