Searching for Dark Matter Sub-structure with HAWC

Simulations of dark matter show a discrepancy between the expected number of Galactic dark matter sub-halos and how many have been optically observed. Some of these unseen satellites may exist as dark dwarf galaxies: sub-halos like dwarf galaxies with no luminous counterpart. Assuming WIMP dark matter, it may be possible to detect these unseen sub-halos from gamma-ray signals originating from dark matter annihilation. The High Altitude Water Cherenkov Observatory (HAWC) is a very high energy (500 GeV to 100 TeV) gamma ray detector with a wide field-of-view and near continuous duty cycle, making HAWC ideal for unbiased sky surveys. We perform such a search for gamma ray signals from dark dwarfs in the Milky Way halo. We perform a targeted search of HAWC gamma-ray sources which have no known association with lower-energy counterparts, based on an unbiased search of the entire sky. With no sources found to strongly prefer dark matter models, we calculate the ability of HAWC to observe dark dwarfs. We also compute the HAWC sensitivity to potential future detections for a given model of dark matter substructure.


JCAP07(2019)022
Abstract. Numerical simulations show that the dark matter halos surrounding galaxies are expected to contain many over-densities or sub-halos. The most massive of these sub-halos can be optically observed in the form of dwarf galaxies. However, most lower mass sub-halos are predicted to exist as dark dwarf galaxies: sub-halos like dwarf galaxies with no luminous counterpart. It may be possible to detect these unseen sub-halos from gamma-ray signals originating from dark matter annihilation. The High Altitude Water Cherenkov Observatory (HAWC) is a very high energy (500 GeV to >100 TeV) gamma ray detector with a wide fieldof-view and near continuous duty cycle, making HAWC ideal for unbiased sky surveys. We perform a search for gamma ray signals from dark dwarfs in the Milky Way halo with HAWC. We perform a targeted search of HAWC gamma-ray sources which have no known association with lower-energy counterparts, based on an unbiased survey of the entire sky. With no sources found to strongly prefer dark matter models, we calculate the ability of HAWC to observe dark dwarfs. We also compute the HAWC sensitivity to potential future detections for a given model of dark matter substructure. Assuming thermal dark matter, we find the corresponding J-factor of a dark dwarf required to reach the HAWC detection criterion is 5.79 × 10 20 GeV 2 cm −5 sr for one particular set of dark matter assumptions. HAWC is found to be able to competitively constrain dark matter annihilation from discovered halos with J-factors on the scale of 10 19 GeV 2 cm −5 sr or greater, with better constraints obtained on dark matter models with > 10 TeV masses and sources that transit overhead.

Introduction
The mass of the known universe is dominated by dark matter. Its effects can be seen in phenomena such as galactic velocity dispersion and gravitational lensing through galaxy clusters [1]. These effects cannot be explained by normal, luminous matter within the known laws of gravity. This leads to the hypothesis that halos of dark matter surround the luminous components, providing the mass necessary to explain these effects [2]. However, the composition of the dark matter remains unknown.
One popular class of dark matter candidates are weakly interacting massive particles or WIMPs. This hypothetical particle interacts with Standard Model particles via a weakscale force [2]. Under this hypothesis, dark matter annihilation can produce Standard Model particles through weak interactions. These particles, in turn, produce gamma rays mainly via pion decay , but also through inverse Compton scattering of photons off produced electrons and positrons pairs. Searching for dark matter through these photonic signals is known as an indirect dark matter search.
Searches for dark matter signals typically focus on regions known to be heavily dark matter dominated, such as the dwarf spheroidal galaxies surrounding the Milky Way [3]. These dwarf galaxies form as a result of dark matter substructure (sub-halos) within the Milky Way halo. These sub-halos form a gravitational nucleus around which stars can form, leading the formation of a dwarf galaxy. However, models of the evolution of the early universe predict there to be far more of these sub-halos than currently observed. Simulations of dark matter show that many more sub-halos are expected to exist than are optically observed [4]. However, lower-mass sub-halos are not expected to have a luminous counterpart due to lack of gravitation power and processes that suppress star formation [5,6]. We refer to this -1 -JCAP07(2019)022 unseen substructure as consisting of "dark dwarfs" -nearby dark matter halos with no luminous counterpart that would escape current optical observation. It should be noted that the inclusion of normal matter in these simulations reduces the amount of expected substructure, but non-luminous sub-halos are still expected [7].
Dark matter annihilation in these unseen dwarfs could produce gamma rays. Therefore, it may be possible to detect dark matter signals from these previously unobserved sub-halos by performing an all-sky survey for gamma-ray signals.
In this paper, we first investigate resolved HAWC gamma-ray sources that might be associated with dark matter sub-halos. Finding no sources that significantly preferred the dark matter hypothesis, we calculate characteristic upper limits on dark matter annihilation for an unbiased sample of the HAWC sky. Using this information in combination with simulations of dark matter substructure formation, we then compute the HAWC sensitivity to detections of dark dwarf signals.

