Gamma-Ray Dark Matter Searches in Milky Way Satellites -- A Comparative Review of Data Analysis Methods and Current Results

If dark matter is composed of weakly interacting particles with mass in the GeV-TeV range, their annihilation or decay may produce gamma rays that could be detected by gamma-ray telescopes. Observations of dwarf spheroidal satellite galaxies of the Milky Way (dSphs) benefit from the relatively accurate predictions of dSph dark matter content to produce robust constraints to the dark matter properties. The sensitivity of these observations for the search for dark matter signals can be optimized thanks to the use of advanced statistical techniques able to exploit the spectral and morphological peculiarities of the expected signal. In this paper, I review the status of the dark matter searches from observations of dSphs with the current generation of gamma-ray telescopes: Fermi-LAT, H.E.S.S, MAGIC, VERITAS and HAWC. I will describe in detail the general statistical analysis framework used by these instruments, putting in context the most recent experimental results and pointing out the most relevant differences among the different particular implementations. This~will facilitate the comparison of the current and future results, as well as their eventual integration in a multi-instrument and multi-target dark matter search.


Introduction
The existence of a dominant non-baryonic, neutral, cold matter component in the Universe, called dark matter, has been postulated in order to explain the kinematics of galaxies in galaxy clusters [1] and stars in spiral galaxies [2], as well as the power spectrum of temperature anisotropies of the cosmic microwave background [3]. In one of the most plausible and thoroughly studied theoretical scenarios, dark matter is composed of weakly interacting particles with mass in the range between tens of GeV and hundreds of TeV, generically referred to as WIMPs [4]. The Standard Model (SM) particles that could result from WIMP annihilation or decay would hadronize, radiate and/or decay, producing detectable stable particles such as photons, neutrinos, proton-antiproton pairs or electron-positron pairs [5]. Looking for unambiguous spectral and/or morphological signatures of dark matter annihilation or decay in the extra-terrestrial fluxes of those particles is usually referred to as indirect dark matter searches.
Gamma rays are promising messengers to search for WIMPs. Since they are electrically neutral, they are not deflected by magnetic fields and point back to their production site, and therefore could be used to determine the underlying dark matter spatial distribution. At non-cosmological scales, gamma rays are also essentially unaffected by energy losses, which would preserve the features expected for dark matter annihilation and/or decay spectra, which depend on the values of the dark matter mass and the branching ratios to the different annihilation/decay channels, which could thus be studied. Finally, the 2 of 27 gamma-ray signal intensity would depend on the annihilation cross-section or the decay lifetime, which could therefore be determined if we measured a signal from an astronomical site for which we have a good estimate of its dark matter content based on independent measurements and/or simulations.
N-body simulations predict the formation of cold dark matter haloes in a hierarchical clustering fashion [6]. dSphs form in dark matter galactic subhalos that contain enough baryonic matter to have activated stellar formation (pure dark matter halos should also exist, but they remain as of yet unidentified). They are irregular satellite galaxies with mass ∼ 10 7 M and the largest known ratios of dark to luminous matter. The extension of the expected gamma-ray emission from the Milky Way dSphs is typically between ∼ 0.1-0.5 • [7], which is of the order of the angular resolution of most of the current-generation gamma-ray telescopes.
Gamma-ray telescopes of the current generation have performed extensive observational campaigns of dSphs in search for dark matter signals. Along the years, gamma-ray telescopes have progressively adopted state-of-the-art statistical analysis techniques for their dark matter searches, optimized to exploit the particular spectral and morphological features expected for the signal. All the instruments have converged into a general statistical analysis framework, albeit with some significant differences among the different implementations. Some of these differences are unavoidable, since they are needed to adapt the analysis to the different experimental scenarios, whereas others rather consist in choices of conventions, approximations, or simplifications. These latter ones include the methods for computing the spectral and morphological models for the expected gamma-ray signal and associated background, their use in the statistical analysis, and the treatment of the related statistical and systematic uncertainties. Understanding both the similarities and the differences among the various analysis implementations is fundamental in view of meaningful comparison and combination of the obtained results.
In this paper, I review the present status of indirect dark matter searches with observations of dSphs with gamma-ray telescopes. In Section 2 I summarize the formalism for the computation of the gamma-ray fluxes expected to be produced by dark matter processes in dSphs. In Section 3, I briefly introduce the current generation of gamma-ray telescopes, their working principles and main features. Section 4 is devoted to the detailed description of the common statistical data analysis framework used by all these instruments in their search for dark matter in dSphs. Finally, in Section 5, I perform a critical comparison of the particular analysis implementations, review and contextualize the latest experimental results published by the different instruments, and show the prospects for their near-future combination.

