Limits to dark matter annihilation cross-section from a combined analysis of MAGIC and Fermi-LAT observations of dwarf satellite galaxies

We present the first joint analysis of gamma-ray data from the MAGIC Cherenkov telescopes and the Fermi Large Area Telescope (LAT) to search for gamma-ray signals from dark matter annihilation in dwarf satellite galaxies. We combine 158 hours of Segue 1 observations with MAGIC with 6-year observations of 15 dwarf satellite galaxies by the Fermi-LAT. We obtain limits on the annihilation cross-section for dark matter particle masses between 10 GeV and 100 TeV - the widest mass range ever explored by a single gamma-ray analysis. These limits improve on previously published Fermi-LAT and MAGIC results by up to a factor of two at certain masses. Our new inclusive analysis approach is completely generic and can be used to perform a global, sensitivity-optimized dark matter search by combining data from present and future gamma-ray and neutrino detectors.

A natural step forward within this collective effort is the combination of data from different experiments and/or observational targets, which allows a global and sensitivityoptimized search [38]. This can be achieved in a relatively straightforward way since, for a given DM particle model, a joint likelihood function can be written as the product of the particular likelihood functions for each of the measurements/instruments. One advantage of such an approach is that the details of each experiment, such as event lists or instrument response functions (IRFs), do not need to be combined or averaged.
In this paper, we present a new global analysis framework for DM searches, applicable to observations from gamma-ray and neutrino instruments, and the results of applying it to the MAGIC and Fermi -LAT observations of dSphs.
DSphs are associated with the population of Galactic DM sub-halos, predicted by the cold DM structure formation scenario and reproduced in N-body cosmological simulations (e.g. [39,40]), that have accreted enough baryonic mass to form stars (other sub-halos may remain completely dark). DSphs have very high mass-to-light ratios, being the most DM-dominated systems known so far [41]. They have the advantage of being free of other sources of gamma-ray emission and have been identified predominantly at high Galactic latitudes, where diffuse astrophysical foregrounds from the Milky Way are lowest. Because they are relatively close, they are expected to appear as point-like or marginally extended sources for gamma-ray and neutrino telescopes, with relatively high DM annihilation fluxes. The stellar kinematics of these systems can be used to determine their DM distribution and its uncertainty using a common methodology [42][43][44]. These measurements enable the combination of dSph observations that constrains models of DM annihilation or decay within the dSph population as a whole.
This article is organized as follows: Section 2 describes the MAGIC and Fermi -LAT instruments, and the data sets used in our study. The global analysis framework is described in Section 3. The results of applying the analysis to MAGIC and Fermi -LAT data are presented in Section 4. Finally, in Section 5 we discuss those results and present our conclusions.

The MAGIC telescopes
The Florian Goebel Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes are located at the Roque de los Muchachos Observatory (28.8 • N, 17.9 • W; 2200 m above sea level), on the Canary Island of La Palma, Spain. MAGIC is a system of two telescopes that detect Cherenkov light produced by the atmospheric particle showers initiated by cosmic particles entering the Earth's atmosphere. Images of these showers are projected by the MAGIC reflectors onto the photo-multiplier tube (PMT) cameras, and are used to reconstruct the calorimetric and spatial properties, as well as the nature of the primary particle. Thanks to its large reflectors (17 meter diameter), plus its high-quantum-efficiency and low-noise PMTs, MAGIC achieves high sensitivity to Cherenkov light and low energy threshold [45]. The MAGIC telescopes are able to detect cosmic gamma rays in the very high energy (VHE) domain, i.e. in the range between ∼50 GeV and ∼50 TeV.
For this work we use MAGIC data corresponding to 158 hours of observations of the satellite galaxy Segue 1 [32], the deepest observations of any dSph by a Cherenkov telescope to date. 9 The data were taken between 2011 and 2013, with observations before, during, and after a major hardware upgrade [46]. This resulted in the necessity of defining 4 different observation periods, each described by a different set of IRFs. Observations were performed in the so-called wobble mode, which allows simultaneous observations of the target and the background control regions. For that purpose, two different positions ∼0.4 • away from the position of Segue 1 were tracked alternatively for 20 min runs. For each of these positions, the spectral shape of the residual background is slightly different, and needs to be modeled independently. Therefore, for MAGIC, we consider 8 independent data sets, each consisting of the gamma-ray candidate events plus the corresponding IRFs and residual background models. The MAGIC data are analyzed with a one-dimensional unbinned likelihood fit to the energy distribution within a signal (or ON) region of radius 0.122 • around the center of Segue 1 (optimized for a source with the angular size of Segue 1). The integral of the DM emission template within the ON region defines the model for the expected signal distribution as a function of energy. See Section 3.2 and Ref. [32] for more details about the MAGIC data analysis.