Dark matter and gamma rays
The expected differential photon flux (per unit energy) from a dark matter halo is described by the following equation: where σv is the velocity-weighted annihilation cross section and M is the dark matter mass. The J-factor is defined as: an integral of the squared dark matter mass density profile over the line of sight and solid angle of the observation. The quantity dN dE is the gamma-ray spectrum from a single dark matter annihilation [8]. We use PYTHIA 8.2 to calculate this function by simulating dark matter annihilation and recording the number of gamma rays produced; our models assume 100 % branching ratios of dark matter into individual standard model particle channels [9]. We include the bb channel since this has been extensively studied in other experiments as well the τ + τ − channel because this is the heaviest solely leptonic channel available. An example of dark matter spectral shapes, showing the characteristic hard energy cut-off at the dark matter mass, is shown in figure 1.

HAWC
The High Altitude Water Cherenkov (HAWC) detector is a gamma-ray observatory located at Sierra Negra, Mexico. Consisting of an array of 300 water Cherenkov detectors and covering an area of 22,000 m 2 , it is able to detect gamma rays by the air showers they produce in the atmosphere. Air showers are created when charged particles or gamma rays interact with Earth's atmosphere, producing cascades of lower-energy particles. HAWC detects these secondary particles by the Cherenkov light they produce while passing through the water Cherenkov detectors, which is then detected with four photomultiplier tubes (PMTs) mounted at the bottom. Timing information is then used to reconstruct angle of arrival, and the distribution of charges is used [10].  . Example dark matter energy spectra for two annihilation channels (bb and τ + τ − ). The spectra assume a mass of 100 TeV and cut off sharply at this energy. These are used in eq. (2.1) and they show the same characteristic shape.
The majority of air showers detected by HAWC come from charged cosmic rays (mostly protons). Therefore, cuts are applied to separate the hadronic events from the gamma rays before the data is analyzed. With cuts, 99% of hadronic events are rejected at the highest energies [10]. To estimate the remaining background after cuts, we use a technique called direct integration. Events are integrated within two hours in right ascension of a point to estimate the associated background [10]. Our signal at each point in the sky is then the number of excess gamma rays above this estimated background level. Due to atmospheric effects HAWC is sensitive to the declination at which a source transits, with the best sensitivity obtained from points that transit overhead [10].
HAWC is sensitive to gamma rays between 500 GeV and >100 TeV and is well-suited for detecting signals from multi-TeV dark matter masses. In addition, HAWC operates on a near-continuous duty cycle with a wide field-of-view that makes it ideal for performing survey-style observations [10]. HAWC does not need to be oriented to observe a source and collects data from 2/3 of the sky each day. With these properties, we perform an unbiased search for dark dwarfs.
The algorithms used in this analysis have relatively broad energy resolution (see figure 2). Therefore, rather than directly reconstruct particle energy, HAWC uses a forward folding method to fit spectra [10]. Events are binned not in energy but in f hit , the fraction of available PMTs that detect signal in a shower event (see figure 2 for bin definitions and their corresponding energy distributions). The observable spectrum in f hit space is sensitive to the energy distribution. To fit energy spectra, we convolve the HAWC f hit response to gamma-ray spectra with different spectral parameters to convert the spectra to f hit space. Using these f hit spectra, we find the energy spectrum which best matches the measured f hit spectrum with maximum likelihood methods [10]. Some examples of fits to such a forward folded spectral fitting are shown in figure 3 for a resolved HAWC source; see section 4.1 for details on the spectra used. Note that future HAWC analyses will use recently developed algorithms to reconstruct energy directly.  HAWC fits the energy spectra of its data using a maximum likelihood ratio test with respect to background. For an assumed model, we fit to maximize the test statistic (TS) defined in eq. (3.1).
where L 0 is the likelihood for the null hypothesis and L is the likelihood for signal. See appendix A of ref. [3] for details on how the likelihoods are evaluated for HAWC data.