Gamma-Ray Signals From Dark Matter Processes in dSphs
dSphs are among the cleanest astronomical targets for indirect dark matter searches. They are thought to be highly dominated by dark matter (mass-to-light ratios of the order of 10 3 [8]), and they harbor no known astrophysical gamma-ray sources that could produce a relevant background. Furthermore, dSphs contain in general no significant amount of dark gas, which allows their dark matter distribution to be inferred with relatively good precision from the stellar motions, enabling in turn robust predictions of the intensity of the associated gamma-ray signals, generally within an accuracy of one order of magnitude [7]. Finally, given how most of the known dSphs sit on relatively clean interstellar environments (i.e., out of the Galactic plane, where the particle densities, cosmic ray fluxes and radiation fields are small), the expected gamma-ray signal would come from well-understood prompt processes. Secondary processes such as inverse Compton scattering of primary or secondary electrons, or gamma-ray cascading processes initiated by their interaction with radiation fields (hence depending on local details of those radiation fields), can be in general ignored when computing the gamma-ray flux expected from dark matter at dSphs. Therefore, since flux predictions rely on relatively few assumptions compared to other typical observational targets 3 of 27 like e.g., the Galactic center or clusters of galaxies, the bounds on the WIMP properties that can be inferred from the presence or absence of a gamma-ray signal are also relatively robust.
If WIMPs (hereafter denoted by χ) concentrate with number density n χ in a dSph, annihilating and/or decaying with a rate Γ χ and an average isotropic gamma-ray spectrum dN γ dE , then the differential flux of gamma rays of energy E observable from Earth coming from directionp, per unit energy and solid angle Ω, is given by the following expression: with l the distance from Earth and the corresponding integral running over the line of sight in the direction p.
As explicitly noted in Equation (1), dE contains all the spectral dependence of the gamma-ray flux, and therefore determines the probability density function (PDF) for the energy of the emitted gamma rays. On the other hand, all the morphological dependence is contained in the line-of-sight integral, which hence determines the PDF for the gamma-ray arrival direction. Given that we can make relatively reliable predictions about these two PDFs, they will constitute key ingredients in the maximum-likelihood data analysis, as we will see below in detail.
The expected primary products of the WIMP annihilation and decay processes are pairs of leptons, quarks or gauge bosons, which would produce secondary gamma-rays (among other stable products) through final-state radiation or hadronization+decay chains. It is straightforward to compute the contribution to dN γ dE from the different annihilation/decay channels, for a given WIMP mass, using standard Monte Carlo simulation packages such as PYTHIA [9]. The spectral energy distribution of the gamma-ray continuum resulting from these processes peaks between one and two orders of magnitude below the WIMP mass, depending on the channel, as shown in Figure 1. The plots show that Fermi-LAT is the most sensitive instrument for searching for WIMPs up to a dark matter mass (m χ ) of few TeV in the case of bb channel and of few 100 GeV for the τ + τ − channel. Cherenkov telescopes dominate the search between those masses and ∼ 100 TeV for bb and few 10 TeV for τ + τ − , and HAWC for even higher WIMP masses. Primary gamma rays like, e.g., those from the χ[χ] → γγ or χ[χ] → γZ processes would be [quasi-]monochromatic. These would constitute the cleanest possible dark matter signal, given how there is no known astrophysical process able to produce such gamma-ray spectral lines, and that backgrounds affecting the measurement could be drastically reduced using spectral criteria. If detected, a gamma-ray line would by itself be considered a clear evidence for the presence of dark matter. However, due to parity conservation, primary gamma rays can only be produced via loop processes, which significantly reduces their associated rate Γ χ .
It is useful to particularize the line of sight integral in Equation (1) for the annihilation and decay cases: • For annihilation, Γ χ = 1 k n χ σv , with σv the average of the product of the WIMP velocity and annihilation cross section. The value of k depends on whether WIMPs are Majorana (k = 2, to take into account that an annihilation involves two identical particles) or Dirac particles (k = 4, reflecting the fact that particles can only annihilate with their-equally abundant-antiparticles). Including this into Equation (1), and writing the WIMP number density n χ in terms of its mass and density (ρ), we obtain:  Figure 1. Expected gamma-ray spectral energy distribution for WIMPs of masses m χ = 0.01, 0.1, 1, 10 and 100 TeV annihilating with σv = 3 × 10 −26 cm 3 s −1 into bb (left) and τ + τ − (right) pairs in a dSph with associated J-factor J ann = 5 × 10 21 GeV 2 cm −5 ; also shown are the sensitivity curves for the instruments considered in this paper. Fermi-LAT sensitivity curve [10] corresponds to observations of a point-like source at Galactic coordinates (l, b) = (120 • , 45 • ) for 10 years, analyzed using the latest (Pass8) data reconstruction tools; HESS [11], MAGIC [12] and VERITAS [13] curves correspond to 50 h of observations of a point-like source at low (Zd 30 • ) zenith distance; HAWC curve [14] is for five years of observations of a point-like source at a declination of +22 • N. The flux sensitivity for 50 h observations with the future Cherenkov Telescope Array [15] is shown for comparison.
where we have defined the annihilation differential J-factor as: • For decay, the rate is given simply by the inverse of the dark matter decay lifetime, i.e., Γ χ = τ −1 χ , since each WIMP particle decays independently of each other. Including this into Equation (1), we get: where we have defined the decay differential J-factor as: The J-factor in a region of the sky ∆Ω is given by: both for J ann and J dec . It is convenient to define the total J-factor for a given dSph as: with ∆Ω tot a region of the sky containing the whole dSph dark matter halo. The differential J-factor can be written as: where dJ dΩ can be interpreted as the PDF for the arrival direction of gamma rays produced by dark matter processes in the dSph halo, since ∆Ω tot dΩ dJ dΩ = 1. Using this notation, the differential gamma-ray flux per energy and solid angle can be written as: (with a being either a ann ≡ 1 4π σv k m 2 χ for annihilation or a dec ≡ 1 4π 1 τ χ m χ for decay). The differential flux per unit energy is given by: The distribution of dark matter within the halo, ρ(p, l), is usually estimated by solving the spherical Jeans equation for the stellar kinematic data [16]. Using this technique, several authors have produced catalogues of J-factors for the different known dSphs. In general, the classical dSphs, with relatively large stellar populations (O(100 − 1000)), have relatively low associated J-factors (typically between 3 × 10 17 and 7 × 10 18 GeV 2 cm −5 within an integrating angle of 0.5 • ), with associated uncertainties also relatively low (typically below 50%), suitable for setting robust limits to dark matter properties. On the other hand, members of the ultra-faint population (those discovered by the Sloan Digital Sky Survey or later, with O(10 − 100) members stellar populations) can have larger estimated J-factors (some above 10 19 GeV 2 cm −5 ) but also larger uncertainties (some above a factor 10), therefore providing better prospects for discovery but less robust constraining power. A detailed review about the expected dark matter content and distribution of the known dSphs can be found elsewhere in this volume.

Gamma-Ray Telescopes
For WIMP indirect searches with gamma-rays, the relevant energy range spans from 100 MeV to 100 TeV (see Figure 1). Photons of these energies interact in the upper layers of the atmosphere, making impossible their direct detection from the ground. Several different experimental techniques have been developed to detect gamma rays, each optimized for a different energy range and hence for different dark matter masses.
At energies below ∼100 GeV, we can efficiently measure gamma rays before their destructive interaction in the atmosphere by direct detection with balloon or satellite-borne detectors. Gamma rays interact within the detector, and convert into e + e − pairs, which are tracked to estimate the direction of the primary particle, and then stopped by a calorimeter to estimate its energy. This method is limited by the relatively small achievable collection area, corresponding essentially to the physical size of the detector. On the other hand, the technique presents the great advantages of ∼100% duty cycle, large field of view, and that the much more abundant charged cosmic rays can be easily identified and therefore vetoed, resulting in virtually background-free gamma-ray measurements. Currently, the most advanced gamma-ray telescope using this detection technique is the Fermi-LAT. It consists of a large-field-of-view (2.4 sr), pair-conversion telescope, sensitive to gamma rays in the energy range between 20 MeV and about 300 GeV [17]. The latest Fermi-LAT source catalogue contains about 5000 sources [18], a third of which remain unassociated. Since its launch in June 2008, the LAT has primarily operated in survey mode, scanning the whole sky every 3 h. The exposure coverage of this observation mode is fairly uniform, with variations below 30% with respect to the average exposure. Thanks to this full-sky coverage, Fermi-LAT will be able to perform dark matter searches using its data archive should new dSphs be discovered in the future. 6 of 27 Above few tens of GeV, gamma-ray fluxes become too low for the relatively small collection area of Fermi-LAT, and it is advantageous to measure them indirectly through the detection of the secondary particles and/or the radiation present in the particle cascade resulting from their interaction in the atmosphere, which greatly increases the effective collection area.
Cherenkov telescopes measure the Cherenkov radiation emitted by the electrons and positrons of the cascade (which travel faster than light in the atmosphere), thus producing an image of such cascade. The intensity, orientation, and shape of Cherenkov images allow for the estimation of the energy and arrival direction of the primary particle, and provide some separation power between gamma rays and charged cosmic rays. Several nearby telescopes observing the same gamma-ray source may image the same cascade from different perspectives, increasing the precision of these measurements. The weak points of this technique are the small duty cycle (about 10-15%, since they operate only during night, with no or relatively dim moonlight and good atmospheric conditions), narrow fields of view of few degrees diameter at most, and the presence of the irreducible background produced by charged cosmic rays. Among its advantages, we find the large collection area, given by the size of the Cherenkov light pool projected on the plane of the telescope reflector (e.g., ∼10 5 m 2 for 1 TeV gamma ray at low zenith distance). The Finally, water Cherenkov particle detectors measure the charged particles present in the cascades initiated by the primary gamma rays when interacting in the atmosphere. The amount of detected particles and their spatial distribution allow to measure the energy of the primary and to discriminate between gamma rays and cosmic rays, whereas the difference of detection time at different detectors allows to estimate the arrival direction. This technique is sensitive to gamma rays and cosmic rays between few hundred GeV and 100 TeV. It has the advantages of 100% duty cycle, plus large effective area and field of view, but a limited separation power between gamma rays and cosmic rays. The currently most advanced water Cherenkov gamma-ray detector is HAWC, composed of 300 water Cherenkov detectors located at an altitude of 4100 m at the Sierra Negra volcano, near Puebla (Mexico), covering 22,000 m 2 . It is sensitive to gamma rays between 500 GeV and 100 TeV, with a field of view of 15% of the sky, and daily coverage of 8.4 sr, or 67% of the sky (a region where dark matter searches using the HAWC data archive will be possible should new dSphs will be discovered in the future). Partial HAWC operations started in 2013, and the full detector was completed in March 2015.

