A Search for Dark Matter in the Galactic Halo with HAWC

The High Altitude Water Cherenkov (HAWC) gamma-ray observatory is a wide field-of-view observatory sensitive to 500 GeV - 100 TeV gamma rays and cosmic rays. With its observations over 2/3 of the sky every day, the HAWC observatory is sensitive to a wide variety of astrophysical sources, including possible gamma rays from dark matter. Dark matter annihilation and decay in the Milky Way Galaxy should produce gamma-ray signals across many degrees on the sky. The HAWC instantaneous field-of-view of 2 sr enables observations of extended regions on the sky, such as those from dark matter in the Galactic halo. Here we show limits on the dark matter annihilation cross-section and decay lifetime from HAWC observations of the Galactic halo with 15 months of data. These are some of the most robust limits on TeV and PeV dark matter, largely insensitive to the dark matter morphology. These limits begin to constrain models in which PeV IceCube neutrinos are explained by dark matter which primarily decays into hadrons.


Heavy Dark Matter
For all the evidence pointing toward the existence of particle dark matter (DM), no observations have yet indicated a preferred mass scale in the dark sector. The prevailing view is that DM should be associated with new physics at the electroweak scale; however, this has so far not been realized experimentally. This motivates a re-examination of the different mass scales where DM could exist.
If the DM in the Universe is composed of a particle with a mass well above the electroweak scale (>1 TeV), then its non-gravitational nature may only be discoverable through indirect detection. At these masses direct production of DM at a collider becomes inaccessible. Furthermore the expected direct detection rate in many models becomes negligible due to the low DM number density and small DM-nucleon cross-sections. Yet if the DM can decay or annihilate to standard model (SM) final states, the stable cosmic rays resulting from these processes can leave detectable imprints in the sky. HAWC, with its ability to search for TeV gamma-rays over large portions of the sky, is ideally placed to search for such signals of heavy DM (mass >1 TeV), and in the absence of a signal set strong limits on the DM parameter space. This is the strategy we pursue to set limits on heavy DM in this work.
Theoretically, heavy DM is also well motivated. As a concrete example, if the dark sector consists only of a confining gauge theory, then there will be a dark glueball in the spectrum with a mass at the scale of confinement, which can be well above the electroweak scale [1][2][3][4][5][6][7].
At these masses it is common to focus on decaying DM, as the simplest models where the DM annihilates cannot produce the observed DM thermal relic abundance without being in violation of the unitarity bound [8,9], see also [10]. Nevertheless, models exist which evade these bounds, see e.g. [11,12], and so limits on DM annihilation at these masses should still be considered.
A further impetus in searching for DM signals at these masses comes from the detection of an astrophysical high-energy neutrino flux at IceCube [13][14][15][16]. At present no point sources have been detected in the data [17][18][19] and whilst a number of purely astrophysical explanations of the flux remain, see e.g. [20] and in particular [21] for a recent summary, the presence of a large contribution from DM cannot at present be excluded using the Ice-Cube data alone. Accordingly, a number of models have been developed which explain the high-energy neutrino flux with DM [7,. Generally these models seek to explain the highest energy neutrinos observed by IceCube-the so-called High Energy Starting Events (HESE)-however some also deal with a putative excess in the lower energy Medium Energy Starting Events (MESE) [44], which has been discussed in the context of DM [45] and astrophysics [46,47]. Nevertheless, in these models the DM decays will generically also produce photons which can be searched for in other experiments. In particular [7] considered the limits obtained with the Fermi-LAT, and demonstrated that models which seek to explain the majority of the IceCube flux appear to be in tension with the gamma-ray data. HAWC, with its ability to study gamma-rays at even higher energy than the Fermi-LAT, can also constrain this parameter space, as we demonstrate in this work.

Using HAWC to Constrain DM
The Galactic center is expected to have the largest DM signal in the sky. However, the Galactic center transits at the edge of the HAWC field-of-view, so we instead consider a region near the Galactic center but extended to higher declinations where HAWC observations are more favorable. In particular, in the work, we consider the region of the Northern Fermi Bubble. The region of the Northern Fermi Bubble is particularly useful in constraining DM annihilations and decays from the Galactic halo for HAWC because: • It is a region which extends outward from the Galactic center; • It has large solid angle for integrating the DM gamma-ray flux; • It extends to declinations of 8 degrees North (culminating at 11 degrees zenith for HAWC), so the analysis is not strongly dependent on the HAWC sensitivity at declinations with low sensitivity; • The systematic errors are well-characterized, including an understanding of the background exclusion mask and well-behaved background significance distribution [48]; and • It is known that HAWC sees no significant gamma-ray excess in this region, which minimizes confusion between DM and standard astrophysical sources.
The DM limits do not strongly depend on the chosen region-of-interest, but the characterized systematics and lack of astrophysical gamma rays in the Northern Fermi Bubble region-ofinterest makes the analysis robust. Additionally, this is an excellent example of how DM limits are based on flux limits, so a calculated flux limit can be used to limit DM annihilation and decay.
In this work, we discuss the region of interest for the analysis as well as the methods used to calculate the background. We use these values to calculate a flux limit on emission from a region toward the Galactic center. We then derive a gamma-ray flux from DM annihilation or decay which could be searched for in this region. Finally, we present the lack of observed signal as a limit on the DM annihilation cross-section and decay lifetime and discuss the implications of these limits in the context of other astrophysical experiments.

Procedure
The present work searches for a gamma-ray signal from the direction of the Galactic center (GC) using data from the HAWC observatory. A search for gamma-ray signals from the Northern Fermi Bubble region has been presented in Ref. [48]. The same analysis procedure, data, analysis bins, and search region are used here. The data for the analysis corresponds to the dates between 2014 November 27th and 2016 February 11th. Since it is required that the detector is stable for a period of 24 hours for this analysis, we use only the 290 days of the data over the 15 month period which meet this criterion. The data are binned according to f Hit,the fraction of functioning photomultipliers triggered in an air shower event. Table  1 in [48] defines the values for each analysis bin. For details about the detector see Ref. [49].
The main challenge of the analysis, as discussed in Ref. [48], is to estimate the background. There are three main factors that need to be considered for this: to distinguish the air showers produced by hadronic cosmic rays from those produced by gamma rays; to find the isotropic flux of hadronic cosmic rays through the direct integration (DI) method [50]; and to remove the effects of the large-scale cosmic-ray anisotropy [51].
For the first factor, gamma-hadron separation cuts are applied to the data. These separators select the gamma-like showers [49]. For more details on the values of these cuts, see Table 2 of Ref. [48]. For the isotropic flux estimation, a 24 hour integration period is used for the DI method. This is to ensure that the analysis is sensitive to large structures. A masking region around known sources is used in order to avoid contamination from the signal from these sources into the estimation of the isotropic background. Figure 1 shows the region used for the masking.
A subtraction procedure is implemented to mitigate the effects of large-scale anisotropy in the data. A more detailed discussion of these analyses may be found in Ref. [48].

HAWC Differential Flux Upper Limits near the Galactic Center Region
Instead of adding the total flux inside the search region as in Ref. [48], the search region is divided in six declination bands of 5 • each, as shown in Figure 2, and a total possible DM flux in each declination band is calculated. This division is done in order to take into account the expected spatial dependence of a DM signal and the declination dependence of the sensitivity of the HAWC detector.
The number of excess events in each declination band is shown in Table 1. Note that the distribution of significances is consistent with Gaussian background distribution. The fluxes in each declination band are calculated using the same procedure as in Ref. [48]. Figures 3,4, and Table 2 show the results of this analysis. The energy bins have a size of ∆ log(E/1 TeV) = 0.5. This size corresponds to the approximate width of the energy Figure 1. Regions masked when calculating the isotropic background, plotted in equatorial coordinates. The masked regions include the Galactic plane, which contains many sources of gamma-ray emission, and the bright HAWC sources Mrk 421, Mrk 501, Geminga, and the Crab Nebula. The cosmic-ray excesses in the anisotropy Regions A and B [51] are masked out to give consistent background calculation. Also, the signal region of the Northern Fermi Bubble is masked out so that any possible signal does not contaminate the background. The Galactic center is marked with a blue point response histograms of the detector (See Figure 8 of [48]). The fluxes in each energy bin are calculated as a weighted average of the overlapping analysis bins. The first energy bin is centered at 2.2 TeV, which is the median energy of analysis bin f 1 . Although the analysis of Ref. [48] included an overflow bin for all energies above 39 TeV, that is not useful in constraining spectral information and is not included here. Note that the highest-declination search bin only goes from 5 to 8 degrees declination because the Northern Fermi Bubble cuts off at higher declinations.   Table 2.
The process of DI background subtraction used in HAWC removes all isotropic backgrounds, primarily hadrons but also any isotropic gamma-ray signal. For a spatially-dependent gamma-ray signal spread over a large region of the sky, such as a signal from the Milky Way DM halo, the portion of the gamma-ray signal that is outside the search region will be subtracted as well. Because of this, we account for only the difference between the amount of signal expected within the search region to that expected outside in our analysis. However, because the DM signal considered in this analysis peaks toward the GC, this residual flux is large enough to set strong constraints.  Table 1. Number of excess events in each declination and analysis bin. See Table 1 of Ref. [48] for detailed description of the analysis bins.

Dark Matter Annihilation and Decay
To calculate the gamma-ray flux from DM annihilation or decay, information both about the particle nature of DM and its astrophysical distribution must be known. The gamma-ray flux is proportional to the square of the DM density for annihilation or the DM density itself for decay, integrated over the line-of-sight through the DM halo from the observer. This integral (with appropriate normalizations) is referred to as the annihilation J-factor or the decay D-factor. These are given by: The distance from the center of the galaxy is given by at a polar angle θ between the observer line-of-sight and the GC. The solid angle of the observation is given by ∆Ω. In this paper, we use the distance from the GC to the Earth as R = 8.2 kpc and a DM density at the Earth of ρ = 0.39 GeVcm −3 as in Ref. [52].
We consider a range of DM masses M χ above 7 TeV for annihilating DM and above 14 TeV for decaying DM. We also consider the spectrum of gamma rays produced in a single DM decay or annihilation dN/dE for several DM annihilation or decay channels. The DM velocity-weighted cross-section for annihilation σv or lifetime for DM decay τ are also needed. With these, the gamma-ray flux from DM annihilation or decay is given by

Dark Matter Halo Profiles
While the presence of DM in the galaxy is well established, the precise density distribution of the DM is not. Therefore, we consider three different DM density profiles in this paper, including profiles peaked toward the GC and those with a density core toward the GC. N-body simulations of DM without baryons are typically fit with a Navarro-Frenk-White (NFW) DM profile [53] with distribution ρ(r) = ρ s (r/r s )(1 + r/r s ) 2 . (3.6) For the Milky Way, the scale radius is on the order of r s = 20 kpc [54]. Normalizing to the local DM density, this gives ρ s = 0.318 GeVcm −3 . A different fit to N-body simulations with a shallower core is given by the Einasto profile [55,56]. For the Milky Way, this profile is given by 7) Figure 5. The three DM profiles considered for the Milky Way. The red curve is the NFW profile, the blue curve is the Einasto profile, and the magenta curve is the Burkert profile. Note that the Burkert profile has a much more cored dark matter profile and does not reach as high a density at the GC that the cusped NFW and Einasto profiles do. However, the orange region shows the range of galactocentric radii considered in this analysis. Our limits here are primarily sensitive at larger radii, at a few kpc, where the three profiles have more similar densities. Therefore, our limits are not strongly dependent on the choice of DM profile.
For the Milky Way, r s = 10 kpc and ρ s = 1.187 GeVcm −3 for the Burkert profile [54]. A comparison of these profiles, including the peaked and cored features, can be seen in Figure 5.
The J-and D-factors for each of these profiles in our signal region are tabulated in Table 3.
Although the three DM profiles considered here differ significantly near the GC, HAWC analysis is most sensitive to regions away from the GC where the three profiles are in general agreement. Therefore, this analysis is not as sensitive to profile shape as other DM annihilation analyses with the GC. DM decay analyses are also less sensitive to the shape of the DM profile, since they depend linearly on the DM density rather than quadratically. Figure 6 demonstrates how much each declination band contributes to the limit for each DM profile. For Burkert annihilation profile and all the decay profiles, the peak sensitivity to DM for HAWC is over 20 • from the GC.