Search method
We begin our analysis by searching 760 days of HAWC data for gamma-ray excesses that could potentially originate from dark dwarfs. We consider all sources with a TS (eq. (3.1)) of 25 or greater in the 2HWC catalog, had a galactic latitude, |b| > 5 • and have no known astrophysical counterpart reported [11]. Because these sources had no known normal-matter association and are so far off the Galactic plane, each one is a potential dark dwarf. Note that although all of these sources had a TS greater than 25 in the 2HWC catalog, three of these sources have TS values less than 25 in our new dataset of 760 days. No additional sources with TS greater than 25 were found in the 760 day dataset. Therefore, we only consider these four sources in this analysis. The sources are listed in table 1, where RA is the right ascension, Dec is the declination, and the radius is the angular size of the disk source hypothesis used in the 2HWC catalog. A radius of zero indicates a point source (smaller than the HAWC angular resolution). Sources with larger angular extent would indicate dark dwarfs that are either close to Earth or have large spatial extent. Because we use a different data set than the 2HWC catalog (which consisted of 507 days of data rather than 760), the maximum TS values for the power law hypothesis will differ from those reported in the catalog consistent with fluctuations from additional data [10]. If the charged cosmic ray background fluctuates to a higher value, but -4 -JCAP07(2019)022 Figure 3. Sample best-fit spectra plotted both in energy (top) and f hit bins (bottom) to one of the unassociated sources (2HWC J1309-054) considered in this analysis for its best-fit simple power law, cut-off power law, and dark matter spectrum. While the spectra appear distinct in energy-space, they are of comparable goodness of fit to the data in f hit bins. The top right plot also shows the statistical uncertainties of the best fits as correspondingly-colored shaded regions, indicating that the range of allowed fits overlaps substantially. We show the f hit bin fits both as terms raw counts (left) and scaled by the estimated error in the data (right). The f hit bin plots (bottom) are calculated summing events over a fixed angular disk in the sky rather than the full spatial-likelihood calculation. These are for visualization purposes only. The shown best-fit spectra and TS (eq. (3.1)) values are those from the full spatial-likelihood calculation.
not the signal, this can lower the TS relative to that reported in the catalog, which is the case for three of our four sources (see table 2 in the following section.) To test these sources for possible dark matter signals, we compare fits of the spectra predicted by eq. (2.1), to those predicted by astrophysical spectra (eq. (4.1) and eq. (4.2)).
Eq. (4.1) is a simple power law where flux is proportional to energy raised to a power γ, while eq. (4.2) is a cut-off power law with the addition of an exponentially decaying factor, with the cut-off energy X. In both cases E 0 is the pivot energy, which we fix to 7 TeV . This choice of parameterization has no effect on the actual spectral shape but does effect the correlation between the analysis bins, and we choose this value in order to minimize the correlation (see  Table 1. Location and spatial data for potential dark dwarfs considered in this analysis. Radius refers to the angular size of the disk source hypothesis (zero corresponding to a point source). All four of these sources were seen in the 2HWC catalog with maximum of TS > 25 [11] and no additional resolved sources were found in the dataset used in this analysis. Only one of these sources (2HWCJ 1040+308) still remains at TS > 25 in the 760 day dataset.
ref. [10]). In the case of the dark matter hypothesis, we do not assume density profiles and therefore do not independently compute J (eq. (2.2)). Instead, we treat σv J as a single free parameter, with dark matter mass being the other. For each of our models, the null hypothesis is that the spectrum is background-only. Each alternative hypothesis is tested separately and the TS measured. We check to see if the dark matter hypothesis produces a TS at least has high as the most favored (highest TS) power law. If this is the case, further analysis with a closer examination of the likelihood distributions will be used to determine if the dark matter hypothesis is statistically favored, and the observed emission a possible detection of a dark dwarf.