Statistical Data Analysis
Advanced searches for dark matter annihilation or decay in dSphs with gamma rays rely on the distinct spatial and spectral features of the expected signals. We expect dark matter signal to be distributed morphologically according to dJ dΩ , and spectrally according to dN γ dE , and those PDFs are in general clearly distinguishable from those expected for background processes.
Regarding the use of the morphological information, the spatial coincidence of the signal with the position of the dSph would provide strong discrimination power, because we do not expect that gamma 7 of 27 rays can be produced at dSphs by any conventional astrophysical process. However, using the information of the morphology of the gamma-ray emission around the position of the dSph is more delicate, because such morphology is in general subject to relatively large uncertainties, and assuming an incorrect shape may bias the result of the search. In addition, the expected size of the dark matter halo is, for many of the known dSphs and for the considered gamma-ray instruments, consistent with point like sources, or at most slightly extended, which means that we can obtain no or little signal/background discrimination power from the use of the morphological information. All this is particularly true for dark matter annihilation, for which, due to the ρ 2 dependence of dJ dΩ , the expected signal is more compact and more affected by uncertainties on the details of the dark matter distribution within the halo. When looking for dark matter decay signal, on the other hand, such dependence is linear with ρ, which leads to less peaked and less uncertain morphologies.
The use of spectral information would be key for univocally attributing a dark matter origin to a detected gamma-ray signal, because in general, the features present in the spectra predicted for dark matter annihilation or decay cannot be produced by other conventional astrophysical processes. For instance, in the most extreme/luckiest case, the detection of gamma-ray spectral lines would be considered as unambiguous prove for the observation of dark matter annihilation or decay. Other processes, like creation of Standard Model particle pairs also produce distinct spectral features providing high discrimination power over backgrounds, such as the existence of sharp kinematic spectral cutoffs (see Figure 1). These considerations are general for all dark matter searches, independently of whether they are performed on dSphs or elsewhere. Searches in dSphs have the additional advantage that dark matter signals are, in principle, universal, any potential detection from a given dSph could be confirmed by looking for the same spectral features in the emission from other dSphs. Contrary to the case of dJ dΩ , uncertainties in dN γ dE can be considered negligible for a given annihilation/decay channel. This is the main reason why gamma-ray instruments utilize the spectral information not only for reinforcing the credibility of an eventual future detection, but also to increase the sensitivity of the search and therefore provide more constraining bounds to the dark matter nature in case of no detection.
Current dark matter searches using gamma rays are based on different implementations of the likelihood-ratio test [19], which we use to quantify the compatibility of the measured data (D) with different hypotheses, in particular with the null hypothesis (i.e., that no dark matter signal is present in D), through the associated p-value. Finding a sufficiently low p-value (by convention in the field p < 3 · 10 −7 ) for the observed data D under the null hypothesis assumption is usually referred to as detecting dark matter. In case of a positive detection, we can use the likelihood function to measure the dark matter physical parameters such as its mass, annihilation cross section, decay lifetime, and branching ratio to the different decay/annihilation channels (collectively represented here by the vector α). Conversely, if the null hypothesis cannot be excluded, we can use the likelihood function to set limits to the parameters α.
The likelihood function can be written in the following general form: where, apart from its dependence on α and D, we have made explicit that L can also depend on other, so-called, nuisance parameters (ν), for which we only know their likelihood function (normally constrained using dedicated datasets). In general, nuisance parameters represent quantities used in the computation of arrival direction present in the signal region, or J. One standard technique to eliminate the nuisance parameters when making statements about α is using the profile likelihood ratio test: whereα andν are the values maximizing L, andν the value that maximizes L for a given α. According to Wilks' theorem −2 ln λ P (α) is distributed, when α are the true values, as a χ 2 distribution with number of degrees of freedom equal to the number of components of α, independent of the value of ν. It is an extended practice in indirect dark matter searches with gamma rays to decrease the n-dimensional vector α of free parameters to a one-dimensional quantity α, by considering that gamma-ray production is dominated either by annihilation (α = σv , i.e., the velocity-averaged annihilation cross section) or by decay (α = τ −1 χ , i.e., the decay rate), and scanning over values of the dark matter particle mass (m χ ) and pure annihilation/decay channels (i.e., considering at each iteration 100% branching ratio to one of the possible SM particle pairs). For each scanned combination, Equation (11) reduces to a likelihood function of just one purely free (i.e., non-nuisance) parameter. In such a case, for instance, 1-sided 95% confidence level upper limits to α are taken as α UL 95 = α 2.71 , with α 2.71 found by solving the equation The data D can refer to N dSph different dSphs, in which case it is convenient to write the joint likelihood function as: where we have factorized the joint likelihood into the partial likelihood functions corresponding to each dwarf, and those subsequently into the parts corresponding to the gamma-ray observations (L γ ) and J-factor measurement (L J ), respectively; J l is the total J-factor (see Equation (7)) of the l-th considered dSph, which, as we have made explicit, is a nuisance parameter degenerated with α in L γ ; µ l represents the additional nuisance parameters different from J l affecting the analysis of the l-th dSph; D γ l represents the gamma-ray data of the l-th dSph, whereas D J l refers to the data constraining J l .
For each dSph, we may have N meas independent measurements, each performed under different experimental conditions, by the same or different instruments. That is, we can factorize the L γ term as: where we have omitted the index l referring to the dSph for the sake of clarity, and with µ k and D γ,k representing the nuisance parameters and data, respectively, referred to the k-th measurement.
For each observation of a given dSph under certain experimental conditions, L γ,k often consists of the product of N E × Np Poissonian terms (P) for the observed number of gamma-ray candidate events (N ij ) in the i-th bin of reconstructed energy and j-th bin of reconstructed arrival direction, times the likelihood term for the µ nuisance parameters (L µ ), with N E the number of bins of reconstructed energy and Np the number of bins of reconstructed arrival direction, i.e.: where the indexes l and k referring to the dSph and the measurement have been removed for the sake of a clear notation. The parameter of the Poissonian term is s ij + b ij , where s ij is the expected number of 9 of 27 signal events in the i-th bin in energy and the j-th bin in arrival direction, computable using αJ as we will see below; and b ij the corresponding contribution from background processes. D µ represents the data used to constrain the values of the nuisance parameters µ. We have made explicit that the uncertainties associated to µ can in principle affect both the computation of the signal and background contributions. For instance, uncertainties in the overall energy scale affect the computation of s ij , whereas uncertainties in the background modeling affect the computation of b ij . However, uncertainties affecting s ij are usually considered to be largely dominated by the uncertainty in the J-factor and the dependence of s ij on µ therefore ignored. Thus, s ij , is given by: where E ,p , E andp are the estimated and true energies and arrival directions, respectively; dΩ and dΩ infinitesimal solid angles containingp andp, respectively; T obs the total observation time; t the time along the observations; and IRF the instrument response function, i.e. IRF(E ,p |E,p, t) dE dΩ is the effective collection area of the detector times the probability for a gamma ray with true energy E and directionp to be assigned an estimated energy in the interval [E , E + dE ] andp in the solid angle dΩ (see more details below), at the time t during the observations. The integrals over E andp perform the convolution of the gamma-ray spectrum with the instrumental response, whereas those over E andp compute the events observed within the i-th energy bin (∆E i ) and the j-th arrival direction bin (∆p j ). It must be noted that, defining several spatial bins within the source produces relatively minor improvement in sensitivity to dark matter searches for not significantly extended sources (i.e., those well described by a point-like source, as it is the case for many dSphs) [20]. For significantly extended sources, on the other hand, using a too fine spatial binning makes the obtained result more sensitive to the systematic uncertainties in the dark matter spatial distribution within the dSph halo. Thus, a realistic optimization of Np based on sensitivity should balance the gain yielded by the use of more spatial information and the loss caused by the increase in the systematic uncertainty.
The IRF can be factorized as the product of the detector collection area A eff (T obs · A eff is often referred to as exposure), times the PDFs for the energy ( f E ) and incoming direction ( fp) estimators, i.e.: where, following the common practice, the (small) dependence of f E withp has been neglected. fp is often referred to as the point spread function (PSF). Finally, the likelihood for the total J-factor is usually written as: with log 10 J obs and σ J the mean and standard deviation of the fit of a log-normal function to the posterior distribution of the total J-factor [21]. Therefore, including L J in the joint likelihood is a way to incorporate the statistical uncertainty of J in the estimation of α. It is worth noting that, because α and J are degenerate, in order to perform the profile of L with respect to J it is sufficient to compute L γ vs α for a fixed value of J, which facilitates significantly the computational needs of the profiling operation (see details in footnote 12 of reference [22]). Including J obs systematic uncertainties is much more complex, since they depend mainly on our choice of the dark matter halo density profile function (e.g., NFW [23], Einasto [24], etc.), and there is no obvious way of assigning a PDF to that choice. Because of this, the impact of that uncertainty in 10 of 27 the bounds in α are usually roughly quantified by performing the likelihood analysis several times, each assuming different fitting functions, and comparing the results obtained for each of them. The PDF of the test statistic −2 ln λ P for the no-dark matter null hypothesis, i.e., when the true value of α is given by α true = 0, is needed for evaluating the significance of a possible signal detection. Computing upper limits to α, on the other hand, consists in finding the value of α true for which the integral of the PDF aboveα corresponds to the required confidence level. Estimating the PDF for −2 ln λ P with fast simulations is feasible (from a computational-demand standpoint) when the involved p-values are high enough so that they can be evaluated with a relatively low number of simulated datasets. In practice, however, results for dark matter searches using gamma-rays are generally computed assuming Wilks' theorem validity, and that −2 ln λ P is distributed as a χ 2 . The adoption of Wilks' theorem by all the experiments allows at least a direct comparison among their results. One should keep in mind, however, that the described statistical framework is also usually affected by the non-fulfillment of the conditions of validity of Wilks' theorem, at least because of two different reasons. First, because α is normally restricted to the physical region (i.e., to non-negative values), which produces over-coverage (i.e., the computed confidence interval contains the true value more often than the quoted confidence level) for negative background fluctuations, i.e., when the likelihood absolute maximum lies at the border of the physical region. This can be avoided by using the correct −2 ln λ P PDF for this situation [25]. Another way commonly used to partially mitigate this problem is to show the obtained result (e.g., the upper limit to α) in comparison to its PDF for the no-dark matter (α true = 0) hypothesis. Such PDF is estimated using fast simulations and/or pure-background datasets (such as those obtained by considering randomly selected directions as potential DM targets), and it is normally characterized by its median (referred to as the sensitivity of the measurement) and the bounds for some predefined (e.g., 68%, 95%, etc.) symmetric containment quantiles. By such comparison one can evaluate whether the obtained result is significantly incompatible with the α true = 0 hypothesis. The second violation of Wilks' theorem validity conditions affects the computation of confidence intervals (i.e., the PDF of the test −2 ln λ P for α true > 0). In this case, however, because α and J are degenerate in the likelihood function, the log-normal shape of the likelihood term L J (see Equation (18)) results in the loss of Gaussianity of the likelihood for α required by the Wilks' theorem.
As we will see in the next Section, the most common simplifications adopted in gamma-ray data analyses consist in ignoring the statistical and/or systematic uncertainties in J or in the background contribution to the signal region. Omitting these relevant uncertainties in general improves artificially the reported sensitivity and bounds obtained by the analysis, which must be taken into account when comparing results obtained under different assumptions.