The Background-subtracted J-and D-factors
The background calculation for the HAWC analysis uses a 24-hour DI technique to determine the cosmic-ray background [48]. For an extended signal, like the DM one we consider here, there is a nontrivial amount of gamma-ray emission in the background region as well. Therefore, the gamma-ray excess limits calculated in section 2 are actually a limit on the gamma-ray excess in the target region minus the gamma-ray excess in the background region DM profile Dec range  Table 3. Table of  (that is, the region outside of the mask). For each declination bin, therefore, we calculate the J-factor and D-factor over both the source and the background region. Our limit is calculated based on the difference of these two values (i.e. J src − J bkg and D src − D bkg .) The values of these J-and D-factors are given in Table 3. HAWC is currently working on background calculation methods which would minimize the amount of gamma-ray flux in the chosen background region, so the limits should correspondingly improve when those methods have matured.

Dark Matter Spectra
When discussing DM spectra, it is typical to assume a dominant annihilation or decay channel (e.g. the DM annihilates into a pair of Z bosons 100 percent of the time). The initial particles created in the channel are then allowed to decay into stable particles (including neutrinos, electrons, positrons, protons, and photons). There are three main sources of gamma rays from heavy DM annihilation or decays: 1. "Prompt" photons originating from the initial interaction occurring within the Milky Way halo. These include all primary photons produced in the initial annihilation or decay and the photons produced in the decays of the initial products. Figure 6. The relative sensitivity of this analysis versus declination. On the y-axis, the DM flux in each declination is normalized by the HAWC sensitivity for that declination (from Ref. [59]). This shows the contribution to the limit from each declination. For the Burkert annihilation profile and all three decay profiles, the peak contribution to the limit stems from declinations > 20 • north of the Galactic center (declination −29 • ) that transit closer to the HAWC zenith.
2. Secondary photons up-scattered from the cosmic microwave background or interstellar radiation fields due to inverse-Compton scattering off the prompt electrons and positrons 3. Secondary photons produced from the cascades that develop due to interactions of the primary photons with background photons (e.g. the cosmic microwave background). This effect is only important over distance scales greater than O(Mpc), much further than the Galactic center.
All three of these contribute to setting a limit in the 200 MeV to 2 TeV energy range of the Fermi-LAT, as shown in Ref. [7]. We follow the prescription in that reference for the various contributions, using the PYTHIA 8.219 software [60][61][62] with electroweak corrections enabled [63]. However, here we will only consider the prompt emission as at TeV energies it is by far the largest component at these distances and energies. Importantly, neglecting the contributions from secondary emission processes is conservative as it decreases our total expected flux and also reduces our sensitivity to the systematics associated with calculating these secondaries. In addition we neglect the attenuation of these prompt gamma rays over Galactic scales. At the energies considered and within our region of interest the attenuation is expected to be a small effect (<5%) [64].