Fit results
For each source, we fit a power law, cut-off power law and dark matter spectrum. The fits are performed using the Multi-Mission Maximum Likelihood framework (3ML) [12,13]. The dark matter fits have two free parameters: the J σv pre-factor and M , the dark matter mass. We vary spectra generated assuming different M and annihilation channels calculate the J σv that maximize the likelihoods in each case. The best fits for the spectra are shown in figure 4 and summarized in table 2.; only the dark matter spectrum with the highest TS is plotted in each case.
We observe no significant improvement in the TS for the dark matter spectra compared to the power laws. In each case, the fits are either comparable or slightly worse (lower TS) for the dark matter hypothesis. It should also be noted that the TS values are nearly indistinguishable between the two power law hypotheses, indicating that the sources lack sufficient statistical power to differentiate even nested spectral shapes. Therefore, we cannot declare the dark matter hypothesis to be favored. However, it is important to note that these results do not exclude the possibility of dark matter signals from the surveyed sources; they simply do not favor the dark matter hypothesis over other astrophysical spectra.

Characteristic upper limits across the sky
With no clear detections of dark matter sub-halos in the HAWC data, we proceed to estimate characteristic upper limits on dark matter annihilation within the HAWC field-of-view. Characteristic upper limits refer to the average dark matter σv which HAWC could exclude at 95% confidence as a function of declination, assuming a dwarf with known J-factor were discovered at that declination and no gamma-ray excess were observed. We estimate our    Table 2. Summary of test statistics (TS) and ∆ TS for the three flux models considered in this paper: power law (PL), cut-off power law (CU), and the best-fit dark matter spectrum (DM). See figure 4 for details on best-fit spectra for each source.
characteristic limits by finding the average flux which could be excluded by HAWC at 95% confidence at each declination. If a dark dwarf were discovered by another instrument, these characteristic limits could be used to give the expected value of the corresponding dark matter annihilation upper limit. Note that for any one given discovery of a dwarf, fluctuations in the background at that location would vary the corresponding limit from the characteristic limit.

JCAP07(2019)022
To calculate the characteristic upper limits we use grid for 0 • < RA < 360 • in steps of one degree and −10 • < Dec < 55 • in steps of five degrees. We chose the declination spacing of five degrees to account for the HAWC detector response to different declinations, which does not substantially change on the scale of five degrees. We chose the right-ascension due to the HAWC point spread function (PSF), which has a width of approximately one degree in the lowest f hit bin [10]. A spacing of one degree ensures each sampled point is roughly independent (little overlap in the PSF) while still obtaining a representative sample of points at each declination. We perform these calculations assuming point source hypotheses since the expected characteristic size of a dwarf is on the scale of the HAWC PSF width, although some dark dwarfs may have slightly larger extent [14,15]. Note that the points analyzed in section 4.1 were drawn from a search of all points in the entire HAWC sky [11], and all possible sources of dark matter have already been identified so there is no chance an additional TS = 25 excess will be missed in this analysis.
To avoid contamination from luminous matter gamma-ray emission, we exclude any points within five degrees of known luminous matter gamma-ray sources reported in the 2HWC catalog [11], as well as the Galactic plane. As the HAWC PSF has a width of approximately one degree at the lowest energies (and narrows at higher energies), we choose five degrees in order to exclude both the main emission and the tails due to the PSF as well as emission from highly extended sources [10]. An example significance distribution is shown in figure 5 and is consistent with background, indicating we have removed significant source contamination and are using a sufficiently unbiased sample of the sky. Then, we fit each of our sampled grid points with spectra from a range of assumed dark matter masses and annihilation channels, assuming point sources. As we expect the dark dwarfs to be relatively low-mass and/or far away, we fit each pixel with a point source model. Since we do not assume any given dark matter density profiles or distances to the satellites, our characteristic limits are on J σv rather than just σv .

