Mapping Extragalactic Dark Matter Annihilation with Galaxy Surveys:
A Systematic Study of Stacked Group Searches
Abstract
Dark matter in the halos surrounding galaxy groups and clusters can annihilate to highenergy photons. Recent advancements in the construction of galaxy group catalogs provide many thousands of potential extragalactic targets for dark matter. In this paper, we outline a procedure to infer the dark matter signal associated with a given galaxy group. Applying this procedure to a catalog of sources, one can create a fullsky map of the brightest extragalactic dark matter targets in the nearby Universe (), supplementing sources of dark matter annihilation from within the Local Group. As with searches for dark matter in dwarf galaxies, these extragalactic targets can be stacked together to enhance the signals associated with dark matter. We validate this procedure on mock Fermi gammaray data sets using a galaxy catalog constructed from the DarkSky body cosmological simulation and demonstrate that the limits are robust, at levels, to systematic uncertainties on halo mass and concentration. We also quantify other sources of systematic uncertainty arising from the analysis and modeling assumptions. Our results suggest that a stacking analysis using galaxy group catalogs provides a powerful opportunity to discover extragalactic dark matter and complements existing studies of Milky Way dwarf galaxies.
I Introduction
Dark matter (DM) annihilation into visible final states remains one of the most promising avenues for discovering nongravitational interactions in the dark sector. While an individual annihilation event is rare, the probability of observing it can be maximized by searching for excess photons in regions of high dark matter density. The center of the Milky Way is potentially one of the brightest regions of DM annihilation as seen from Earth, but the astrophysical uncertainties associated with the baryonic physics at the heart of our Galaxy motivate exploring other targets. Gammaray studies of DMdominated dwarf galaxies in the Local Group currently provide some of the most robust constraints on the annihilation cross section Albert et al. (2017); Ackermann et al. (2015a). However, many more potential targets are available beyond the Local Group. This paper proposes a new analysis strategy to search for DM emission from hundreds more DM halos identified in galaxy group catalogs.
A variety of methods have been used to study gammaray signatures of extragalactic DM annihilation, including modeling potential contributions to the Isotropic GammaRay Background Bengtsson et al. (1990); Bergstrom et al. (2001); Ullio et al. (2002); Bottino et al. (2004); Bertone et al. (2005); Bringmann and Weniger (2012); Ajello et al. (2015); Di Mauro and Donato (2015); Ackermann et al. (2015b); Feng et al. (2017), measuring the Fermi autocorrelation power spectrum Ackermann et al. (2012); Fornasa et al. (2016); Ando et al. (2007); Ando and Komatsu (2013), and crosscorrelating the Fermi data with galaxy counts Branchini et al. (2017); Xia et al. (2011); Ando (2014); Ando et al. (2014); Xia et al. (2015); Regis et al. (2015); Cuoco et al. (2015); Ando and Ishiwata (2016), cosmic shear Camera et al. (2015); Tröster et al. (2016); Choi et al. (2016); Camera et al. (2013); Shirasaki et al. (2015, 2014, 2016) and lensing of the Cosmic Microwave Background Fornengo et al. (2015); Feng et al. (2017). These methods typically rely on using a probabilistic distribution of the DM annihilation signal on the sky. Our approach is more deterministic in nature. In particular, we treat a collection of known galaxies as seeds for DM halos. The properties of each galaxy—such as its luminosity and redshift—enable one to deduce the characteristics of its associated halo and the expected DMinduced gammaray flux from that particular direction in the sky. In this way, we can build a map of the expected DM annihilation flux that traces the observed distribution of galaxy groups.
In certain ways, our approach resembles that used in previous studies of DM annihilation from individual galaxy clusters. For example, most recently the Andromeda galaxy Ackermann et al. (2017) and Virgo cluster Ackermann et al. (2015c) have been the subject of dedicated study by the Fermi Collaboration. Other work has inferred the properties of the DM halos associated with galaxy clusters detected in Xrays Ackermann et al. (2010); Ando and Nagai (2012); Ackermann et al. (2014a); Anderson et al. (2016); Ackermann et al. (2016); Ahnen et al. (2016); Liang et al. (2016). Most of these studies focused on a small number of galaxy clusters and obtained DM sensitivities weaker than those from dwarf galaxies.
Recent advancements in the development of galaxy group catalogs allow us to now build a fullsky map of the nearby galaxies that should be the brightest DM gammaray emitters. Catalogs based primarily on the 2MASS Redshift Survey (2MRS) Huchra et al. (2012) provide an unprecedented amount of information regarding a group’s constituents and halo properties Tully (2015); Kourkchi and Tully (2017); Lu et al. (2016). This information allows us to build a list of the brightest extragalactic DM targets on the sky and to perform a stacked analysis for gammaray emission from their centers. A gammaray line search using this methodology was recently performed by Ref. Adams et al. (2016). Our focus is on continuum DM signatures, which carry considerably more complications in terms of the treatment of astrophysical backgrounds.
In a companion paper Lisanti et al. (2017), we present results implementing a stacked analysis of the group catalogs from Ref. Tully (2015); Kourkchi and Tully (2017) on Fermi data and show explicitly that this method yields competitive sensitivity to the dwarf searches. Here, we present the full details of the analysis method and a thorough discussion of the systematic uncertainties involved in deducing the DMinduced flux associated with a given galaxy group. To fully understand these uncertainties, we apply these methods on mock data where it is possible to compare the inferred DM properties to their true values. For this purpose, we use the DarkSky cosmological body simulation Skillman et al. (2014); Lehmann et al. (2017) and an associated galaxy catalog from Ref. Lehmann et al. (2017). We emphasize that, while we illustrate the analysis method on gammaray data, it can also be applied to other wavelengths and even other messengers, such as neutrinos.
This paper is organized as follows. In Sec. II, we describe how to build DM annihilation flux maps starting from a galaxy group catalog and discuss the associated systematic uncertainties. Sec. III presents a detailed description of the statistical methods that we follow to implement the stacking. We show the results of applying the limitsetting and signal recovery procedures on mock data in Sec. IV and conclude in Sec. V. Appendix A provides a detailed discussion of the factor expressions used in the main text. Appendix B discusses the validity of approximations made in the profile likelihood procedure. Appendix C shows the results of performing a fullsky template study, as a contrast to the stacking results presented in the main text.
Ii Tracing Dark Matter Flux with Galaxy Surveys
In this Section, we describe how to construct catalogs of extragalactic DM targets starting from a list of galaxy groups. We begin by reviewing the properties of the galaxy group catalogs and then describe how to predict the DM signal from a given galaxy group and quantify the systematic uncertainties of this extrapolation.
ii.1 Galaxy and Halo Catalogs
The approach that we use throughout this work relies on galaxy surveys as an input. Different galaxy catalogs span a range of redshifts and luminosities. Optimal catalogs for DM searches should cover as much of the sky as possible (to increase statistics) and sample low redshifts. The strength of the DM signal increases at lower redshifts due to accretion of mass at late times, affecting both the halo mass distribution and substructure Ando (2014). In contrast, the integrated gammaray flux of standard astrophysical sources, such as Active Galactic Nuclei and starforming galaxies, is expected to peak at higher redshifts between 0.1 and 2 depending on the specific source class and model for its unresolved contribution Ando (2014); Xia et al. (2015).
The Two Micron AllSky Survey Extended Sources Catalog (2MASS XSC) Bilicki et al. (2013); Huchra et al. (2012) satisfies the criteria listed above and has been used extensively in past crosscorrelation studies Ando et al. (2014); Ando (2014); Ando and Ishiwata (2016); Cuoco et al. (2015); Regis et al. (2015); Xia et al. (2011, 2015). The XSC is an allsky infrared survey that consists of approximately one million galaxies up to a limiting magnitude of mag. Several redshift surveys based on the 2MASS XSCmap the redshifts associated with these galaxies. The 2MRS Huchra et al. (2012), for example, samples about 45,000 galaxies in the 2MASS XSC with redshifts to a limiting magnitude of mag. This corresponds to a nearly complete galaxy sample up to redshifts of , which is ideal for DM studies.
Galaxies from large surveys such as 2MASS can be organized into group catalogs. A group of gravitationallybound galaxies shares a DM host halo. The brightest galaxy in the group is referred to as the central galaxy; the additional galaxies are bound satellites surrounded by their own subhalos. The total luminosity of the galaxies in the group is a good predictor of the mass of the DM host halo. A variety of group finders have been developed and applied to the 2MASS data set Tully (2015); Lu et al. (2016); Kourkchi and Tully (2017), using the 2MRS which adds information in the redshift dimension. The groups in these catalogs range from cluster scales with 190 members and associated halo masses of 10 M, down to much smaller systems with only a single member. Galaxy group catalogs are especially relevant for the present study, since halo properties tend to be correlated with properties of galaxy groups rather than those of individual galaxies.
While in our companion paper Lisanti et al. (2017) we use information from the 2MASS group catalogs in the analysis of Fermi data, we focus on a catalog of simulated galaxies and halos here. We use the DarkSky400 cosmological body simulation (version ds14_i) Skillman et al. (2014); Lehmann et al. (2017) and an associated band galaxy catalog. Using the code 2hot Warren (2013), DarkSky400 follows the evolution of particles of mass M in a box 400 Mpc per side. Initial perturbations are tracked from to today, assuming . The halo catalog was generated using the Rockstar halo finder Behroozi et al. (2013); Lehmann et al. (2017). Crucially, the simulation covers the relevant redshift space for DM studies.^{1}^{1}1The snapshot of the simulation analyzed in this work is taken at , but we will refer to distance using redshift because that is the more appropriate language when applied to real data. In particular, an observer at the center of the simulation box has a complete sample of galaxies out to , with the furthest galaxies extending out to . In our work, we only consider groups located within , which is the approximate redshift cutoff of the catalogs in Ref. Tully (2015); Kourkchi and Tully (2017); Lu et al. (2016). We include only wellresolved halos in our analysis by imposing a lower cutoff of M on the mass of included host halos. The associated galaxy catalog is generated using the abundance matching technique following Ref. Behroozi et al. (2010); Reddick et al. (2013) with luminosity function and twopoint correlation measurements from the Sloan Digital Sky Survey (SDSS). Specifically, the model from Ref. Lehmann et al. (2017) is used, which was shown to provide the best fit to SDSS twopoint clustering. The DarkSky galaxy catalog contains the same information that would be found in, e.g. the 2MASS galaxy catalog and associated group catalogs, such as individual galaxy luminosities and sky locations.
Figure 1 shows a sky map of the galaxy counts in DarkSky up to for an observer at the center of the simulation box. It is a HEALPix Gorski et al. (2005) map with resolution nside=128. To first approximation, the galaxies are isotropically distributed throughout the sky. However, regions of higher and lower galaxy density are clearly visible. Note that this is shown for a particular sky realization and placing the observer in different parts of the DarkSky box would change the regions of contrasting galaxy density.
ii.2 Dark Matter Annihilation Flux Map
One can predict the DM annihilation flux associated with a halo that surrounds a given galaxy group. This requires knowing the halo’s properties, including its mass and concentration. In this subsection, we discuss how to determine the flux when the halo’s properties are known exactly. Then, in the following subsection, we consider how to generalize the results to the more realistic scenario where the halo properties have to be inferred.
Each halo in DarkSky is fit by the Rockstar halo finder with a NavarroFrenkWhite (NFW) distribution Navarro et al. (1996) of the form
(1) 
where is the scale radius and is the normalization. The NFW parameters are determined from the parameters that are provided for each DM halo—specifically, its redshift , virial mass , virial radius , and virial concentration parameter .
In the simplest scenarios, the annihilation flux factorizes as
(2) 
where is the photon energy and () encodes the particle physics (astrophysical) dependence. The particle physics contribution is given by
(3) 
where is the DM mass, is its annihilation cross section, is its branching fraction to the annihilation channel, is the photon energy distribution in this channel, which is modeled using PPPC4DMID Cirelli et al. (2011), and is the redshift. We consider the case of annihilation into the channel as a generic example of a continuum spectrum. Of course, the exact limits will vary for different spectra, and one should consider a range of final states when applying the method to data, or use model independentapproaches (see, e.g., Ref. Elor et al. (2015, 2016)).
The factor is defined as the integral along the lineofsight of the squared DM density of the observed object:^{2}^{2}2As defined, the factor has units of []. This definition is convenient for extragalactic objects, but beware because another common definition of the factor involves dividing out by a solid angle to remove the units of . A detailed discussion of the units is provided in Appendix A.
(4) 
where is the lineofsight distance and is the socalled boost factor. The boost factor accounts for the enhancement in the flux due to the annihilation in subhalos. For the case of extragalactic objects, one can obtain a closed form solution that is an excellent approximation to the integral in Eq. 4, which is proportional to
(5) 
where is the comoving distance (a function of redshift, ), is the critical density, and is the concentration. In our analysis, we calculate the factor exactly, but the scaling illustrated in Eq. 5 is useful for understanding the dependence of on the halo mass and concentration. The derivation of the factor expression is reviewed in detail in Appendix A, where we also show the result for the Burkert profile.
Figure 1 illustrates the truth factor map associated with DarkSky, obtained by putting the observer in the center of the simulation box. This map is constructed by applying Eq. 4 to all host halos in the DarkSky catalog and using the boost model from Ref. Bartels and Ando (2015) to describe the contribution from substructure. Once the factors are known, the expected photon counts per pixel can be determined using Eq. 2 and Fermi’s exposure map. This is also shown in Fig. 1, assuming a DM particle with GeV that annihilates to with cms. Not all the pixels that contain one or more galaxies correspond to significant regions of DM annihilation. The DM annihilation flux is largest for the most massive, concentrated, and/or closest galaxy groups.
Note that when constructing Fig. 1, we perform the angular integrals in Eq. 4 as a function of angular extent, . In doing so, we implicitly assume that the boost factor is simply a multiplicative factor. In reality, the boost factor likely broadens the angular profile, because the subhalo annihilation should extend further away from the halo center. However, since the angular extent of the annihilation in most halos is small compared to the instrument pointspread function (PSF), we do not model this extension here. Some nearby halos may have significantly larger angular extent, as would be expected for the Andromeda galaxy. Nevertheless, such considerations need to be made case by case and are discussed in detail in our companion paper Lisanti et al. (2017), where we choose to exclude Andromeda due to its size.
Figure 2 is a heatmap representing the average factor, for a given and , of the DarkSky halos in the above configuration. The halos span a wide range of masses and redshifts, with factors averaging over several orders of magnitude from 10. The largest factors are observed for the most massive, clustersized halos at 0.01–0.02, as well as for lessmassive halos at smaller redshifts ().
ii.3 Uncertainties in Halo Modeling
Now, we consider more carefully the systematic uncertainties associated with modeling the halo properties. A halo with an NFW density profile has a factor dictated by its parameters as given in Eq. 5. In addition to the distance, the factor also depends on the virial mass and concentration.^{3}^{3}3Note that uncertainties on the halo redshift also feed into the factor. However, we consider this uncertainty to be subdominant for spectroscopically determined redshifts. For nearby halos, where the relation between distance and redshift is nontrivial, the uncertainty on the distance can be noticeably larger, and as high as 5% Tully et al. (2016). Nonetheless, even such uncertainties are considerably smaller than those associated with the mass and concentration, and so we do not consider them. Therefore, any uncertainty in the determination of these halo properties is propagated through to the uncertainty on the DM annihilation flux. Up until now, we have taken the halo mass and concentration directly from DarkSky, but in practice these parameters need to be inferred from properties of the observed galaxy groups.
Within DarkSky, the halo mass can be inferred from the absolute luminosity of its associated galaxy group. We obtain a deterministic relation following a procedure similar to that in Ref. Vale and Ostriker (2006), which derived a relation between the band galaxy luminosity and the mass of its DM halo. The left panel of Fig. 3 shows the true masses for the DarkSky halos, as a function of central galaxy luminosity (green) or the total luminosity, which includes the luminosity of the satellite galaxies (red). The DarkSky catalog provides the associations for all galaxies, central and satellite, so we include all satellites that are associated to the group when calculating the total absolute luminosity. This is similar to what is done in published group catalogs Tully (2015); Kourkchi and Tully (2017); Lu et al. (2016), where they account for the loss in luminosity of satellite galaxies that are farther away.
From Fig. 3, we see that the spread in the associated halo mass increases above L, up to the brightest galaxy at L, when the central galaxy luminosity is used. In contrast, the spread is significantly smaller when the total luminosity is used, making it a better predictor for the halo mass. As demonstrated in the right panel of Fig. 3, including the satellite luminosities allows one to better reconstruct the halo mass. Therefore, we use the median relation thus obtained as our fiducial case to infer the central mass estimate, and we use the spread in the relation to infer the uncertainty on the mass. Note that the relation shown in Fig. 3 is constructed by binning the DarkSky data in luminosity and calculating the 16, 50, and 84 percentiles in ; different results would be obtained by binning in and then constructing the percentiles from the luminosity distributions. This procedure is similar to that adopted by galaxy group catalogs to infer the halo mass Tully (2015); Lu et al. (2016); Kourkchi and Tully (2017). Using this relation, we can infer the halo mass and uncertainty for each galaxygroup host halo in DarkSky.
DM halos of the same mass can have very different characteristics, usually reflecting their distinct formation history and environment. One such characteristic is the halo’s virial concentration . The scale radius is the relevant quantity to compare to as it indicates an isothermal slope for the density profile, which is required for a flat rotation curve. The virial radius corresponds to the spherical volume within which the mean density is times the critical density of the Universe at that redshift. We use with in accordance with Ref. Bryan and Norman (1998). The cosmology associated with the DarkSky simulation is used throughout, with , and .
In general, the concentration correlates strongly with halo mass due to the dependence of halo formation time on mass—on average, lower mass halos tend to be more concentrated because they collapsed earlier, when the Universe was denser. For the same reason, the concentration is sensitive to the cosmology, which determines how early halos start to assemble. The concentration of field halos has been extensively studied and several concentrationmass relations have been proposed in the literature, usually based on body simulations or physically motivated analytic approaches Duffy et al. (2008); Prada et al. (2012); SánchezConde and Prada (2014); Diemer and Kravtsov (2015); Correa et al. (2015); Moliné et al. (2017); Klypin et al. (2016); Dutton and Macciò (2014). In the left panel of Fig. 4, we show the median value of the concentrationmass relation derived directly from the DarkSky simulation, as well as the middle 68 and 95% spread. The middle 68% scatter in the relation is typically in the range 0.140.19 across the halo mass range considered. For comparison, we also show several concentration models that are commonly used in the literature. As is standard in the literature Hütten et al. (2016); SánchezConde and Prada (2014), we model the uncertainty in the concentration, for a given virial mass, as a lognormal distribution around its median value.
To summarize, it is possible to infer the halo mass from the luminosity of the galaxy group and to then obtain the concentration. The final remaining property that is needed to solve for the factor in Eq. 4 is the boost factor, which depends on the distribution and minimum cutoff of the subhalos’ mass. The boost factor encapsulates the complicated dependence of the subhalo mass distribution on both the particle physics assumptions of the DM model as well as the dynamics of the host halo formation. A variety of different boost models typically used in the literature are illustrated in the right panel of Fig. 4. As our fiducial case, we adopt the boost model of Ref. Bartels and Ando (2015) (labeled as ‘Bartels Boost Model’), which selfconsistently accounts for the concentrationmass relation of subhalos (compared to field halos) as well as the effects of tidal stripping. Specifically, in the subhalo mass function , we use a minimum subhalo mass cutoff of M and slope that varies selfconsistently with host halo mass while accounting for evolution effects (see Ref. Bartels and Ando (2015) for details).
We have now built up a framework that allows us to determine the expected DM annihilation flux map associated with a catalog of galaxy groups. Next, we show how to use this information to search for signals of DM from hundreds of galaxy groups.
Iii Statistical Methods
In this work, we introduce and study a statistical procedure to search for gammaray signals from DM by stacking galaxy groups. All analyses discussed here are run on mock data, which is based on the expected astrophysical contributions to the real Fermi data set. When building this mock data set, we include contributions from (1) the diffuse emission, for which we use the Fermi Collaboration’s p7v6 model; (2) isotropic emission; (3) emission from the Fermi Bubbles Su et al. (2010); and (4) emission from point sources in the Fermi 3FGL catalog Acero et al. (2015). The overall flux normalization for each component must be known a priori to create the mock data. To obtain this, we fit spatial maps of (1)–(4) above to the actual Fermi data. We use 413 weeks of UltracleanVeto (all PSF quartile) Pass 8 data collected between August 4, 2008 and July 7, 2016. We break the data into 40 equally logspaced energy bins between 200 MeV and 2 TeV, applying the recommended data cuts: zenith angle , DATA_QUAL , and LAT_CONFIG . To minimize the Galactic contamination in this initial fit, we mask the region as well as the 68% containment radius for the 300 brightest and most variable sources in the 3FGL catalog. We emphasize that these masks are only used when creating the mock data and not in the stacked analysis. The fitting procedure described here provides the expected astrophysical background contribution from the real data. Monte Carlo (MC) is then generated by summing up these contributions and taking a Poisson draw from the resulting map. In the following discussion, we will show how results vary over different MC realizations of the mock data as a demonstration of Poisson fluctuations in the photon distribution.
We now describe in detail the statistical procedure we employ to implement the stacking analysis on the mock data. We perform a templatefitting profile likelihood analysis in a 10 regionofinterest (ROI) around each group. Template studies amount to describing the sky by a series of spatial maps (called templates). The normalization of each template is proportional to its relative gammaray flux. We use five templates in our study. The first four are associated with the known astrophysical sources (1)–(4) described above. Within of the halo center, we independently float the normalization of each 3FGL source.^{4}^{4}4The results do not change when floating all the point sources together as one combined template. This can potentially cause problems when implemented on data, however, because the 3FGL normalizations can be erroneous in certain energy bins. Allowing the normalizations of the sources to float separately helps to mitigate this potential problem. Sources outside this region may potentially contribute within the ROI because of the tails of the Fermi PSF. Therefore, between 10 and 18 of the halo center, we float the sources as a single template. The fifth and final template that we include is associated with the expected DM annihilation flux for the halo, which is effectively a map of the factor and is described in Sec. II. Note that all templates have been carefully smoothed using the Fermi PSF. The diffuse model is smoothed with the Fermi Science Tools, whereas other templates are smoothed according to the instrument response function using custom routines. Mismodeling the smoothing of either the point sources or individual halos can potentially impact the results.
A given mock data set, , is divided into 40 logspaced energy bins indexed by . Each energy bin is then spatially binned using HEALPix Gorski et al. (2005) with nside=128 and individual pixels indexed by . In this way, the full data set is reduced to a twodimensional array of integers describing the number of photons in energy bin and pixel . For a given halo, indexed by , only a subset of all the pixels in its vicinity are relevant. In particular, the relevant pixels are those with centers within 10 of the object. Restricting to these pixels leaves a subset of the data, which we denote by . Template fitting dictates that this data is described with a set of spatial templates binned in the same way as the data, which we label as , where indexes the different templates considered. The number of counts in a given pixel, energy bin, and region consists of a combination of these templates:
(6) 
Here, represents the set of model parameters. For Poissonian template fitting, these are given by the normalizations of the templates , i.e., . Note that the template normalizations have an energy but not a spatial index, as the templates have an independent degree of freedom in each energy bin as written, but the spatial distribution of the model is fixed by the shapes of the templates themselves. In principle, we could also remove this freedom in the relative emission across energy bins, because we have models for the spectra of the various background components, and in particular DM. Nevertheless, we still allow the template normalizations to float independently in each energy bin for the various backgrounds. This is more conservative than assuming a model for the background spectra, and in particular we can use the shape of the derived spectra as a check that the dominant background components are being correctly modeled. The spectral shape of the DM forms part of our model prediction, however, and once we pick a final state such as annihilation to two quarks, we fix the relative emission between the energy bins.
As we assume that the data comes from a Poisson draw of the model, the appropriate likelihood in energy bin and ROI is
(7) 
Of the templates that enter this likelihood, there are some we are more interested in than others. In particular, we care about the the DM annihilation intensity, which we denote as . We treat the normalizations of the templates associated with the known astrophysical emission as nuisance parameters, . Below, we will describe how to remove the nuisance parameters to reduce Eq. 7 to a likelihood profile that depends only on the DM annihilation intensity, but for now we have .
Importantly, the nuisance parameters have different values between ROIs, but the DM parameters do not. This is because the DM parameters, such as the DM mass, annihilation rate, and set of final states, are universal, while the parameters that describe the astrophysical emission can vary from region to region. We do, however, profile over the factor uncertainty in each ROI. Explicitly, each halo is given a model parameter , which is described by a lognormal distribution around the central value with width , both of which depend on the object and hence ROI considered. The factor error, , is determined by propagating the errors associated with the mass and concentration of a given halo. To account for this, we append the following addition onto our likelihood as follows:
(8)  
A detailed justification for this choice of the likelihood is provided in Appendix A.4. Note that this procedure does not account for any systematic uncertainties that can bias the determination of the factor.
The nuisance parameter can now be eliminated via the profile likelihood—see Ref. Rolke et al. (2005) for a review. Unlike for the other nuisance parameters, the value of does not depend on energy and so we eliminate the energydependent parameters first:
(9) 
The full implementation of the profile likelihood method as suggested by this equation requires determining the maximum likelihood for the template coefficients, for every value of . Nevertheless, an excellent approximation to the profile likelihood, which is computationally more tractable, is simply to set the nuisance parameters to their maximum value obtained in an initial scan where all templates are floated.^{5}^{5}5The DM template is only included for energy bins above 1 GeV. At lower energies, the large Fermi PSF leads to confusion between the DM, isotropic and point source templates, which can introduce a spurious preference for the DM template. We follow this procedure throughout in calculating likelihoods and discuss its validity in Appendix B.
Using this approach to determine the likelihood in Eq. 9, we can build a total likelihood by combining the energy bins. Once this is done, the likelihood depends on the full set of DM intensities , which are specified by a DM model , cross section , mass , and factor via Eq. 2. Explicitly:
(10) 
and recall that unlike the other parameters on the left hand side, the factor not only determines the , but also enters the likelihood through the expression in Eq. 8. We emphasize that in this equation, the DM model and mass specify the spectra, and thereby the relative weightings of the , whereas the cross section and factor set the overall scale of the emission.
The remaining step to get the complete likelihood for a given halo is to remove , again using profilelikelihood:
(11) 
This provides the full likelihood for this object as a function of the DM model parameters. The likelihood for the full stacked catalog is then simply a product over the individual likelihoods:
(12) 
Using this likelihood, we define a test statistic (TS) profile as follows:
(13)  
where is the cross section that maximizes the likelihood for that DM model and mass. From here, we can use this TS, which is always nonpositive by definition, to set a threshold for limits on the crosssection. When searching for evidence for a signal, we use an alternate definition of the test statistic defined as
(14) 
Iv Analysis Results
In this Section, we present the results of our analysis on mock data using the DarkSky galaxy catalog. We begin by describing the sensitivity estimates associated with this study, commenting on the impact of statistical as well as systematic uncertainties and studying the effect of stacking a progressively larger number of halos. Then, we justify the halo selection criteria that are used by showing that we can recover injected signals on mock data.
iv.1 Halo Selection and Limits
We now discuss the results obtained by applying the halo inference pipeline described in Sec. II and the statistical analysis described in Sec. III to mock gammaray data. We focus on the top 1000 galaxy groups in the DarkSky catalog, as ranked by the inferred factors of their associated halos, placing ourselves at the center of the simulation box. In addition, we mask regions of the sky associated with seven largescale structures that are challenging to model accurately: the Large and Small Magellanic Clouds, the Orion molecular clouds, the galaxy NGC5090, the blazar 3C454.3, and the pulsars J1836+5925 and Geminga.
While we start from an initial list of 1000 galaxy groups, we do not include all of them in the stacking procedure. A galaxy group is excluded if:

it is located within ;

it is located less than from the center of another brighter group in the catalog;

it has TS and ,
where is the bestfit cross section at any mass and is the bestfit limit set by any halo at the specified DM mass. Note that the second requirement is applied sequentially to the ranked list of halos, ordered by factor. We now explain the motivation for each of these requirements separately. The first requirement listed above removes groups that are located close to the Galactic plane to reduce contamination from regions of high diffuse emission and the associated uncertainties in modeling these. The second requirement demands that the halos be reasonably wellseparated, which avoids issues having to do with overlapping halos and accounting for multiple DM parameters in the same ROI. The nonoverlap criterion of 2 is chosen based on the Fermi PSF containment in the lowest energy bins used and on the largest spatial extent of gammaray emission associated with the extended halos, which collectively drive the possible overlap between nearby halos.
The final requirement excludes a galaxy group if it has an excess of at least 3 significance associated with the DM template that is simultaneously excluded by the other galaxy groups in the sample. This selection is necessary because we expect that some galaxy groups will have true cosmicrayinduced gammaray emission from conventional astrophysics in the real data, unrelated to DM. To identify these groups, we take advantage of the fact that we are starting from a large population of halos that are all expected to be bright DM sources in the presence of a signal. Thus, if one halo sets a strong limit on the annihilation rate and another halo, at the same time, has a large excess that is severely in conflict with the limit, then most likely the large excess is not due to DM. The worry here is that we could have misconstructed the factor of the halo that gave the strong limit, so that the real limit is not as strong as we think it is. However, with the TS and criteria outlined above, this does not appear to be the case. In particular, we find that the criteria very rarely rejects halos due to statistical fluctuations. For example, over 50 MC iterations of the mock data, halos (out of 1000) remain after applying the TS and cross section cuts alone, and the excluded halos tend to have lower factors, since there the requirement is more readily satisfied.
We expect that this selection criteria will be very important on real data, however, where real excesses can abound. In addition, as we will describe in the next subsection, injected signals are not excluded when the analysis pipeline is run on mock data. In an ideal scenario, we would attempt to understand the origin of these excesses by correlating their emission to known astrophysics either individually or statistically. In the present analysis, however, we take the conservative approach of removing halos that are robustly inconsistent with a DM signal and leave a deeper understanding of the underlying astrophysics to future work.
We apply the procedure outlined in Secs. II and III to the mock data to infer the 95% confidence limit on the DM annihilation cross section. The resulting sensitivity is shown by the blue dashed line in the left panel of Fig. 5, which uses the boost factor from Ref. Bartels and Ando (2015). For comparison, we also show the limit assuming no boost factor (red dashed line); note that the boostfactor model that we use provides a modest improvement to the limit. Because the limit can potentially vary over different MC realizations of the mock data, we repeat the procedure for 100 MCs (associated with different Poisson realizations of the map); the blue band indicates the middle 68% spread in the limit associated with this statistical variation.
To see how the limit depends on the observer’s location within the DarkSky simulation box, we repeat the procedure described above over nine different locations.^{6}^{6}6The nine locations we used are at the following coordinates Mpc/h in DarkSky: , , , , , , , , . The first listed location is our default position, and any time we use more than one location they are selected in order from this list. At each location, we perform 20 MCs and obtain the median DM limit. The green band in the left panel of Fig. 5 denotes the middle 68% spread on the median bounds for each of the different sky locations. In general, we find that the results obtained by an observer at the center of the DarkSky box are fairly representative, compared to random locations. Note, however, that this bound does not necessarily reflect the sensitivity reach one would expect to get with actual Fermi data. The reason for this is that the locations probed in DarkSky do not resemble that of the Local Group in detail. We will come back to this point below, when we compare the factors of the DarkSky halos to those from galaxy catalogs that map the local Universe.
The orange line in the left panel of Fig. 5 shows the limit obtained by requiring that the DM emission from the groups not overproduce the measured isotropic gammaray component Ackermann et al. (2015d). This should not be compared to the published DM bounds obtained with the Fermi Isotropic GammaRay Background Ackermann et al. (2015b) because that study accounts for the integrated effect of the DM annihilation flux from halos much deeper than those we consider here. The inclusion of these halos results in a total flux that can be greater than those from our sample by over an order of magnitude. Nevertheless, this gives an idea of how much we gain by resolving the spatial structure of the local DM population and knowing the locations of the individual galaxy groups.
The right panel of Fig. 5 shows the effect of propagating uncertainties associated with inferring the halo properties. The green line indicates how the limit improves when no uncertainties are assumed, i.e., we can perfectly reconstruct the virial mass and concentration of the halos. The sensitivity reach improves by roughly a factor of two in this case. We further show the effect of individually reducing the error on (dashed purple line) and (purple line) by 50%. The reductions in the uncertainties provide only marginal improvements to the overall sensitivity, still far below the level of systematic uncertainty associated with extragalactic analyses in general.
It is interesting to study how the limit scales with the number of halos, , included in the stacking procedure. This result is shown in Fig. 6 for and GeV, for four different observer locations in the simulation box. The dashed red line indicates the median 95% confidence limit. The red bands are the 2.5, 16, 84 and 97.5 percentiles on the limit, obtained from 100 MC realizations of the mock data. We observe that the limit typically improves continuously for the first 10 halos. As more halos are included in the stacking, the gains diminish. For some sky locations, the limit simply remains flat; for others we see some marginal improvements in the bounds. These results are consistent, within uncertainties, between the DM masses and the different sky locations of the observer.
We emphasize that the scaling on can be very different on applicaton to real data, because the distribution of factors in the random DarkSky locations is not representative of our own environment in the Local Group and also some halos can have residuals that are not related to DM but rather to mismodeling or real cosmicray–induced emission from the galaxy groups. The former point is demonstrated in Fig. 7, where we histogram the top 1000 factors associated with the baseline DarkSky analysis (blue line/band). For comparison, we also show the distributions corresponding to 2MRS galaxy group catalogs, specifically the Tully et al. Tully (2015); Kourkchi and Tully (2017) (green line) and the Lu et al. Lu et al. (2016) (red line) catalogs. We see that the distribution of factors for the 2MRS catalogs is skewed towards higher values compared to that from DarkSky. (Note that the cutoff at low factors is artificial and is simply a result of including 1000 halos for each catalog.)
The differences in the factor distributions can be traced to the redshift distribution of the galaxy groups, as illustrated in Fig. 8. We see specifically that the mass function of the top 1000 DarkSky halos in each of the random sky locations sampled is roughly consistent with that observed in the 2MRS catalogs. In contrast, the actual catalogs have more groups at lower than observed in the random DarkSky locations.
While a random location in the DarkSky box does not resemble our own Local Group, we can try to find specific locations in the simulation box that do. Therefore, we place the observer at ten random Milky Way–like halos in the simulation box, which have a mass M. More specifically, we select halos with mass and at least 100 Mpc from the box boundaries. The distribution of the top 1000 factors is indicated by the orange line/band in Fig. 7, while the corresponding mass and redshift distributions are shown in Fig. 8. We see that the redshift—and, consequently, factor—distributions approach the observations, though the correspondence is still not exact. A more thorough study could be done assessing the likelihood that an observer in DarkSky is located at a position that closely resembles the Local Group. However, as our primary goal here is to outline an analysis procedure that we can apply to actual data, we simply conclude that our own local Universe appears to be a richer environment compared to a random location within the DarkSky simulation box, which bodes well for studying the actual Fermi data.
iv.2 Signal Recovery Tests
It is critical that the halo selection criteria described in the previous section do not exclude a potential DM signal if one were present. To verify this, we have conducted extensive tests where we inject a signal into the mock data, pass it through the analysis pipeline and test our ability to accurately recover its cross section in the presence of the selection cuts. Figure 9 summarizes the results of the signal injection tests for two different observer locations in the DarkSky simulation box (top and bottom rows, respectively). We inject a signal in the mock data that is associated with annihilation for three different masses ( GeV) that traces the DM annihilation flux map associated with DarkSky. The dashed line in each panel delineates where the injected cross section, , matches the recovered cross section, .
The green line shows the 95% onesided limit on the cross section found using Eq. 13, with a TS threshold corresponding to . The green band shows the 68% containment region on this limit, constructed from twenty different MC realizations of the mock data set. Importantly, the limit on roughly follows—but is slightly weaker than—the injected signal, up until the maximum sensitivity is reached and smaller cross sections can no longer be probed. This behavior is generally consistent between the three DM masses tested and both sky locations. We clearly see that the limit obtained by the statistical procedure never excludes an injected signal over the entire cross section range.
Next, we consider the recovered cross section that is associated with the maximum test statistic, TS, in the total likelihood. The blue line in each panel of Fig. 9 shows the median value of over 20 MCs of the mock data. The blue band spans the median cross sections associated with . The inset plots show the median and 68% containment region for TS as a function of the injected cross section. The maximum test statistic is an indicator for the significance of the DM model and as such the distributions are only influenced by the data at high injected cross sections where TS has begun to increase. At lower injected cross sections, the distributions for are not meaningful.
Two issues are visible in Fig. 9: (i) at high injected cross sections, the bestfit recovered cross sections are systematically around 1 too high, and (ii) at high and low DM masses and nearzero injected cross sections, the distribution of TS deviates from the chisquare distribution. The first issue stems from the way we model the factor contribution to the likelihood, while the second arises from the approximations we make to perform the profile likelihood in a computationally efficient manner. Appendix B discusses these issues and ways they may be mitigated, in more detail.
V Conclusions
In this paper, we introduced a procedure to build a fullsky map of extragalactic DM targets based on galaxy surveys and demonstrated this methodology using the DarkSky cosmological simulation. Starting from the galaxies in the DarkSky catalog, we inferred the properties of their respective DM halos using the galaxyhalo connection. In so doing, we identified the halos that are the brightest sources of extragalactic DM annihilation and which act as the best annihilation targets. This procedure allows us to account for the fact that not all galaxy groups are expected to be bright DM emitters; the most massive, concentrated, and/or most nearby galaxies dominate the signals. By building a map of extragalactic DM targets, we can focus our search for DM annihilation on the most relevant regions of sky. This philosophy contrasts with that of crosscorrelation studies, which treat all galaxies as equally good targets for DM.
With a list of extragalactic DM halos in hand, as well as their inferred factors, we performed a stacked analysis to search for gammaray signatures of DM annihilation in mock data. We described the likelihood procedure for the stacking analysis in detail. There are two clear advantages to this approach over, say, a fullsky template study.^{7}^{7}7 Appendix C includes a more detailed discussion of using fullsky DM annihilation flux templates. First, focusing on smaller regions around each halo significantly reduces the sensitivity to mismodeling of the foregrounds. Second, uncertainties on the predicted DM annihilation flux can be straightforwardly included in the likelihood function. In particular, we outlined how uncertainties in the factors, which arise from the determination of the virial mass and concentration, are marginalized over in the analysis.
We presented limits on the DM annihilation cross section for mock data and, most importantly, demonstrated that the analysis procedure robustly recovers injected signals. We found that the sensitivity improves by nearly two orders of magnitude when the structure of extragalactic DM emission on the sky is accounted for, rather than simply assuming an isotropic distribution. Typically, the limit is dominated by the brightest halos in the stacking, though this varies depending on the location in the simulation box. The factor distribution of nearby groups in our own Galaxy differs from the random locations sampled in the DarkSky box, which can change the number of halos that dominate the limit. In actuality, one would want to continue adding halos to the analysis—ranked starting from the brightest factors—until the gains in the limit are observed to level off.
One advantage of using the DarkSky simulation in this initial study is that the truth information for all the halos is known. We can therefore study how the DM limits improve when the virial mass and concentration of the halos are known precisely. For this ideal scenario, we find that that the limits improve by roughly 50% over those obtained by marginalizing over uncertainties. This suggests that a concrete way to improve the bounds on DM annihilation is to reduce the uncertainties on and for the brightest halos in the catalog.
The substructure boost factor remains one of the most difficult systematics to handle. In this work, we use recent boostfactor models that account for tidal stripping of subhalos. This boost factor changes the limit by an factor, which is more conservative than other models sometimes used in extragalactic DM studies. While the boostfactor enhancement is fairly modest, it is still the dominant systematic uncertainty over the halo mass and concentration.
The analysis outlined in this paper can be repeated on Fermi data using published galaxy group catalogs. In particular, the Tully et al. catalogs Tully (2015); Kourkchi and Tully (2017) and the Lu et al. catalog Lu et al. (2016) provide a map of the galaxy groups in the local Universe within . Both catalogs are based primarily on 2MRS, but use different clustering algorithms and halo mass determinations. Taken together, they provide a way to estimate the systematic uncertainties associated with the galaxy to halo mapping procedure. Previous cluster studies on Fermi data Ackermann et al. (2010); Ando and Nagai (2012); Ackermann et al. (2014a); Anderson et al. (2016); Liang et al. (2016) used the extended HIghest Xray FLUx Galaxy Cluster Sample (HIFLUGCS) Reiprich and Boehringer (2002); Chen et al. (2007), which includes 106 of the brightest clusters observed in Xray with the ROSAT allsky survey. These clusters cover redshifts from ; the distribution of their factors, masses, and redshifts are shown in Fig. 7 and 8. In general, the 2MRS catalogs provide a larger number of groups that should be brighter in DM annihilation flux, so we expect a corresponding improvement in the sensitivity to annihilation signatures.
The recent advancement of galaxy catalogs based on 2MRS and other nearby group catalogs allows us for the first time to map out the most important extragalactic DM targets in the nearby Universe. This, in turn, enables us to perform a search that focuses on regions of sky where we expect the DM signals to be the brightest outside the Local Group. We present the complete results of such an analysis, as applied to data, in our companion paper Lisanti et al. (2017).
Acknowledgments
We thank S. Ando, N. Bahcall, R. Bartels, J. Beacom, P. Behroozi, F. Calore, W. Coulton, A. DrlicaWagner, D. Hooper, S. Horiuchi, A. Kravtsov, T. Linden, Y. Mao, K. Murase, L. Necib, J. Ostriker, A. Peter, T. Slatyer, B. Tully, C. Weniger, and S. Zimmer for helpful conversations. This research made use of the DarkSky Simulations, which were produced using an INCITE 2014 allocation on the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory. We thank S. Skillman, M. Warren, M. Turk, and the DarkSky collaboration for their efforts in creating these simulations and for providing access to them, and Y. Mao for providing the galaxy catalogs used herein. This research made use of the Astropy Astropy Collaboration et al. (2013), IPython Pérez and Granger (2007), Minuit James and Roos (1975) and NPTFit MishraSharma et al. (2017) software packages. ML is supported by the DOE under contract DESC0007968, the Alfred P. Sloan Foundation and the Cottrell Scholar Program through the Research Corporation for Science Advancement. NLR is supported by the DOE under contracts DESC00012567 and DESC0013999. BRS is supported by a Pappalardo Fellowship in Physics at MIT. RHW received support from the U.S. Department of Energy under contract number DEAC0276SF00515. This work was performed in part at Aspen Center for Physics, which is supported by NSF grant PHY1607611.
Appendix A  and factors for Extragalactic Sources
In this Appendix, we derive the factor relations used in the main text. We also derive the corresponding factor relations, which apply to the case of decaying DM. Although we do not make use of the decay results in the main text, we include these results for completeness because much of our main analysis can be extended to the decaying case. This Appendix is broken into three subsections. In the first of these, we detail the units and conventions used in our definition of the  and factors. After this, we derive an approximate form of the astrophysics factors for different DM density profiles and discuss the accuracy of the approximations made. We conclude with a discussion of error propagation in the factors. Note that several of the details presented in these appendices have been discussed elsewhere, see e.g., Ref. Abdo et al. (2010); Charbonnier et al. (2011, 2012); Evans et al. (2016).
a.1 Units and Conventions
a.1.1 Dark Matter Flux
We begin by carefully outlining the units associated with the  and factors. The flux, , associated with either DM annihilation or decay factorizes into two parts:
(A1)  
where is the photon energy and the ‘ann.’ (‘dec.’) superscripts denote annihilation (decay). The particle physics factors are given by:
(A2)  
where is the velocityaveraged annihilation cross section, is the DM mass, Br is the branching fraction into the channel, is the photon energy distribution associated with this channel, and is the DM lifetime. The annihilation factor assumes that the DM is its own antiparticle; if this were not the case, and assuming no asymmetry in the dark sector, then the factor would be half as large. The particle physics factors carry the following dimensions:
(A3)  
where ‘counts’ refers to the number of gammarays produced in the interaction and the is associated with the in the particle physics factors. Note that some references include this in the definition of the  or factors, but this is not the convention that we follow here.
The  and factors are defined as follows:
(A4)  
where is the subhalo boost factor. The  and factors carry the following units:
(A5)  
Combining these with Eq. A3, we find that
(A6) 
for both the annihilation and decay case. This means that is given in units of counts per experimental effective area [] per experimental run time []. In this work, we study extragalactic objects with small angular extent. So long as each object is centered on the regionofinterest (ROI), we expect that all of its flux will be contained within the ROI as well. This means that the photon counts obtained by integrating Eq. A4 over the entire sky corresponds to the total counts expected from that object in the ROI. The situation is different when treating objects with a large angular extent that exceeds the size of the ROI—e.g., when looking for emission from the halo of the Milky Way. In such cases, it is more common to divide the  and factors by the solid angle of the ROI () such that both they, and consequently , are averages rather than totals.
a.1.2 Halo Mass and Concentration
We briefly comment here on different mass and concentration definitions (virial and 200) as relevant to our analysis. Boostfactor models, concentrationmass relations, and masses are often specified in terms of 200 quantities, which must be converted to virial ones. In order to do this, we use the fact that
(A7) 
for the NFW profile Navarro et al. (1996), where is the normalization of the density profile, is the critical density, is the concentration parameter, and is the critical overdensity. For virial quantities, with in accordance with Ref. Bryan and Norman (1998), while for 200 quantities, . Therefore, Eq. A7 can be equated between the 200 and virial quantities and solved numerically to convert between definitions of the concentration.
For different mass definitions, we have
(A8) 
where the concentration definitions on the righthand side depend on and and may have to be converted between each other and we have suppressed the redshift dependence for clarity. Solving this numerically, we can convert between the two mass definitions.
a.2 Approximate  and factors
For an extragalactic DM halo, the astrophysical factors in Eq. A4 can be approximated as:
(A9)  
where the integrals are performed in a coordinate system centered on the halo, and is the comoving distance, which is a function of redshift for a given cosmology. The aim of this subsection is to derive Eq. A9 from Eq. A4 and to quantify the error associated with this approximation.
To handle the  and factors simultaneously, we consider the following integral over all space:
(A10) 
with . Here, is playing the role of a radius in a spherical coordinate system centered on the Earth. Therefore, we can rewrite the measure as
(A11) 
Next, we transform to a coordinate system (denoted by primed quantities) that is centered at the origin of the halo described by . Because this change of coordinates is only a linear translation, it does not induce a Jacobian and . Assuming that the Earth is located at a position from the halo center and the DM interaction occurs at position , then and
(A12) 
where we take and .
Eq. A12 can be simplified by taking advantage of several properties of the halo density. First, it is spherically symmetric about the origin of the primed coordinate system. Second, it only has finite support in . In particular, it does not make sense to integrate the object beyond the virial radius, . This allows us to rewrite the integral as follows:
For extragalactic objects, . As a result, we can take advantage of the following expansion:
(A14) 
where . It follows that the leadingorder approximation to Eq. A.2 is
(A15) 
We can calculate the size of the neglected terms in Eq. A14 to quantify the accuracy of this approximation. We take the parameters of the halo with the largest factor in the catalog to estimate the largest error possible amongst the DarkSky halos. For this halo, the fractional correction to the factor of the first neglected term in the expansion is for either an NFW or Burkert profile (described below), whilst for the factor it is . These values are significantly smaller than the other sources of uncertainty present in estimating these quantities and so we conclude that the approximations in Eq. A9 are sufficient for our purposes.
a.3 Analytic Relations
Starting from the approximate forms given in Eq. A9 and specifying a DM density profile , the  and factors can often be determined exactly. We will now demonstrate that the final results only depend on the distance, mass, and concentration of the halo—for a given substructure boost model and cosmology.
As a starting point, consider the NFW profile:
(A16) 
The parameter is the scale radius and dictates how sharply peaked the core of the DM distribution is. Starting from this distribution, the volume integral in the factor evaluates to
(A17)  
where is the virial concentration. To remove the normalization factor from this equation, we can write the virial mass of the halo as
(A18)  
which, when combined with Eq. A17, gives
(A19)  
Stopping here, we would conclude that the factor scales as . However, for a given and cosmology, is not an independent parameter. Using the results of Ref. Bryan and Norman (1998), we can write:
(A20) 
where is the critical density and
(A21)  
This relation can then be used to remove from the volume integral and we conclude that
(A22)  
We see the additional mass dimension required from the fact this scales as not is carried by . The dependence highlights that the annihilation flux is critically dependent upon how sharply peaked the halo is. To summarize, Eq. A22 demonstrates that the factor is fully specified by three halo parameters for a given substructure boost model and cosmology: the redshift , mass , and concentration .
The basic scalings and dependence shown above are not peculiar to the NFW profile, but are in fact more generic. To demonstrate this, we can repeat the above exercise for the cored Burkert profile Burkert (1996):
(A23) 
which is manifestly nonsingular as unlike the NFW profile. Here, and are the Burkert analogues of and in the NFW case, but they are not exactly the same. Indeed, following e.g., Ref. Bartels and Ando (2015), by calculating physically measurable properties of halos such as the radius of maximum rotational velocity for both the NFW and Burkert cases and setting them equal, we find
(A24) 
We will replace with a concentration parameter . Following the same steps as for the NFW profile, we arrive at:
(A25)  
from which we see that .
For the case of decaying DM, the approximate integral given in Eq. A9 can be evaluated independent of any choice for the halo profile. Specifically:
(A26) 
where the second equality follows from the fact that the volume integral gives the virial mass exactly. For DM decays in relatively nearby halos, the emission can be quite extended, as the flux is not as concentrated towards the center of the halo as in the annihilation case. As such, it is often useful to have a version of the extragalactic factor where one only integrates out to some angle on the sky from the center of the halo, or equivalently to a distance . In this case:
(A27)  
for the NFW profile, where we have made explicit the fact that is a function of . When , this reduces to the simple result in Eq. A26.
a.4 Propagating factor Uncertainties
We now discuss the propagation of uncertainties from inferred halo parameters into an overall uncertainty on the factor. Specifically, we will justify the form of the lognormal distribution for the factor that was used in Eq. 8 of the main text. While we focus our attention on the case of the factor for the NFW profile, the logic carries over straightforwardly to the Burkert profile, or even to the factor.
Our starting point is the approximate form of the NFW factor given in Eq. A22, which shows that the factor is determined by the redshift , mass , and concentration , for a given substructure boost model and cosmology. Therefore, the errors on these three parameters need to be propagated to determine the total error on the factor. When the redshift is determined spectroscopically, the uncertainty on is subdominant to the uncertainties on the factor that are induced by the mass and concentration. As such, we neglect the uncertainty in the redshift. If one were using photometric redshifts, however, these uncertainties would also need to be accounted for. Additionally, for nearby halos in particular, the relation between redshift and distance can have further uncertainties, since this relation is affected by local peculiar velocities.
From our studies of DarkSky, we see that the and distributions are wellapproximated as lognormal. If the factor simply scaled as the product of several lognormal distributions (), then would also be lognormally distributed. However, the dependence of the concentration parameter and boost factor on the virial mass mean that will not be exactly distributed in this way. Nevertheless, by explicitly calculating the full distribution, we confirm that it is very accurately lognormally distributed. The reason for this is that the mass dependence of the boost factor in the Bartels substructure model Bartels and Ando (2015) is very mild and additionally the dependence is subdominant in dictating the form of , beyond the dependence. Because the lognormal is considerably simpler than the full distribution, we adopt it for the factor.
The form of the likelihood that we use for the factor in Eq. 8 is the same as that used in Ref. Ackermann et al. (2015a), but stands in contrast to Ref. Ackermann et al. (2011, 2014b) with the substitution of the nominal central factor in the denominator instead of the marginalized value. The interpretation of the factor as a lognormal likelihood Ackermann et al. (2015a) rather than a lognormal prior Ackermann et al. (2011, 2014b) ensures proper normalization for all values of and, when interpreted as a maximum likelihood estimator for signal recovery, coincides with the true underlying signal strength. This is illustrated in Fig. A1, where we show the test statistic as a function of DM cross section obtained for an injected signal of at GeV after factor marginalization in the two different cases for a single halo. The left panel shows the traditional lognormal prior form of the likelihood. The maximum of this, indicated by the vertical dotted line, is offset from the true value of the underlying cross section. This discrepancy can be substantially amplified when stacking a large number of objects. Using the lognormal likelihood form (right panel), we see that the maximum test statistic associated with the recovered signal coincides with the injected signal at the true value, .
Appendix B Validity of Likelihood Approximation
In this Appendix, we address issues pointed out in the main text where our analysis procedure appears to induce small systematic discrepancies. First, we discuss the incorporation of the factor uncertainties, and then we discuss issues related to the profile likelihood procedure itself.
b.1 Jfactor Likelihood
In Fig. 9 of the main text, the bestfit cross sections are systematically higher than their injected values at high where DM detection is significant. While this could be consistent with statistical fluctuations in the factor distributions, it likely results from the fact that the assumed factor distributions, and the methods we have for calculating the central values and uncertainties, are not an exact representation of the actual body data. To demonstrate this, Fig. A2 shows the TeV injected signal plot for an observer at Location 1, where all factors are fixed to their true values. In this case, the bestfit cross section exactly matches the injected cross section at high values of , confirming that it is indeed the assumed factor distributions that induce the bias seen in Fig. 9.
b.2 Profile Likelihood Approximation
As described in Sec. III, we use an approximation to the full profile likelihood procedure to make the analysis computationally tractable. This approximation comes into play when removing the nuisance parameters associated with the astrophysical templates:
(A28) 
following the notation of Sec. III. To briefly review, the full implementation of the profile likelihood requires maximizing the for all and . Instead, we set the to their maximum values, as obtained in an initial scan where all the template normalizations are floated. In general, we find that this approximation works very well, except at high photon energies ( GeV) and at low photon energy ( MeV). We now describe the challenges incurred in more detail.
At the highest energies, the number of photons becomes statisticslimited, and it is likely that there are very few photons in a given ROI. If one of these photons happens to fall near the expected halo center, then the DM template will pick up flux in the initial scan, but all the other (astrophysical) templates will not. In this case, the bestfit values for the astrophysical templates can be near zero, even though much larger normalizations for these templates would still be consistent with the data. The problem arises when we set the to their values from the initial scan and determine the DM intensities, because there will be evidence for a signal where there should be none, since the astrophysical templates are not allowed to adjust from their nearzero values. This problem is not present in the full profile likelihood method, because in that case, one maximizes the likelihood over the when constructing the likelihood profiles as functions of the DM intensities. In particular, in the full profile likelihood method one obtains a better fit at low DM intensities, compared to that obtained in our approximation, because the astrophysical template normalizations are allowed to be higher. We stress that this is not a concern at moderate energies, where the photon counts in each ROI are large enough such that the normalizations of the astrophysical templates are always welldetermined.
The behavior described above can be observed directly in the mock data (with no injected signal). For examples, the inset plot with TeV in Fig. 9 shows the maximum test statistic as a function of the injected cross section. At low injected cross sections, the data is described by the null hypothesis and so the TS should follow a chisquare distribution. In particular, this means that the 84 percentile, which is given by the upper boundary of the red bands in Fig. 9, should asymptote to a value , while the lower part of the band should be consistent with zero. The discrepancy between the MC results and the chisquare expectation are most pronounced at high masses where the high energy bins are relatively more important.