Results
None of the different gamma-ray telescopes has obtained a significant detection in their search for dark matter signals from dSphs. Therefore, they provide results in the form of upper limits to the annihilation cross section or lower limits to the decay lifetime. In this section, I summarize the results obtained by the different considered instruments. In addition, I highlight and motivate the main analysis choices adopted by the different experiments as well as the differences with respect to the general framework described in Section 4, also summarized in Table 1.
In their 6-year-data search, the Fermi-LAT Collaboration applied their most developed data (re-)analysis, known as Pass 8. They subsequently searched for gamma-ray signals individually in 25 dSphs (including the classical and the ultra-faint ones discovered by the Sloan Digital Sky Survey [38]), and combined the 15 targets with better determined dark matter content. The dark matter distribution in each dSph was parameterized using the Navarro-Frenk-White (NFW) profile [23], constrained using the prescription by Martinez (2015) [39]. The dN γ dE average spectra for the different considered channels, on the other hand, were obtained from the PHYTIA-based [9] DMFIT package [40]. Table 1. Summary of dark matter searches with gamma-ray instruments. From left to right, columns show: Bibliographic reference; Instrument; Targets; Investigated decay and/or annihilation channels; dN/dE source; J-factor source; whether the following aspects have been included in the analysis: J-factor uncertainty, morphology of the source, restriction of α to physical region, statistical and systematic background uncertainties, determination of the true −2 ln λ P PDF; other relevant differences of the analyses with respect to the general framework (1: In Equation 19, assuming dφ/dE ∝ E −2 and energy resolution and bias disregarded; 2: In Equation (16), A eff dependence onp disregarded; 3: In Equation (16), effect of angular resolution disregarded (i.e., fp → δ(p −p )); 4: In Equation 28, f s assumed radially symmetric with respect to the center of the dSph). See main text for more details.
Cembranos et al. [42] Martinez [39] assuming NFW [23] and Burkert [43] × × × Fermi-LAT data statistical analysis follows Equations (13) and (14), with N dSph = 15 observed dSphs and N meas = 4 referring to the four independent datasets, each containing events with one of the four possible event direction reconstruction quality level, and hence each described by different IRF. They consider bins of reconstructed energy in the range between 500 MeV and 500 GeV and bins of incoming direction in a region of interest of 10 • × 10 • around the position of each dSph. The dominant background is produced by gamma rays from nearby sources (whose estimated energy and incoming direction are consistent, due to the finite angular resolution of the instrument, with being originated at the dSph dark matter halo), or by the diffuse gamma-ray component resulting from the interaction of cosmic rays with the interstellar medium or from unresolved sources of Galactic and extragalactic origin, depending on the particular dSphs being considered. The analysis does not explicitly treat the relevant background parameters µ in Equation (15) as nuisance parameters. Instead, the spectral parameters (e.g., normalization, photon index, etc.) of the different background sources are fixed using the following simplified method. The flux normalizations of the different background components are determined by means of a maximum-likelihood fit to the spacial and spectral distributions of the observed events, with the rest of spectral parameters fixed to the values listed in the updated third LAT source catalog [52]. Then, it is checked that the values of the background normalization factors obtained using this method do not change significantly by including an extra weak source at the locations of the dSph, which shows that the background are well-constrained by this procedure. Studies showed that the effect of the background uncertainty from this procedure contributed at a few percent of statistical uncertainty of the signal and are therefore safe to neglect.
In order to produce a result valid for arbitrary spectral shapes (i.e., arbitrary value of m χ and of the branching ratios to the different annihilation/decay channels), the Fermi-LAT Collaboration computes, for each observed dSph and bin of estimated energy ∆E i , the value of: as a function of Φ i , that is, the sum over the spatial bins of the L γ likelihood values within ∆E i . In order to obtain a set of generic L γ i values, they compute s ij (Φ i ) using Equation (16) assuming a power-law gamma-ray spectrum ( dΦ dE ∝ E −Γ ) of spectral index Γ = 2. The spatial distribution of gammas (which does not depend on the energy) is considered known and fixed, and given by the dJ dΩ curves obtained from the fit to the stellar kinematics to the different dSphs. Equation (15) can then be written in terms of the L γ i factors as: The values of L γ i vs Φ i for each of the analyzed dSphs were computed, tabulated and released by the Fermi-LAT Collaboration [53]. This allows any scientist to compute L γ for the dark matter model of their choice by just selecting the corresponding values of α, m χ , the total dN γ dE and J, and computing the corresponding Φ i values as: with dΦ dE given by Equation (10). We note that this approach allows to compute bounds on αJ with associated confidence level known only to a certain (unquantified) precision that depends on how similar are the investigated spectral shape and the one assumed when computing the values of L γ i (i.e., a power-law with Γ = 2 in the Fermi-LAT case). In addition, it should be stressed that such precision depends also on 14 of 27 the bb and τ þ τ − channels with expectation bands derived from the analysis of 300 randomly selected sets of blank fields. Sets of blank fields are generated by choosing random sky positions with jbj > 30°that are centered at least 0.5°from 3FGL catalog sources. We additionally require fields within each set to be separated by at least 7°. Our expected limit bands are evaluated with the 3FGL source catalog based on four years of PASS7 REPROCESSED data and account for the influence of new sources present in the six-year PASS8 data set.
Comparing with the results of Ackermann et al. [13], we find a factor of 3-5 improvement in the limits for all channels using six years of PASS8 data and the same sample of 15 dSphs. The larger data set as well as the gains in the LAT instrument performance enabled by PASS8 both contribute to the increased sensitivity of the present analysis. An additional 30%-40% improvement in the limit can be attributed to the modified functional form chosen for the J factor likelihood (3). Statistical fluctuations in the PASS8 data set also play a substantial role. Because the PASS8 six-year and PASS7 REPROCESSED four-year event samples have a shared fraction of only 20%-40%, the two analyses are nearly statistically independent. For masses below 100 GeV, the upper limits of Ackermann et al. [13] were near the 95% upper bound of the expected sensitivity band while the limits in the present analysis are within 1 standard deviation of the median expectation value. FIG. 1 (color). Constraints on the DM annihilation cross section at the 95% CL for the bb (left) and τ þ τ − (right) channels derived from a combined analysis of 15 dSphs. Bands for the expected sensitivity are calculated by repeating the same analysis on 300 randomly selected sets of high-Galactic-latitude blank fields in the LAT data. The dashed line shows the median expected sensitivity while the bands represent the 68% and 95% quantiles. For each set of random locations, nominal J factors are randomized in accord with their measurement uncertainties. The solid blue curve shows the limits derived from a previous analysis of four years of PASS7 REPROCESSED data and the same sample of 15 dSphs [13]. The dashed gray curve in this and subsequent figures corresponds to the thermal relic cross section from Steigman et al. [5].
FIG . 2 (color). Comparison of constraints on the DM annihilation cross section for the bb (left) and τ þ τ − (right) channels from this work with previously published constraints from LAT analysis of the Milky Way halo (3σ limit) [57], 112 hours of observations of the Galactic center with H.E.S.S. [58], and 157.9 hours of observations of Segue 1 with MAGIC [59]. Pure annihilation channel limits for the Galactic center H.E.S.S. observations are taken from Abazajian and Harding [60] and assume an Einasto Milky Way density profile with ρ ⊙ ¼ 0.389 GeV cm −3 . Closed contours and the marker with error bars show the best-fit cross section and mass from several interpretations of the Galactic center excess [16][17][18][19]. the PDF of the energy estimator and that, therefore, the range of investigated spectral shapes for which we can establish bounds within a certain precision using this technique is different for different instruments.
No significant gamma-ray signal from dSphs was found in the Fermi-LAT data, either individually in each dSph (the largest deviation from the null hypothesis is found for Sculptor, with −2 ln λ P = 4.3), or in the combined analysis (−2 ln λ P = 1.3). Some of the obtained exclusion limits are shown in Figure 2. This work represents the most constraining search for WIMP annihilation signals for the dark matter particle mass range below ∼1 TeV. As shown in the figure, the limits exclude the thermal relic cross section for m χ < 100 GeV in the case of annihilation into bb or τ + τ − pairs.
These results were combined with MAGIC observations of Segue 1, into the first coherent search for dark matter using several gamma-ray instruments [22]. Details about this work are provided below.
In a later work, the Fermi-LAT and the Dark Energy Survey (DES) collaborations also used the data from 6 years of observations to look for dark matter signals over a sample of 45 stellar systems consistent with being dSphs [55]. The search was performed shortly after the discovery of 17 of the considered dSph candidates, for which no reliable estimate of the dark matter content was available at the time. Because of this, all considered candidates were assumed to be point-like sources, and the J-factors for the non-confirmed dSphs estimated from a purely empirical scaling relation based on their heliocentric distance. For four of the examined dSphs, a 2σ discrepancy with the null hypothesis was found, which does not contradict significantly such hypothesis, particularly once the number of investigated sources, channels and masses is considered. Overall, the strategy of observing a set of not fully confirmed dSphs candidates, for which no reliable estimate of the J-factor exists yet is justified since a solid positive gamma-ray signal from any of the observed targets would have been considered a strong experimental evidence of dark matter annihilation or decay. In absence of such signal, however, the obtained limits are less robust than those from the 15 confirmed dSphs described above, which remain the reference in the field for the sub-TeV mass range.