Characteristic limits
Since the HAWC sensitivity is highly declination-dependent, we will report our characteristic upper limits as both a function of dark matter mass and declination, rather than showing the individual limits for each sampled point. We estimate our characteristic 95% confidence level upper limits using the distribution of fitted fluxes. Each best-fit J σv value in a given declination is placed in a histogram, and we select our characteristic upper limit on J σv as that which is greater than 95% of the best-fit values. In practice, this is equivalent to averaging the individual limits obtained above. Our characteristic limits are plotted in figure 6.
With these characteristic limits, should a new dwarf galaxy be discovered and its J-factor measured, one could find the expected σv upper limit by scaling the limits in figure 6 by the corresponding J-factor. An example is shown in figure 7 assuming a fairly massive dwarf galaxy (with J = 10 19 GeV 2 cm −5 ) were discovered at each sampled declination. Assuming an even more massive or close dwarf with J = 8.7 × 10 19 GeV 2 cm −5 were discovered at our best declination (20 • ) our expected upper limits would become competitive with the HAWC combined dwarf-spheroidal upper limits. In this case, for assuming bb, M = 10 TeV spectrum, our upper limit on σv is 3 × 10 −23 cm 3 s −1 comparable to the corresponding value in the dwarf-spheroidal analysis [3]. In order to become competitive with the corresponding Magic upper limits from observations of Segue 1 (4.33 × 10 −24 cm 3 s −1 ) [16], the J-factor would need to be at least 6 × 10 19 GeV 2 cm −5 .  . Significance (as defined in the 2HWC catalog [11]) distribution of the grid points from all strips of declination used in the characteristic limit calculations (Observed) superimposed over the distribution expected from background-only (Expected). The likelihoods used for this plot were calculated assuming a 10 TeV, bb channel dark matter spectrum. The shape is roughly consistent with background, indicating no excess remains after resolved sources from the catalog are removed, and that the sample is an unbiased representation of the remaining sky. Note that the actual upper limit derived from a discovered dwarf will depend on fluctuations in the HAWC data at the dwarf location (see the analysis in ref. [3]). The limits shown here characterize the sensitivity of HAWC to dark matter emission at different declinations and show how constraining we can expect these future analyses to be.