The Fermi -LAT
The Fermi -LAT is a pair-conversion telescope that is sensitive to gamma rays in the energy range from 20 MeV to more than 300 GeV [47]. With its large field-of-view (2.4 sr), the LAT is able to efficiently survey the entire sky. Since its launch in June 2008, the LAT has primarily operated in a survey observation mode that continuously scans the entire sky every 3 hours. The survey-mode exposure coverage is fairly uniform with variations of at most 30% with respect to the average exposure. The LAT point-source sensitivity, which is dependent on the intensity of diffuse backgrounds, shows larger variations but is relatively constant at high Galactic latitudes (|b| > 10 • ). More details on the on-orbit performance of the LAT are provided in Ref. [48].
The Fermi -LAT data set used in this work corresponds to 6 years of observations of 15 dSphs. We analyzed the data with the latest (Pass 8) event-level analysis [34], using the Fermi Science Tools version 10-01-01 and the P8R2 SOURCE V6 IRFs. DM signal morphological and spectral templates are used in a three-dimensional likelihood fit to the distribution of events in reconstructed energy and direction within the 10 • × 10 • dSph region-of-interest (ROI). The model for each ROI contains the dSph DM intensity template, templates for Galactic and isotropic diffuse backgrounds, and point sources taken from the third LAT source catalog (3FGL) [49] within 15 • of the ROI center. A broadband fit in each ROI is performed fitting the normalizations of the Galactic and isotropic diffuse components and 3FGL sources that fall within the ROI boundary. After performing the broadband fit, a set of likelihoods is extracted for each energy bin by scanning the flux normalization of a putative DM source modeled as a power law with spectral index 2 at the location of the dSph. Tables with likelihood values versus energy flux for each energy bin are produced for all considered targets. The likelihood tables used for the present work are taken from Ref. [34] and can be found in the corresponding online material. 10 These tables allow the computation of jointlikelihood values for any gamma-ray spectrum, and are used as input to the present analysis (see Section 3.2 for more details).

Dark Matter annihilation flux
The gamma-ray (or neutrino) flux produced by DM annihilation in a given region of the sky (∆Ω) and observed at Earth is given by: where m DM is the mass of the DM particle, σv the thermally-averaged annihilation crosssection, dN/dE the average gamma-ray spectrum per annihilation reaction (for neutrinos this term includes the oscillation probability between target and Earth), and is the so-called astrophysical factor (or simply the J-factor), with ρ being the DM density, and the integrals running over ∆Ω and the line-of-sight (l.o.s.) through the DM distribution.
Using the PYTHIA simulation package version 8.205 [50], we have computed the average gamma-ray spectrum per annihilation process (dN/dE) for a set of DM particles of masses between 10 GeV and 100 TeV (i.e. in the WIMP range), annihilating into SM pairs bb, W + W − , τ + τ − and µ + µ − . For each channel and mass, we average the gamma-ray spectrum resulting from 10 7 decay events of a generic resonance with mass 2 × m DM into the specified pairs. For each simulated event, we trace all the decay chains, including the muon radiative decay (µ − → e −ν e ν µ γ, not active in PYTHIA by default), down to stable particles. We take J-factors around the analyzed dSphs from Ackermann et al. [34], who follow the approach of Ref. [42]. The DM distributions in the halos of the dSphs are parameterized following a Navarro-Frenk-White (NFW) profile [51]: where r s and ρ 0 are the NFW scale radius and characteristic density, respectively, and are determined from a fit to the dSph stellar surface density and velocity dispersion profiles. The properties of the dSphs used in our analysis, including the J-factors and their uncertainties, are summarized in Table 1. We quote the measured J-factors (J obs ) for a reference integrated radius of 0.5 • from the halo center in all cases, which encompasses more than 90% of the annihilation flux for our dSph halo models (which have halo scale radii between 0.1 • and 0.4 • ). We note that Table 1 together with Equations 3.2 and 3.3 allow the computation of the J-factors for any other considered ∆Ω, and therefore the DM emission templates. In Ref. [34] it has been shown how the Fermi-LAT limits can vary by up to ∼ 35% (for a 100 GeV DM mass, and decreasing with the mass value) by assuming different parameterizations for the DM density profile in dSphs. Using Equation 3.1 together with the values of dN/dE and J obtained as detailed in the previous paragraphs, we compute morphological and spectral intensity templates for the DM emission in each dSph. Folding these templates with the response of the MAGIC and LAT instruments, we compute the expected count distribution as a function of the measured energy and position within the observed field of each dSph.  Table 1: DSphs used in the present analysis and their main properties: Name, Galactic longitude and latitude, distance to Earth, angular size of the DM halo scale radius, and J-factor (with statistical uncertainty) assuming an NFW density profile and integrated to a radius of 0.5 • from the dSph center.