Combining Limits in Declination and Energy
We use the J-and D-factors for each morphological assumption along with the spectra described in section 3.4 to calculate the expected gamma-ray flux from DM via equations 3.4 and 3.5. These values are calculated in each declination bin for each energy bin. With the spectrum and DM mass fixed, we have only the cross-section or lifetime as a global free scaling parameter on the dark matter flux in each bin. We perform a χ 2 test comparing these dark matter fluxes with the measured gamma-ray flux values for each declination/energy bin (from Table 2), following Appendix B of reference [65]. This gives us the 95% CL limit on either the dark matter lifetime or cross-section for that mass and channel. We then repeat this over several masses and channels to get the results shown below.

Dark Matter Limits
Following the methods of the previous section, we calculate limits on the DM annihilation cross-section and decay lifetime for HAWC observations near the GC. We show the limits for ten DM annihilation and decay channels: uū, bb, tt, W + W − , ZZ, hh, e + e − , µ + µ − , τ + τ − , and νν (an average of the three neutrino species). The annihilation limits are shown in Figures 7-8 and the decay limits are shown in Figures 9-10. The data in the region-of-interest chosen show a slight deficit compared to background, at roughly 1-1.5 sigma. Therefore, the limits shown are slightly better than would be predicted based purely on the HAWC sensitivity. Therefore, we show this predicted limit in Figure 11 as a black line. These are shown for the bb and τ + τ − DM profiles. Surrounding the black line is a green band and yellow band. These correspond to the 68-and 95-percent fluctuations away from the black line. This data has under-fluctuations for this region on the sky. In addition to this statistical uncertainty, it should be noted that HAWC has a systematic uncertainty of 50% at these high energies [49].