Statistical variations and expected limits
To verify our characteristic limits, we also calculate the expected statistical variation in the upper limits. To do so, we re-calculate the upper limits at each declination using a series of pseudo-maps generated using our measured background. Each pseudo-map is created by injecting a simulated signal drawn from a Poisson distribution about the measured background in each bin. We then histogram the pseudo-map limits and find the boundaries that contain 68% and 95% of the distribution, which demonstrates the possible statistical variation of true limits compared to the characteristic limits. We use simulated pseudo-maps rather than our actual limit distribution because there are not enough independent points in each declination bin to rigorously find the statistical variation. The characteristic limits we calculated in section 5.1 are consistent with the expected (median) pseudo-map limits, as would be expected from a background-dominated distribution.
The resulting statistical variations and median (expected) limits are shown in figure 8, superimposed over the corresponding calculated limits from figure 6. The size of our statistical error bars was found to be roughly independent of the spectral assumptions and declinations. Therefore, we show only a sample of error bars for select declinations and assuming a 2 TeV mass, bb dark matter spectrum.   . Characteristic upper limits on σv for the bb and τ + τ − channels assuming a dwarf galaxy with J = 10 19 GeV 2 cm −5 sr were discovered at each declination. These are obtained from scaling the limits in figure 6 by this J-factor. For reference, the canonical thermal σv value is 2.2 × 10 −26 cm −3 s −1 [17]. Individual upper limits will depend on statistical fluctuations at the location of a discovered dwarf and its measured J-factor.  Our characteristic limits as a function of declination reflect the HAWC sensitivity curve [10], with the most constraining limits found at points that transit directly overhead. This analysis gives an estimate of the HAWC sensitivity specifically to dark matter spectra. In the following section, we use this estimate in conjunction with simulations of substructure to compute the HAWC sensitivity to detection of flux from an ensemble of dark dwarfs. It should be emphasized that the results in this section show characteristic limits consistent with expected limits within a declination band. The individual upper limits on flux from newly discovered dwarf galaxies will have statistical fluctuations as shown in figure 8.

HAWC sensitivity to modeled dark dwarfs
The HAWC sensitivity estimates of the section 5 can be applied to ensembles of dark dwarfs, not just individual dark dwarfs. Here we calculate the HAWC ability to detect dark dwarfs based on models of dark matter sub-structure across the sky. This will give an estimate of how likely HAWC is to observe a dark dwarf in its field-of-view.

Modeling dark matter
We use the clumpy (version 2015.06) package to obtain dark matter substructure models. clumpy models both the main (smooth) dark matter halo as well as dark matter subhalos (clumps) and computes the corresponding J-factor at each modeled point [18]. As the substructure distribution is not perfectly constrained by current observations, we use a set of characteristic halo models sample of 100 Monte Carlo trials in each case. For each trial,  Table 3. Parameters used in the assumed dark matter density profiles (eq. (6.1) and eq. (6.2)). R and ρ are respectively the distance from the sun to the Galactic center and the local dark matter density of the solar system and determine the scale density ρ s . The scale radius r s is the radius such that ρ(r s ) = ρ s and α determines the slope of the Einasto profile (the Burkert profile does not set the slope as a free parameter).
clumpy creates a healpix [19] map of the J-factor associated with each point given a user defined integration angle, which we use to determine the expected halo locations. The exact behavior of dark matter density profiles towards the center of a halo is not well-constrained, so we use two different parameterizations consistent with numerical simulations and observational data. Many fits to numerically simulated halos favor the Einasto density profile shown eq. (6.1), which is characterized by a cuspy shape towards the halo center [20,21].
Here, ρ s is a normalization constant on the dark matter mass density determined by the total halo mass. r s is the characteristic scale radius of the halo and α determines the profile's curvature. Note that in the case of the Galactic halo clumpy determines this constant internally using the distance from the sun to the galactic center (R ) and the local dark matter density (ρ ) [18]. All parameter values are reported in table 3. Emperical measurements of dark matter density profiles through observation of the baryonic component of the halos may favor a Burkert profile [22], shown in eq. (6.2) and characterized by a cored center. ρ(r) = ρ s (1 + r/r s )(1 + (r/r s ) 2 ) (6.2) Here, ρ s and r s are again the density normalization and scale radius and are given the same values as in the Einasto profile for the main halo component (see table 3). For the simulated sub-halos, ρ s and r s are determined as a function of the total subhalo mass M sub using the concentration parameter formalism defined in ref. [23]. In this framework, the halo shape parameters are related to the total halo mass by the concentration parameter c vir defined as: where r −2 is the radius at which the logarithmic slope of the halo is equal to -2 and R vir (M sub ) is the virial radius of the halo. The virial radius is given as the following function of M sub : where the values of the constants are given in ref. [23]. Various models exist for the functional form of c vir (M sub ) and we use the one provided by ref.  Table 4. Parameters used in clumpy simulation of the Galactic substructure mass function and spatial distribution. M min and M max are respectively the maximum and minimum subhalo masses considered (in units of the solar mass M ) and n is the slope of the power law assumed for the mass function. The spatial distribution is modeled by eq. (6.1) with the r s and α values reported here and ρ s choses to normalize to unity. The values reported here are taken from ref. [26] and ref. [28].
where M is the solar mass, z is the halo redshift, and the coefficients C i are taken from eq. 3 of ref. [24]. It should be noted that to eq. (6.5) was fitted in a simulation assuming the dark matter density follows the NFW profile [24], and the coefficients could differ had the simulation assumed an Einasto or Burkert profile. Recent work indicates subhalos may be more concentration than the main component [25]. The model we report here is therefore relatively conservative. For a given halo of mass M sub , R vir and c vir are calculated and used to determine r −2 using eq. (6.3). The corresponding r s is then obtained from r −2 depending on the functional form of the density profile. For the Einasto profile, r s = r −2 ; for the Burkert profile, r s ≈ r −2 /1.52. Once r s is determined for a given halo mass, ρ s is determined by normalizing the integrated density profile to the total mass. Put in equation form: where f (r) is the functional form of the density profile (eq. (6.1) or eq. (6.2)). Following the Aquarius simulation, we assume the number distribution of sub-halo masses follows a power law form ( dN dM sub ∼ M −n sub ), where n is the power law index, and set the minimum and maximum subhalo masses considered accordingly [26]. The parameters used here are summarized in table 4 and the resulting total fraction of dark matter mass contained in substructure is then 18%. It should be noted that the values determined here from Aquarius were obtained using a set of cosmological parameters that differ from those most recently obtained by the Plank experiment [27].
The spatial probability distribution of the dwarfs is modeled by the Einasto profile, where ρ s is chosen such that the profile is normalized to unity [28]. See table 4 for the values used. Note that in all cases, we only consider dwarfs with simulated J-factors > 10 16 GeV 2 cm −5 sr. In both cases, the sub-halo density profiles have the same functional form as the main halo profile [18].
A sample simulation is shown in figure 9. This map was generated assuming an Einasto profile for the main halo and sub-halos and is one of the 100 trials used. Under this profile assumption, we expect to observe roughly 50% more high J-factor (> 10 19 GeV 2 cm −5 sr) dwarfs than have been observed.