Cherenkov Telescopes
Dark matter searches with Cherenkov telescopes have evolved from simple event-counting analyses to more complex maximum-likelihood analyses of optimized sensitivity thanks to the inclusion of the expected spectral and morphological features of the dark matter signals [56].
In the most basic version of the likelihood function, the nuisance parameters µ (see Equation (15)) are the b ij factors themselves. They are constrained by measurements in signal-free, background-control (or Off) regions with τ times the exposure of the signal (or On) region. A more complete analysis also includes the treatment of τ as a nuisance parameter, given that the latter is normally affected by significant statistic (σ τ,stat ) and systematic (σ τ,sys ) uncertainties. The statistical uncertainty comes from the fact that τ is often estimated from the data themselves (comparing the events observed in regions adjacent to the On and Off ones). The systematic uncertainty takes into account the residual differences of exposure between the Off and On regions, and it is normally assumed to be of the order of 1% for the current generation of Cherenkov telescopes [12]. It can be shown that this systematic uncertainty is the limiting factor to the sensitivity of the event-counting analyses for N On (τ+1)τ ∆τ 2 , i.e., between ∼ 10 4 and ∼ 2 × 10 4 events for τ in the typical range between 1 and 10, and 1% systematic uncertainty in τ (i.e., σ τ,sys = 0.01τ). Once we reach this number of observed events in the signal region, increasing the statistics of the dataset does not longer contribute to improve the sensitivity of the search.
The gamma-ray likelihood function for Cherenkov telescopes can thus be written as the product of Poisson likelihoods for the On and Off region times a Gaussian likelihood for τ, i.e.: with N On,ij and N Off,ij the number of observed events in the On and Off regions, respectively, in the i-th bin of reconstructed energy and the j-th bin of reconstructed arrival direction; and G an (often neglected) Gaussian PDF with mean the measured value τ obs and width σ τ = σ 2 τ,stat + σ 2 τ,sys . The considered energy range depends on the instrument (e.g., larger reflectors provide lower thresholds) and the dSph observation conditions (e.g., higher zenith angle observations imply higher threshold). For the current instruments and observed dSphs, the lowest energy bin starts between 80 and 800 GeV, whereas the highest one can reach up to between 10 and 100 TeV.
In the analysis of Cherenkov telescope data, the convolution of d 2 Φ dEdΩ · A eff with the PSF function fp needed to compute s ij according to Equation (16) is usually performed numerically through the analysis of Monte Carlo simulated events. We note that Equation (16) can be written as: withĀ eff,j the signal morphology-averaged effective area within spatial bin j, defined as:

of 27
A eff,j depends on the morphology of the gamma-ray emission ( dJ dΩ ), although not on its intensity, hence not on J. Therefore, for point-like sources observed with constant IRF at a given fixed direction,Ā eff,j is only a function of the energy. As a matter of fact, what normally is referred to as the effective area of a given Cherenkov telescope is the value ofĀ eff,j (E) for a circular spatial bin centered at the position of a point-like source (observed at low zenith distance under dark and good weather conditions), with radius optimized to maximize the signal-to-noise ratio. In practice,Ā eff,j (E) is computed with Monte Carlo simulations: Using a sample of simulated gamma rays with arrival directions distributed according to dJ dΩ and trajectories impacting uniformly in a sufficiently large area (A tot ) around the telescope pointing axis,Ā eff,j is computed scaling A tot by the ratio between events detected within spatial bin j and the total number of generated events. For reasons of economy of computational resources,Ā eff,j is computed in some of the analyses described here approximating dJ dΩ by a point-like source (i.e., by a delta function), even for the analysis of moderately extended dSphs. This approximation is less accurate the more extended the source. The bias introduced inĀ eff,j becomes relevant when the source extension is comparable to or bigger than the region for which A eff may be considered flat.