Discussion
The limits in Figures 7-8 are shown for three different DM profiles, the peaked NFW and Einasto profiles and the cored Burkert profile. However, because this analysis is not sensitive to the DM flux directly at the GC, these limits change by less than a factor of three between the three DM profiles. In this way, these limits are robust and much less dependent on the DM profile than Galactocentric limits done in the immediate vicinity of the GC, such as the limits from HESS [68].
As discussed in Ref. [65], the limits obtained using this analysis are approximate, due to the finite energy range of each energy bin. These effects dominate near the edge of each energy bin. Particularly for extremely hard spectra, such as monochromatic gamma-ray lines, the spectral limits can change significantly as a spectral cutoff shifts to energies outside of each energy bin. In order to correctly account for these abrupt unphysical discontinuities, we treat these as systematic errors, which we show as shaded bands around our limits.
The systematic error bands for the curves are small for softer spectra, including the hadronic and bosonic channels. Even for the fairly hard spectra of τ + τ − and νν channels, the systematic band is small. However, for the hardest channels, µ + µ − and e + e − , the spectrum cuts off sharply, making it fairly monochromatic and resulting in a wider systematic error band. This is a feature of all analyses of monochromatic (or nearly-monochromatic) spectra analyzed using energy bin-based analyses, particularly those with large energy uncertainties. As the energy uncertainty decreases and the energy bins become narrower, this systematic error is reduced.