Likelihood analysis
For each considered annihilation channel and DM particle mass, we compute the profile likelihood ratio as a function of σv (see, e.g. [52]): with D representing the data set and ν the nuisance parameters. 11 σv andν are the values maximizing the joint likelihood function (L), andν the value that maximizes L for a given value of σv . The likelihood function can be written as: where the index i runs over the different targets (dSphs in our case); J i is the J-factor for the corresponding target (see Equations 3.2 and 3.3); µ i denotes any nuisance parameters additional to J i ; and D i is the target-related input data. J is the likelihood for J i , given measured log 10 (J obs,i ) and its uncertainty σ i [34]: The likelihood function for a particular target (L i ) can in turn be written as the product of the likelihoods for different instruments (represented by the index j), i.e.: where µ ij and D ij represent the nuisance parameters and input data set for the given target i and instrument j. Equations 3.4, 3.5 and 3.7 are generic, i.e. they are valid for any set of instruments and observed targets. 12 In our case, the likelihood for a given target i observed by the Fermi -LAT (j ≡ F ) is computed as: with k running over energy bins, and The values of L iF k vs. EΦ corresponding to 6 years of observations of each of the considered dSph are tabulated and released by the Fermi -LAT Collaboration [33,34]. For the case of MAGIC (j ≡ M ), the likelihood corresponding to a given target i can be written as: with the index k running over N = 8 different data sets (each described by a different IRF, see Section 2.1). The likelihood for a given data set follows the method described in Refs. [38] and [32] and can be written as (target, experiment and data set indices are omitted for 12 It is also worth noting that the values of Lij ultimately depend on the flux of DM-induced gamma rays, hence on the product σv · Ji (see Equation 3.1). Therefore, in order to compute Li and L (and its profile with respect to Ji) it is enough to know (in addition to J ) the values of Lij vs. σv , for fixed values of Ji (e.g. for J obs,i ), since: This is a particularly useful property, since it allows a significant reduction of the computing time requested for the profile of L over Ji, which can be explicitly written as: where the values of Li vs. Ji can be computed using Equations 3.7 and 3.8. In addition, this allows the combination of results from different instruments and targets, starting from tabulated values of Lij vs. σv , for a fixed value of Ji and profiled with respect to µij. These values can be produced and shared by different experiments without the need of releasing or sharing any of the internal information used to produce them, such as event lists or IRFs.
the sake of clarity): where g is the expected number of gamma rays detected with reconstructed energy E in the telescope range [E min , E max ] and an observation time T obs , i.e.: 14) τ and b (nuisance parameters) are the ratio of exposures between the OFF (background control) and ON (signal) regions and the expected number of background events in the OFF region, respectively. h and f are, respectively, the probability density functions (PDFs) for measured OFF and ON events with reconstructed energy E . h is obtained by fitting or interpolating the measured differential event rate from one or several additional background control regions observed simultaneously to the ON and OFF regions (see Ref. [53] for details). Finally, f can be written as: A(E) is the telescope effective area computed for a gamma-ray source with the morphology expected for Segue 1 according to Equations 3.1 to 3.3, after analysis cuts (including the selection of events in the ON region with radius of 0.122 • around Segue 1 center). G(E; E ) is the PDF for the energy estimator (E ) for a given true energy (E). Note that p is also a PDF and therefore does not depend on T obs , J or σv .