H.E.S.S
The first dark matter searches using observations of dSphs with the H.E.S.S telescopes were based on an event-counting analysis, with no attempt to use the expected spectral and morphological signatures in the search [57,58]. H.E.S.S also performed early searches on non-confirmed dSphs like Canis Major [59] or even globular clusters [60]. Their most recent searches use state-of-the-art analysis techniques like the one described in Section 3, and are based on observations of Sagittarius (∼90 h), Coma Berenices (∼9 h), Fornax (∼6 h), Carina (∼23 h) and Sculptor (∼13 h) dSphs, where H.E.S.S has searched for both continuum [41] and line-like [44] dark matter spectra. We will concentrate in these two latter works.
There are significant differences in the high-level maximum-likelihood analyses used by H.E.S.S in their searches for continuum and spectral lines. In the search for continuum spectra the dN γ dE was taken from analytical parameterizations [42], generally valid up to dark matter mass of 8 TeV; and the J-factors were estimated using the prescription by Martinez (2015) [39], assuming alternatively cuspy NFW [23], and cored Burkert [43] profiles. In the search for spectral lines, on the other hand, the dN γ dE is trivial, and the J-factors were taken from the work by Geringer-Sameth et al. (2014) [7], which assumes a Zhao-Hernquist dark matter density profile [61]. The case of Sagittarius dSph was treated separately in both works, given that this galaxy is likely affected by tidal disruption [62], and therefore the J-factor calculation is subject to comparatively larger systematic uncertainties, not included in the likelihood analysis.
There are also slight differences in the likelihood function used by H.E.S.S in the continuum and spectral line dark matter searches. In both cases they use the likelihood function of Equation (22) without including the term accounting for the uncertainty in τ. In the case of continuum spectra, only one (Np = 1) circular spatial bin centered at the position of each dSph was considered, whereas for spectral lines Np = 2 or 3 (depending on the size of the considered dSph), concentric 0.1 • -width ring-like spatial bins were used. The reason for this difference must be purely historical (given how the expected and measured spatial information is essentially common to both searches), probably in an attempt to increase the sensitivity by including more information in the likelihood analysis. The drawback of this approach, has already been discussed in Section 4: It can introduce a bias in s ij , with an unquantified effect in the final sensitivity to α. In both analyses, for the computation of s ij following Equation (16), the dependence of the effective area withp within the signal region and the effect of the PSF are ignored. We note that these two simplifications require opposite conditions: A eff can be better approximated by a constant value for smaller signal regions, i.e., smaller dSphs, whereas the effect of the angular resolution in the distribution of measured events These exclusion curves are subject to uncertainties also in the particle physics side. In order to illustrate the particle physics uncertainties, Fig. 7 shows the limits obtained from the Sgr alone for different annihilation channels, assuming the NFW halo profile as reference. The strongest bounds stronger. Similarly, one can interpret the two-muon final state channel constraints: although weaker than the τ þ τ − at high mass, since the photon yield is lower, it becomes comparable at low mass. Note that this channel, which has the hardest spectrum among the considered channels,  Poisson realizations of both NON,ij and NOFF,ij. The mean sensitivity as well the statistical 68% (±1 ) and 95% (±2 ) containment bands are also plotted. A limit of h vi . 3 ⇥ 10 25 cm 3 s 1 is reached in the mass range of 400 GeV to 1 TeV. The combination of all five galaxies allows an improvement in the constraints up to a factor of 2 around 600 GeV with respect to individual galaxies.
Note that, at certain DM masses, the combined limit becomes worse than some individual limits, and the combined limit without Sagittarius becomes more constraining than the combined one that includes Sagittarius. This is due to the statistical effect of adding an individual data set with large negative fluctuations (or excesses) and large expected signal around those energies (case of Sculptor at ⇠ 350 GeV, Coma Berenices at ⇠ 2 TeV, and Sagittarius at ⇠10 TeV), to large data sets with positive fluctuations and smaller or similar expected signal (case of Carina and Sagittarius at ⇠ 350 GeV, Carina, Sculptor and Sagittarius at ⇠ 2 TeV, and Coma Berenices at ⇠10 TeV). The limits of the individual data sets will be highly overestimated, while the combination with the large data sets will push the combined limits to less constraining values.

Limits on pure WIMP models
In this section 'pure WIMP' models are briefly introduced and the distinctive shapes of their gamma-ray annihilation spectra are discussed. The results for the limits on this class of models are then presented.  is smaller for larger dSphs. The effect in the final result of adopting these two simplifications is not quantified.
H.E.S.S found no significant gamma-ray signal in the observed dSphs, considered either individually or collectively, for any of the assumed emission spectra (continuum or spectral line). The maximum observed deviations from the null hypothesis are ∼2.6σ for the continuum spectra search in Fornax, and ∼1.2σ for the spectral-line search in Sagittarius. The exclusion limits for the annihilation cross-section for continuum spectra (see Figure 3-left) peak at dark matter masses of around 1-2 TeV, depending on the considered channel. Assuming a NFW density profile, the strongest constraint is provided by Sagittarius dwarf, with σv UL 95 ∼ 2 × 10 −23 cm 3 s −1 for a combination of W + W − and ZZ annihilation channels. The bounds resulting from the combination of all the observed dSphs are only marginally better because Sagittarius has, under the NFW-profile assumption, the largest by far J-factor among the considered dSph, and because it has been observed by H.E.S.S for significantly longer time than the rest of the dSphs. However, given that the value of the J-factor for Sagittarius is affected by large systematic uncertainties (on account of the possibility that the system is affected by tidal disruption), H.E.S.S has also provided constraints obtained from the combination of all the other dSphs, which results in the limit σv UL 95 ∼ 10 −22 cm 3 s −1 , for the same annihilation channel. In the case of the search for spectral lines, limits do not depend significantly on the inclusion or not of Sagittarius (see Figure 3-right), since with the newer approach in the evaluation of the J-factor used in this work, the limits are dominated by Coma Berenices results. In the mass range between 400 GeV and 1 TeV, the obtained limit to the velocity-averaged cross section is σv UL 95 ∼ 3 × 10 −25 cm 3 s −1 .