Dark Matter Annihilation Limits
The DM annihilation limits are shown in Figures 7-8. Although the GC peaks at a zenith angle of 48 degrees, a large zenith angle for HAWC, the large expected gamma-ray flux from DM from the GC still provides strong limits on the DM annihilation cross-section. While this analysis is only calculated for energies below 39 TeV, the limits are within a factor of  2 of the anticipated HAWC sensitivity plots of Ref. [69] at 1-10 TeV DM masses. This is consistent with those plots being 5-year sensitivity curves and these limits made with less than two years of data.
For energies higher than tens of TeV, the current HAWC analysis has large systematic errors, so we limit the results here to energies below 39 TeV. More precise energy estimators are being vetted which will significantly improve the understanding of HAWC energies above this value. The precise HAWC sensitivity above these energies is currently under study, but the limits at higher masses should improve by at least a factor of 3 with HAWC flux limits up to 80 TeV. A more optimized analysis of the region-of-interest surrounding the DM signal and more optimized choice of background region should improve these limits further as well.
Below hundreds of TeV in DM mass, these annihilation limits are not as strong as those from the HESS collaboration study of the GC [68]. However, our limits are much less sensitive to the DM profile, as discussed above. Therefore, these limits complement current limits for cases where the DM halo is cored, such as a Burkert scenario. Additionally, these Figure 9. 95% confidence level lower limits on the DM decay lifetime for the hadronic and bosonic DM annihilation channels considered in this analysis. Limits for a Burkert (bottom, green), Einasto (top, orange), and NFW (middle, blue) profile are shown. Segue 1 dwarf galaxy limits from MAGIC [66] (lower left) are dashed in red. HAWC limits from dwarf galaxies (upper left) are shown dashed in purple. IceCube limits [67] are shown in brown (on right). The red star is a model in which the IceCube neutrinos are produced in decaying DM [7]. conservative limits are within factor of 3 of HESS limits at 100 TeV for leptonic DM channels. For hadronic and bosonic DM channels, our limits are within a factor of 30 of HESS limits at 100 TeV and should be within a factor of three of the HESS limits above 1 PeV, extrapolating the published HESS limits. Therefore, at these high masses, the DM is being well-constrained regardless of DM profile. However, it should be noted that our limits are still only probing cross-sections 2-3 orders of magnitude larger than for thermal DM.