Results
We compute one-sided, 95% confidence level upper limits to σv by numerically solving the equation −2 ln λ P ( σv 2.71 | D) = 2.71, with σv restricted to the physical (≥0) region. This prescription is the one used by the Fermi -LAT Collaboration in Refs. [33,34], and differs slightly from the one used by MAGIC in Ref. [32]. This has consequences on the comparison of the results presented here with previous MAGIC results, which will be discussed below. Fermi -LAT individual limits are σv UL = 1.2 × 10 −24 cm 3 s −1 and 6.8 × 10 −25 cm 3 s −1 , respectively, whereas the combined analysis yields σv UL = 4.6×10 −25 cm 3 s −1 . Considering no uncertainties in J produces the combined limit σv UL = 3.4 × 10 −25 cm 3 s −1 . Figures 2 and 3 show our main results: the 95% confidence level limits to σv for DM particles with masses between 10 GeV and 100 TeV annihilating into SM pairs bb, W + W − , τ + τ − and µ + µ − , as obtained from the combination of Fermi -LAT and MAGIC observations of 15 dSphs and Segue 1 alone, respectively. We also show the observed limits obtained with MAGIC and Fermi -LAT data considered individually. The combined limits are compared to the median, and the 68% and 95% containment bands, expected under the null (H 0 : σv = 0) hypothesis. These are estimated from the distributions of limits obtained by applying our analysis to 300 independent H 0 realizations. Each realization consists of data from Fermi -LAT observations of one empty field per considered dSph, combined with fast simulations of MAGIC Segue 1 observations. The empty fields constituting the LAT realizations were selected by choosing random sky positions with |b| > 30 • centered at least 0.5 • away from a source in the 3FGL catalog. MAGIC fast simulations consist of a set of event energies randomly generated according to the background PDF, h (see Section 3.2 and Ref. [32] for details), for both ON and OFF regions. In all cases, we assume the same exposures on the different considered dSph as for the real data, and J-factors randomly selected according to the PDF in Equation 3.6.
We observe no evidence for DM annihilation in our data set (which would appear as a positive deviation of the limit from the null hypothesis). The maximum such deviation has a (local) significance of about 0.3 σ, found for m DM 10 GeV in the bb and τ + τ − annihilation channels, and for m DM 3 TeV in the µ + µ − and τ + τ − annihilation channels. We also observe a negative 2 σ fluctuation at m DM 5 TeV in the bb and W + W − annihilation channels and m DM 500 GeV in the µ + µ − and τ + τ − annihilation channels. A deviation of  (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 µ + µ − ), with a maximum improvement of the combined limits with respect to the individual ones by a factor ∼2. MAGIC individual limits shown in this work are stronger by up to a factor ∼4 than those presented in Ref. [32], which needs a dedicated explanation, provided in Appendix A.
Systematic uncertainties in the determination of the J-factors would weaken the limits on σv that can be inferred from the MAGIC and Fermi -LAT observations. In this work we take the J-factors and associated uncertainties from Ref. [42], which are largely consistent with the independent analysis of Ref. [43]. The analysis presented in Ref. [44], using a more flexible parameterization for the stellar velocity distributions and more stringent criteria for stellar membership, produces substantially larger J-factor uncertainties for the ultrafaint dSphs used in the present study (0.4-2.0 dex versus 0.2-0.3 dex). They find the J-factor of Segue 1 to be particularly uncertain due to the ambiguous membership status of a few stars, and when excluding stars with membership probability less than 95% they find log 10 (J obs /(GeV 2 cm −5 )) =17.0 +2.1 −2.2 . This value is more than two orders of magnitude smaller than the J-factor used in the present analysis (log 10 (J obs /(GeV 2 cm −5 )) =19.5±0.3) and would imply substantially weaker limits from the MAGIC observations of Segue 1. Although a full examination of the J-factor uncertainties is beyond the scope of the present work, we note that these uncertainties do not diminish the power of the joint likelihood approach and can be fully accounted for in our analysis procedure by updating the uncertainty parameter in the J-factor likelihood (Equation 3.6).