MAGIC
MAGIC performed early dark matter searches using observations of the dSphs Draco [63], Willman 1 [64] and Segue 1 [65]. These searches had a relatively poor sensitivity, due to the fact that they were based on one-telescope observations and simple event-counting data analysis. With the addition of the second telescope, MAGIC dark matter search strategy was based on deep observations (∼160 h) of the   the analysis of 1000 realizations of the null hypothesis (h v i = 0 (for both ON and background regions) generated from backg exposures as for the real data, and J-factors assumed as nu likelihood function. All bounds are consistent with the no-det  dSph with the highest J-factor known at that moment, namely Segue 1 [45], and the use for the first time by Cherenkov telescopes, of advanced maximum-likelihood analysis techniques. In addition, MAGIC Segue 1 observations were part of the aforementioned first multi-instrument combined search, together with data from Fermi-LAT [22], a work that I will discuss later in more detail. After that, MAGIC initiated a diversification of observed targets, starting by ∼100 h of observations of the Ursa Major II dSph [46].
MAGIC dark matter searches in Segue 1 and Ursa Major II follow essentially the same data selection, calibration and processing procedures, but contain significant differences in several elements of their high-level analysis. The gamma-ray average spectra per annihilation reaction ( dN γ dE ) were obtained from the parameterization by Cembranos et al. (2011) [42] in the case of Segue 1, and the PPPC 4 DM ID computation [47] in the case of Ursa Major II. The spectra provided by these two works do not differ significantly for the considered energy range. The J-factor for Segue 1 was computed solving the Jeans equation assuming an Einasto density profile [24], and for Ursa Major II was taken from Geringer-Sameth et al. (2014) [7].
The likelihood function used by MAGIC for dark matter searches has also evolved over the years. For the observations of Segue 1 they used, instead of Equation (22), the following unbinned likelihood function: where the uncertainty on τ is ignored, only one spatial bin is considered, and the energy-wise product of Poisson terms is substituted by a global Poisson term for the total number of observed events (N On ), times the joint likelihood for the observed values of estimated energies. The latter is computed as the product of the PDF for the reconstructed energy f s+b (E ) evaluated at each observed E , where f s+b = 1 s+b (s f s + b f b ), with f s and f b the PDFs for the reconstructed energies for signal and background events, and s (the free parameter) and b (a nuisance parameter) the total expected number of signal and background events, respectively. f s is the normalized convolution of the gamma-ray spectrum with the IRF, i.e.: withĀ eff (E) computed following Equation (24). f b , on the other hand, is modeled using the data from one or several Off regions. This approach presents the drawback of neglecting the statistical and systematic uncertainties in the background spectral shape. In comparison, in the binned version of the likelihood function (Equation (22)) the statistical uncertainty is taken into account by the inclusion of the nuisance parameters b ij and τ. This unbinned analysis hence typically produces results that are several tens of percent artificially more constraining than the binned one. Another important difference of the MAGIC Segue 1 analysis with respect to the general framework is that it does not include statistical uncertainties in the J-factor. This was justified by the fact that the bounds to α scale with 1/J, and therefore the provided results allow the computation of the limits for any other J-factor value (provided dJ dΩ is kept fixed). This argument is valid only for single-target observations, but not for results obtained combining observations from different dSphs with different J values and uncertainties. Another main difference between this analysis and the general framework is the treatment of the cases when the valueα maximizing the likelihood lies outside the physical region, i.e.,α < 0. For those cases, the 95% confidence limit on α was computed as α UL 95 = α 2.71 −α, with α unrestricted (i.e., allowed to take negative values) during the likelihood maximization process. With this prescription, the limit obtained for any negative fluctuation in the number of excess events is equal to the limit for zero excess events (i.e., the sensitivity), at the expense of some over-coverage (i.e., the bounds are conservative).
In the analysis of Ursa Major II data, MAGIC used the general analysis framework described in Section 4, with binned likelihood, statistical uncertainties in the J-factor considered, and α restricted to positive values. In addition, for the first time in the analysis of Cherenkov telescope data, the Off/On exposure ratio τ in Equation (22) was considered a nuisance parameter, taking into account both its statistical and systematic (σ τ,sys = 1.5%) uncertainties, thus providing more realistic results.
MAGIC found no significant gamma-ray signal in the observations of Segue 1 or Ursa Major II. This was translated into limits to the dark matter annihilation cross section (and decay lifetime), assuming different dark matter induced gamma-ray production mechanisms. Using Segue 1, MAGIC carried out a systematic search for annihilation and decay processes, looking for the continuum spectra from production of bb, tt, µ + µ − , τ + τ − , W + W − and ZZ pairs, spectral lines from γγ and γZ channels, and other spectral features such as those produced by virtual internal bremsstrahlung emission (XXγ) and gamma-ray "boxes" (ΦΦ → γγγγ). With Ursa Major II data, the searches were limited to annihilation into bb, µ + µ − , τ + τ − and W + W − pairs. Figure 4 shows the results for annihilation into bb pairs obtained from each of the observed dSphs (there is no MAGIC-only combined result). The obtained limits are in general within the 68% containment region expected for the null hypothesis, except for the low mass range m χ 300 GeV in the case of Segue 1, where they stay nevertheless within the 95% containment region. 95% confidence level upper limits to the annihilation cross-section of dark matter particles into bb pairs reach σv UL 95 ∼ 5 × 10 −24 cm 3 s −1 for m χ ∼ 2 TeV in the case of Segue 1, and σv UL 95 ∼ 2 × 10 −23 cm 3 s −1 in the case of Ursa Major II. Segue 1 observations were also used to constrain the lifetime of m χ ∼ 20 TeV particles decaying into bb pairs to be larger than τ  simple event counting analysis approach. More recently, they analyzed their full datasets and combined them using advanced analysis techniques [48].
In this latter work, the average gamma-ray spectra ( dN γ dE ) for the investigated dark matter annihilation channels were taken from the PPPC 4 DM ID computation [47], and the differential J-factors from Geringer-Sameth et al. (2014) [7]. For the high-level, statistical data analysis, VERITAS used a test statistic equivalent to the ratio of the following likelihood function [68], namely: This likelihood function is similar to the one used by MAGIC in the Segue 1 analysis (Equation (25)). They are both unbinned simplified versions of the general likelihood function for Cherenkov telescopes shown in Equation (22). With respect to the MAGIC Segue 1 likelihood function, in Equation (27) the external Poisson term for the total number of observed events is omitted, and the event-wise term consists in the evaluation of the 2-dimensional PDF for the measured energy E and the angular separation θ between the measured arrival direction and the dSph center. We remind the reader that f s+b = 1 s+b (s f s + b f b ). In the 2-dimensional case, assuming that the convolution of the gamma-ray distribution with the IRF is radially symmetric with respect to the center of the dSph (i.e., the dependence onp reduces to a dependence on θ ), then f s (E , θ ) is given by: Only events in an On region defined by a maximum distance of θ cut = 0.17 • from the center of the dSphs are considered, and the dependence of the effective area A eff on the arrival directionp for events passing such cut is ignored. The dependence of f b on E is modeled by smearing the distribution of E measured for events of the background-control (Off) region, whereas the spatial distribution is assumed to be uniform within the On region. Both b and f b are fixed during the likelihood maximization, i.e., no statistical or systematic uncertainties in the background estimation are considered. Moreover the J-factor uncertainty is not included in the likelihood. Instead, the effect of the uncertainty in J is quantified by repeating the limit calculation over an ensemble of dark matter halo realizations using, for each dSph, halo parameter values randomly chosen from their inferred PDFs, and reporting the 68% confidence level containment quantiles of the obtained distribution of results 1 . However, the main reported result in this case is still the median of such distribution, which is only sensitive to the central J-factor and not to its uncertainty, producing limits a factor ∼ 2 more constraining than if J was considered a nuisance parameter.
A possible advantage of the use of the likelihood function of Equation (27) is that it allows a relatively simple estimation of the PDF for the associated −2 ln λ P test statistic [68] directly from the data and without relying on the validity of the conditions of the Wilks' theorem. This is so, because −2 ln λ P can be expressed as the sum of two random variables (those corresponding to the signal and background contributions to −2 ln λ P , respectively), which, for the likelihood function of Equation (27), are distributed according to a compound Poisson distribution. VERITAS results are hence robust in the sense that have a well determined confidence level under the assumption that the likelihood function was correct.
VERITAS has not found evidence of dark matter signals from neither of the four considered dSphs individually, or combined in a joint analysis. The null-hypothesis significance is well within the ±2σ quantile, for all considered targets, annihilation channels (uū, dd, ss, bb, tt, e + e − , µ + µ − , τ + τ − , W + W − , ZZ and hh) and m χ values, except for m χ ≥ 5 TeV dark matter particles annihilating into γγ in Draco dSph. In this latter case, a negative fluctuation slightly below −2σ is observed, which is not incompatible with purely statistical fluctuations, or could be alternatively explained by unaccounted systematic uncertainties in the background estimation. Figure 5 shows VERITAS limits to the annihilation cross-section into bb and τ + τ − pairs, compared with other limits from dSph observations by other gamma-ray instruments. The constraints reach σv UL 95 ∼ 10 −23 cm 3 s −1 at m χ ∼ 1 TeV for bb, and σv UL 95 ∼ 3 × 10 −24 cm 3 s −1 at ∼300 GeV for τ + τ − annihilation channels, respectively.

HAWC
HAWC has searched for dark matter annihilation and decay signals in 15 dSphs observed during 507 days between November 2014 and June 2016 [49]. They computed the average gamma-ray spectra per annihilation or decay event ( dN γ dE ) using the PYTHIA v8.2 simulation package [50], and the J-factors using the CLUMPY software package [51], assuming NFW [23] dark matter density profiles. The searches were carried out using the binned likelihood function described in Equation (15). Data were binned in reconstructed energy E (referred to as f hit in HAWC publications [14]) covering the range between 500 GeV and 100 TeV, and in reconstructed arrival directionp , covering an area of 5 • radius around each of the analyzed dSphs. The computation of the signal events s ij in each bin was performed using Monte Carlo simulations of the whole observations, assuming point-like sources and a reference value of α, and scaling the result for any other needed value, which is equivalent to using Equation (16).
No nuisance parameters accounting for uncertainties in the background estimation were considered, i.e., no L µ term was included in the Equation (15) likelihood function. The values b ij were estimated from the measured number of events in the same bin of local (or detector) coordinates at times when such coordinates do not correspond to any of the analyzed dSphs or any known HAWC sources. Measured background rates at each local spatial bin were then normalized using the all-sky event rate measured in 2-hour intervals. Using this method, the statistics used for background estimation correspond to an larger than θ max , where the DM halo is assumed to end. We impose this physically motivated constraint on the J-and D-factor uncertainty calculations, resulting in a one-side uncertainty. For the combined limit uncertainties, we use the uncertainties corresponding to Segue 1 (42% for annihilation cross-section limits and 38% for decay lifetime limits) since it is one of the strongest sources that is driving the limits. Though it would have been better to calculate and use these uncertainties for Triangulum II, the required information is not yet available. larger than θ max , where the DM halo is assumed to end. We impose this physically motivated constraint on the J-and D-factor uncertainty calculations, resulting in a one-side uncertainty. For the combined limit uncertainties, we use the uncertainties corresponding to Segue 1 (42% for annihilation cross-section limits and 38% for decay lifetime limits) since it is one of the strongest sources that is driving the limits. Though it would have been better to calculate and use these uncertainties for Triangulum II, the required information is not yet available. 8 Figure 6. The 95% confidence level upper limits to the annihilation cross-section of dark matter particles annihilating into bb (left) and τ + τ − (right) pairs, from HAWC observations of dSphs (black solid line). Results from other gamma-ray instruments are also shown (see legend for details), as well as the median and 65% and 95% symmetric quantiles of the distribution of limits obtained under the null hypothesis. Figure reproduced with permission from reference [49], ©AAS.
Off/On exposure ratio factor of τ = 30-300 [70], and the related statistic uncertainties (included in the case of Cherenkov telescopes by the second Poisson term in Equation (22)), can therefore be safely neglected. However, the effect of the systematic uncertainty associated to this method is not quantified or taken into account in the analysis. In addition, similarly to the case of VERITAS, HAWC does also not include in the maximum likelihood analysis the statistical uncertainty in the J-factor, i.e., they ignore the L J term in Equation (13). They do quantify the impact on the limits caused by the consideration of the dSphs as point-like sources and by several detector effects not perfectly under control in the Monte Carlo simulations used for calibrating the detector.
HAWC has not found gamma-rays associated to dark matter annihilation or decay from the examined dSphs, considered either individually or collectively. The significance of rejection of the null hypothesis for all considered targets, channels (bb, tt, τ + τ − , W + W − and µ + µ − ), and m χ values (between 1 and 100 TeV) is within 2σ, except for few marginally larger negative fluctuations. Figure 6 shows the limits to the annihilation cross section obtained by HAWC for the bb and τ + τ − annihilation channels, compared to limits obtained by other gamma-ray instruments. Limits reach σv UL 95 ∼ 10 −23 cm 3 s −1 at m χ ∼ 3 TeV for bb, and σv UL 95 ∼ 2 × 10 −24 cm 3 s −1 at ∼1 TeV for τ + τ − annihilation channels, respectively. For decay, lower limits to the decay lifetime were set to τ LL 95 χ ∼ 3 × 10 26 s for the 100 TeV mass dark matter particle decaying into bb pairs or τ LL 95 χ ∼ 10 27 s for decaying into τ + τ − pairs.

Multi-Instrument Searches
Following Equations (13) and (14), MAGIC and Fermi-LAT have computed a multi-target, multi-instrument, joint likelihood, producing the first coherent joint search for gamma-ray signals from annihilation of dark matter particles in the mass range between 10 GeV and 100 TeV [22]. The data used in this work correspond to the Fermi-LAT 6-years [21] and the MAGIC Segue 1 [45] observations discussed earlier in Sections 5.1 and 5.2.2, respectively. MAGIC analysis was slightly adapted to match LAT conventions, in the following aspects: (i) The determination of the J-factor; (ii) the treatment of the statistical uncertainty of J through the L J term in Equation (13); and (iii) the treatment of the cases in which the limits lie outside the physical (α ≥ 0) region.
The MAGIC/Fermi-LAT combined search for dark matter did not produced a positive signal, but it allowed setting global limits to the dark matter annihilation cross section and, for the first time, a  (Table 1) are considered as described in Section 3.2. The thin-dotted line, green and yellow bands show, respectively, the median and the symmetrical, two-sided 68% and 95% containment bands for the distribution of limits under the null hypothesis (see main text for more details). The red-dashed-dotted line shows the thermal relic cross-section from Ref. [54].
this magnitude would be expected in 5% of the experiments under the null hypothesis and is therefore compatible with random fluctuations.
As expected, limits in the low and high ends of the considered mass range are dominated by Fermi -LAT and MAGIC observations, respectively, and the combined limits coincide with the individual ones. The combination provides a significant improvement in the range between ⇠1 and ⇠100 TeV (for bb and W + W ) or ⇠0.2 and ⇠2 TeV (for ⌧ + ⌧ and µ + µ ),  (Table 1) are considered as describ line, green and yellow bands show, respectively, the median and the sym containment bands for the distribution of limits under the null hypothesi The red-dashed-dotted line shows the thermal relic cross-section from Ref this magnitude would be expected in 5% of the experiments u is therefore compatible with random fluctuations.
As expected, limits in the low and high ends of the conside by Fermi -LAT and MAGIC observations, respectively, and t with the individual ones. The combination provides a significa between ⇠1 and ⇠100 TeV (for bb and W + W ) or ⇠0.2 and ⇠ -9 -  Figure 7 shows the 95% confidence level limits to the cross-section of dark matter particles of mass in the range between 10 GeV and 100 TeV annihilating into bb and τ + τ − pairs. The obtained limits are the currently most constraining results from dSphs, and span the widest interval of masses, covering the whole WIMP range. In the regions of mass where Fermi-LAT and MAGIC achieve comparable sensitivities, the improvement of the combined result with respect to those from individual instruments reaches a factor ∼ 2.
This approach is applicable to all the high-energy gamma-ray instruments (and also to high energy neutrino telescopes, with slight modifications in Equation (16) to account for the oscillations). The so-called Glory Duck working group has initiated an activity aimed at the combination of all dark matter searches performed with Fermi-LAT, H.E.S.S, MAGIC, VERITAS and HAWC using observations of dSphs [71]. Each collaboration will analyze their own datasets and will provide the likelihood values as a function of the free parameter α (i.e., the terms L γ,k in Equation (14)) for the different considered annihilation channels and m χ values, for their combination and J-factor profiling through Equation (13). Likelihood values from the different instruments will be computed using the same conventions for the computation of the gamma-ray spectra and the J-factors, as well as the same statistical treatment of the data, most notably a common consideration of all relevant uncertainties by the inclusion of the corresponding nuisance parameters in the likelihood functions. While in principle foreseen only for the combination of gamma-ray data in the search of annihilation signals, this work could pave the path for other combined searches, such as searches for decay signals, the inclusion of other kinds of targets or even extending the searches to include also results from neutrino telescopes. This approach will ensure that all the combined individual results will be directly comparable among them, and will produce the legacy result of the dark matter searches using the current generation of gamma-ray instruments.