Dark Matter Decay Limits
The DM decay limits are shown in Figures 9-10. For DM decay, which has a more spatiallyextended profile than annihilation, these HAWC limits are some of the most constraining. More pointed instruments typically cannot observe the wide fields-of-view needed to do these limits with the GC and have instead focused on Galaxy clusters or dwarf galaxies for DM decay limits. Compared to one of the strongest of these, the MAGIC observations of the Segue 1 dwarf galaxy [66], the HAWC limits are stronger over the entire mass range studied, above 14 TeV DM mass. Figure 11. Demonstration of the statistical errors compared to the limits for bb and τ + τ − DM annihilation and decay channels considered in this analysis. Limits for an Einasto (orange) profile are shown as in Figures 7-10. Segue 1 dwarf galaxy limits from MAGIC [66], GC limits from HESS [68], HAWC limits from dwarf galaxies, and the model for IceCube neutrinos from decaying DM [7] are also shown as in Limits based on IceCube observations of the Galactic halo also can constrain PeV-level DM masses. In particular, Ref. [67] showed that neutrino limits are competitive with gammaray limits on DM for DM masses above 100 TeV. In Figures 9-10, these are shown as the brown lines. Below 100 TeV, the IceCube sensitivity falls quickly, so our limits are much better. However, our limits are similar to the IceCube limits from 1-10 PeV in DM mass for the hadronic and bosonic DM channels. For leptonic and neutrino dark matter channels, many neutrinos are produced, so the IceCube limits are stronger than those in our analysis.
A study of Fermi-LAT data that includes Inverse Compton emission from charged particles emitted by the DM decays has also shown strong DM decay limits in this mass range [7]. Most of the limits in this paper are 1-2 orders of magnitude worse than the Fermi-LAT limit, but our analysis only includes the prompt gamma-ray emission from the decays. Our analy-sis is insensitive to the details of particle propagation and the low-energy tails of the decay spectra. Our limits can be viewed as more robust, more conservative limits to complement theirs. Also, at the highest masses, the addition of higher energy bins, more high-energy data, and better analysis in the future should improve the HAWC limits significantly and make these conservative limits even more constraining.
Also included in these plots is an interpretation of the IceCube neutrino excess coming from DM [7]. For the hadronic channels, all three DM profiles are in tension with this interpretation of the IceCube excess. For bosonic profiles, the more peaked dark matter profiles do constrain a DM interpretation of the IceCube excess, while the less-peaked Burkert profile does not quite. In all three DM profiles, the decay to Z bosons is currently allowed, though the HAWC limits are within tens of percent of testing that hypothesis. For the leptonic channels, which readily produce neutrinos, the HAWC limits do not yet constrain the IceCube excess. However, the limits on decay into τ leptons and e + e − will test those hypotheses soon.

Conclusions
In this analysis, we searched for DM annihilation and decay signals coming from near the GC with HAWC. Due to the lack of significant excess in this region, it has been possible to constrain the DM annihilation cross-section and decay lifetime for TeV-PeV DM masses. These limits were robust across different DM profiles, allowing us to make morphologyindependent constraints on the DM. Particularly for DM decay, we provide some of the strongest and most robust limits on the DM decay lifetime. With improved understanding of the HAWC high-energy sensitivity and optimized DM analysis, HAWC DM searches from the GC should be even stronger in the future.