Detection thresholds
For each clumpy trial, we compute the HAWC sensitivity to gamma-ray signals from the generated J-factor map of the sky. As the HAWC sensitivity is dependent on declination, we bin the sky sensitivity into declination bands of five degrees. For each dark dwarf we generate the corresponding TS (eq. (3.1)) profile for a given mass and channel as a function  of assumed σv , using the observed background from a randomly chosen pixel within the corresponding declination band.
We then calculate σv such that at least one of these simulated sub-halos' TS value changes by 25 from the background-only case. This definition matches that used in the 2HWC catalog (i.e. a source is detected if TS ≥ 25 [11].) See appendix A of ref. [3] for additional details on evaluating TS values with HAWC. We repeat the process described above for each of the one-hundred clumpy trials and obtain an ensemble of σv sensitivity values. In figure 10, we report the median of these values as a function of dark matter mass and annihilation channel. The data in figure 10 should then be interpreted as the σv value such that HAWC would have a 50% chance of detecting dark dwarf emission at TS = 25 within the ensemble of possible realizations of dark matter substructure.
It should be noted that these results are not intended to be as constraining as the 95% upper limits shown for a hypothetical new dwarf in figure 7. Resolving a source at TS = 25 requires a much larger σv than can be excluded at TS = 2.7 (the corresponding number for 95% one-sided limits).
In the background-only case, ignoring effects from fluctuations, TS is approximately proportional to σv 2 , which can be seen by performing a Taylor expansion on the TS as a function of scale factor (eq. A4 of ref. [3]). Therefore, one would expect the sensitivity for this analysis to be a factor of approximately √ 25/ √ 2.7 ≈ 3.08 less constraining than the 95% CL threshold if considered with the same declination and J-factor.
To verify this expectation we examine the results for one of our one-hundred clumpy trials, assuming a dark matter mass of 10 TeV. In this trial the σv threshold was 6.06 × 10 −22 cm −3 s −1 , and resulted in a TS = 25 detection of a simulated dwarf with J = 10 19.23 GeV 2 cm −5 sr at a declination of five degrees. The corresponding characteristic upper limit from figure 6, scaled by this J-factor gives 95% CL limit on σv of 1.94 × The TS = 25 case is a factor of 3.1 times less constraining than the 95% CL case, matching the estimated change.