Discussion and Conclusions
This work presents, for the first time, limits on the DM annihilation cross-section from a comprehensive analysis of gamma-ray data with energies between 500 MeV and 10 TeV. Using a common analysis approach (both in the applied statistical methods and in the determination of the J-factors) we have combined the MAGIC observations of Segue 1 with Fermi-LAT observations of 15 dSphs. This allowed the computation of meaningful global DM limits, and the direct comparison of the individual results obtained with the different instruments. Our results span the DM particle mass range from 10 GeV to 100 TeV -the widest range covered by a single gamma-ray analysis to date.
We find no signal of DM in our data set. Consequently, we set upper limits to the annihilation cross-section. The obtained results are the most constraining bounds in the considered mass range, from observations of dSphs. For the low-mass range, our results (dominated by Fermi -LAT) are below the thermal relic cross-section σv 2.2 × 10 −26 cm 3 s −1 . In the intermediate mass range (from few hundred GeV to few tens TeV, depending on the considered annihilation channel), where Fermi -LAT and MAGIC achieve comparable sensitivities, the improvement of the combined result with respect to the individual ones reaches a factor ∼2. In addition, we present limits to DM particle masses above 10 TeV (dominated by MAGIC) that have not been shown before.
Our global analysis method is completely generic, and can be easily extended to include data from more targets, instruments and/or messenger particles provided they have similar sensitivity to the considered DM particle mass range. Of particular interest is the case of a global DM search from dSphs including data from all current gamma-ray (Fermi -LAT, MAGIC, H.E.S.S, VERITAS, HAWC) and neutrino (Antares, IceCube, SuperKamiokande) instruments, and we hereby propose a coordinated effort toward that end. Including results obtained from other types of observational targets like the Galactic Center, galaxy clusters or others is formally also possible, but a common approach to the J-factor determination remains an open question. In the future, this analysis could include new instruments like CTA [58], GAMMA-400 [59], DAMPE [60] or Km3Net [61]. 13 Our global approach offers the best chances for indirect DM discovery, or for setting the most stringent limits attainable by these kind of observations. 13 The combination with results from direct searches or accelerator experiments following a similar approach is, in principle, also formally possible, but it would necessarily be model-dependent. This possibility should be however regarded as the culminating step of our proposal. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

A Comparison with previous MAGIC results
The data, IRFs and likelihood functions (Equations 3.12 to 3.16) used in this work are the same as for previous MAGIC Segue 1 publication [32]. Aside from enlarging the explored DM mass range, the only differences between the two works are: the assumed J-factor central value, the treatment of the J-factor statistical uncertainties, and the prescription used for cases when the maximum of the profile likelihood is found in the non-physical ( σv < 0) region. We motivate each of these changes and comment on their effect in the following paragraphs.
• In this work, J-factors are obtained following Ref. [42] and assuming an NFW DM density profile, whereas previous MAGIC results were obtained assuming an Einasto profile [31,55]. This change has been introduced to homogenize MAGIC and Fermi -LAT computation of J-factors, and produces a factor 2 lower (stronger) MAGIC limits with respect to the previously published ones.
• Previous MAGIC results did not include statistical uncertainties in the J-factor. This was justified by the fact that σv UL scales with 1/J, and therefore the provided results allow to compute limits for any J-factor value. This argument is true only for single-target limits, but not for results obtained combining observations from different targets with different J values and uncertainties. In this study, we include the statistical uncertainties on J for all targets as described in Section 3.2. In consequence, MAGIC limits increase (degrade) by a factor between ∼1.4 and ∼1.8 (depending on the considered annihilation channel and mass) with respect to the previously published ones.
However, for the scalability argument given before, it is also interesting to provide single-target (Segue 1) results with and without uncertainties on the J-factor, which are shown in Figure 3.
• The requirement of producing upper limits in the physical ( σv ≥ 0) region normally leads to ad-hoc recipes implying over-coverage (see e.g. Ref. [52]). The problem arises when σv < 0, and it is aggravated when σv 2.71 < 0. One possible prescription to deal with these cases is to restrict −2 ln λ P ( σv | D) to σv ≥ 0 [56]. This is the solution adopted by Fermi -LAT in Refs. [33,34], and the one we have followed in this work.
Another, more conservative (i.e. producing larger over-coverage) prescription, was followed by MAGIC in Ref. [32], consisting in computing the 95% confidence limits as σv svt = σv 2.71 − σv , whenever σv < 0. σv svt corresponds to what the authors in Ref. [57] define as the "sensitivity" of the measurement. Using this prescription, σv svt is the lowest possible value of the upper limit, irrespective of the presence of arbitrarily intense negative background fluctuations.
The latter approach cannot be applied to Fermi -LAT analysis since, due to low bin statistics, the likelihood function can be undefined for negative σv values (see Equation 4 in Ref. [34]). We therefore homogenize the treatment of limits close to the physical boundary for our whole data set by limiting σv to non-negative values. The difference between old and new prescriptions in MAGIC individual limits goes, depending on the annihilation channel and DM particle mass, from none (when σv ≥ 0) up to a factor ∼4 (for ∼2σ background fluctuations, see Figures 2 and 3).