J-factor sensitivity
Another possible interpretation of the results in section 6.2 is to determine the minimum J-factor required for a dark dwarf to be detected (assuming some σv value and dark matter model). Rather than associating each likelihood profile with a simulated J-factor, one could instead fix the value of σv and compute J such that the TS reaches 25. We show a sample calculation here for the τ + τ − channel assuming a dark matter mass of 10 TeV (the best-constrained channel according to figure 6), and a velocity-weighted cross-section of 10 −24 cm 3 s −1 (roughly the highest value consistent with the upper limits set in ref. [3]).
Using TS profiles generated at a declination of 20 (our most sensitive declination bin), we find the average J-factor required for a TS of 25 to be 5.79 × 10 20 GeV 2 cm −5 sr. For comparison the Segue 1 dwarf spheroidal galaxy was found to have a J-factor of 1.1 × 10 19 GeV 2 cm −5 sr [16]. Given that Segue 1 is at a distance of 23 kpc, and that the Jfactor is approximately proportional to the square of distance, a dwarf of comparable mass and scale would need to at most ∼ 3 kpc away to guarantee a detection by HAWC.

Conclusion
We were able to use the excellent survey capabilities of HAWC to perform an unbiased search for dark matter annihilation originating from substructure. We observed no significant excesses that could not be equally well explained emission from luminous matter. Though this does not rule out the possibility of the signals originating from dark matter, the HAWC data does not significantly favor the dark matter hypothesis over other astrophysical hypotheses.
We then calculated characteristic limits on flux from dark matter as a function of declination -the sensitivity of HAWC to set limits from discovered dwarf galaxies. Since we did not assume dark matter density profiles for the hypothetical dwarfs, these characteristic limits were done for J σv rather than σv . These characteristic limits match what is expected from the HAWC sensitivity, with the most constraining limits coming from the points which transit directly overhead of HAWC. As new dwarf galaxies are discovered in our survey region, we will be able to estimate expected constraints on σv by scaling our current limits by the measured J-factors. (Dark matter limits from any individual dwarf are dependent on statistical background fluctuations at the dwarf location.) With these calculations, we computed the HAWC sensitivity to gamma ray emission from dark dwarfs. Using simulated realizations of dark matter sub-halos, we found the σv value which would give an expectation of at least one dark dwarf at TS = 25 50% of the time. Although the ensuing detection thresholds appear less sensitive than the 95% CL limits shown in previous HAWC analyses, they are consistent with what is expected for the more stringent TS = 25 criterion used. A future analysis could improve this sensitivity by taking a combined likelihood approach, finding a total change in TS of 25 above background, summed over all sub-halo contributions. The combined additional significance above background could reach TS = 25 at a much lower σv . In addition, HAWC is continuously taking additional data and improving its reconstruction algorithms. With these additional tools, it will be possible to gain even more sensitivity to dark matter substructure in the gamma-ray regime.  . Estimated HAWC detection threshold σv values for a sample of dark matter simulated substructure. Substructure models for a range of dark matter particle masses, annihilation channels and halo mass models. are shown. The plotted curves show the σv such that 50% of possible realizations yielded a TS of 25 or greater for at least one dark dwarf.