Constraints to dark matter annihilation from high-latitude HAWC unidentified sources

The $\Lambda$CDM cosmological framework predicts the existence of thousands of subhalos in our own Galaxy not massive enough to retain baryons and become visible. Yet, some of them may shine in gamma rays provided that the dark matter (DM) is made of weakly interacting massive particles (WIMPs), that would self-annihilate and would appear as unidentified gamma-ray sources (unIDs) in gamma-ray catalogs. Indeed, unIDs have proven to be competitive targets for DM searches with gamma rays. In this work, we focus on the three high-latitude ($|b|\geq 10^\circ$) sources present in the 2HWC catalog of the High Altitude Water Cherenkov (HAWC) observatory with no associations at other wavelenghts. Indeed, only one of these sources, 2HWC J1040+308, is found to be above the HAWC detection threshold when considering 760 days of data, a factor 1.5 more exposure time than in the original 2HWC catalog. Other instruments such as Fermi-LAT or VERITAS at lower energies do not detect this source. Also, this unID is reported as spatially extended, making it even more interesting in a DM search context. While waiting for more data that may shed further light on the nature of this source, we set competitive upper limits on the annihilation cross section by comparing this HAWC unID to expectations based on state-of-the-art N-body cosmological simulations of the Galactic subhalo population. We find these constraints to be particularly competitive for heavy WIMPs, i.e., masses above $\sim 25$ (40) TeV in the case of the $b\bar{b}$ ($\tau^+\tau^-$) annihilation channel, reaching velocity-averaged cross section values of $2\cdot10^{-25}$ ($5\cdot10^{-25}$) $cm^3s^{-1}$. Although far from the thermal relic cross section value, the obtained limits are independent and nicely complementary to those from radically different DM analyses and targets, demonstrating again the high potential of this DM search approach.


Introduction
As of today, we believe that about 85% of all the matter in the Universe is of non-baryonic nature [1][2][3]. Despite its large abundance, the ultimate nature of this so-called dark matter (DM) remains unknown, being at present one of the most puzzling questions in modern physics.
N-body cosmological simulations reveal that DM structures form hierarchically in a bottom-up scenario: the DM particles first collapse into small gravitationally bound systems (known as halos), and then form more massive structures through a complex history of merging and accretion. As a result, DM halos contain a very large number of smaller subhalos [4,5].
For DM candidates with weak-scale masses and interactions, as those preferred by supersymmetric particle physics theories [6], subhalos with masses between approximately 10 −11 − 10 −3 M [7][8][9] (depending on the specific DM particle model) up to roughly 10 10 M are expected to exist in a Milky Way-like galaxy. The vast majority of the Galactic DM subhalos are not expected to be massive enough to host baryons and therefore remain completely dark, with no visible counterparts [10][11][12][13][14]. Yet, many of these light subhalos or dark satellites will lie much closer to the Earth compared to the most massive ones, given both their much larger number density and higher survival probability at small Galactocentric radii against tidal disruption [15,16].
If the DM is made of the so-called Weakly Interacting Massive Particles (WIMPs) DM [17,18], these small subhalos may outshine in gamma rays. WIMPs can achieve the correct relic DM abundance (the so-called "WIMP miracle"), through self-annihilation in the early Universe. This process gives rise to a Standard Model (SM) particle-antiparticle pair which, among other possible subsequent by-products, may yield gamma-ray photons [19,20]. Small subhalos may be potentially interesting for this kind of indirect dark matter search since, as mentioned, many of them may be located close enough to Earth to yield high enough gamma-ray fluxes [21][22][23][24]. Interestingly, an important fraction of objects (typically between 30 − 40%, depending on the catalog) in the very high energy (VHE, E > 100 GeV) sky are unidentified sources (unIDs), i.e., objects with no clear single association or counterpart. Some of these unIDs may actually be DM subhalos [22,23,[25][26][27][28][29][30][31][32][33]. In our work, we will explore this possibility focusing on the unIDs recently detected in the VHE regime by the High Altitude Water Cherenkov Observatory (HAWC) [34,35].
One of the greatest challenges when using dark satellites to search for DM is to come up with a reliable characterization of the low-mass subhalo population that would allow for robust predictions of their DM annihilation fluxes. Indeed, there is no N-body cosmological simulation at present able to resolve the entire Galactic subhalo population all the way down to the expected minimum subhalo mass. In this paper, we will adopt the results in [24], where the authors devised an algorithm to repopulate the original Via Lactea II (VL-II) [15] N-body cosmological simulation with low-mass subhalos below its resolution limit, of about ∼ 10 5 M . To do so, they first studied what found for the abundance and distribution of resolved subhalos in the simulation, and extrapolated the relevant quantities to smaller subhalo masses. In addition, they adopted state-of-the-art models to describe the subhalo structural properties [36]. With the predictions of [24] at hand, and in the absence of an obvious DM subhalo candidate among the pool of HAWC unIDs, which is actually reduced to just one source in the end, 2HWC J1040+308, we will show in this work that it is possible to place competitive constraints on the WIMP DM parameter space following the procedure described in [22,23].
The structure of this paper is as follows. In Section 2 we briefly describe the observational status of the VHE sky, paying special attention to the HAWC observatory. We also address our specific target selection in this section. Section 3 is devoted to the computation of the DM constraints taking into account the different instrumental parameters and the theoretical predictions from N-body simulations. We conclude in Section 4, where we also make explicit the caveats in our analysis, compare our results to previous work and make qualitative assessments for future unID-based DM search strategies and experiments.

HAWC unidentified sources
In this section, we briefly review the HAWC observatory and discuss the target selection we made for this work.

HAWC and the TeV sky
The High Altitude Water Cherenkov Observatory (HAWC) is a VHE gamma-ray detector with a one-year sensitivity of 5 − 10% of the flux of the Crab Nebula. Unlike imaging atmospheric Cherenkov telescopes (IACTs), such as H.E.S.S. [37], MAGIC [38] or VERITAS [39], which detect the Cherenkov light emitted by the extensive air showers produced by the gamma rays in the atmosphere, HAWC detects particles of these air showers that reach ground level with water-based Cherenkov detectors. Therefore, HAWC is able to operate continuously and to achieve a field of view (FoV) larger than 1.5 sr [40]. HAWC is located on the flanks of the Sierra Negra volcano near Puebla, Mexico, and consists on a large pond of water tanks, each of 7.3m of diameter, located at 4100 m elevation. The pond contains 900 photomultiplier tubes (PMTs), while each tank contains 200,000 liters of purified water [41]. Starting operations in 2012 with 30 tanks, they were progressively increased to 300 in 2015, detecting TeV gamma and cosmic rays with a high duty cycle (>95%) [42,43].
Currently, there are about 200 known VHE gamma-ray sources detected by a handful of observatories, compiled within the TeVCat 1 catalog [44]. A Galactic map in Hammer-Aitoff projection with all TeVCat sources and their class is plotted in figure 1. From these, roughly 60 sources remain with no clear association. Most of them are probably astrophysical sources of Galactic origin, since many were discovered after the Galactic plane survey by H.E.S.S [45]. Yet, there are several unIDs located at high Galactic latitudes, in particular three of them reported in the 2HWC catalog -the second HAWC catalog [46] comprising 507 days of instrument operation. In this work, we will focus on these sources.

Target selection
We chose to work only with the high-latitude unIDs listed in the 2HWC catalog mainly for two reasons: • Why HAWC? -HAWC is currently the only VHE observatory able to survey a significant fraction of the whole sky. Indeed, the combination of its 1.5 sr FoV, Earth rotation and location in Mexican soil, 19 • above the Equator, makes it possible to observe 2/3 of the sky every day. 2 In contrast, IACTs are pointing telescopes, i.e., they need to know the location of the sources in advance and, thus, just point towards specific targets. For the purposes of this work, it becomes critical to have observations of a large fraction of the sky. This is so because we adopt the methodology introduced in [23] to set constraints on the DM annihilation parameter space using unIDs. The method is based on a comparison between the number of catalogued unID sources and subhalo predictions derived from N-body cosmological simulations. Large fractions of the sky, observed with a nearly uniform exposure, will allow for a more accurate comparison to 1 tevcat2.uchicago.edu 2 In fact, HAWC applies a zenith cut of 45 • , which for a daily exposure translates in 8.4sr, corresponding exactly to 66.85% of the sky.
N-body simulations, as we can derive a more robust statistical determination of the subhalo annihilation fluxes for significantly large sky areas. We remind that these fluxes are proportional to the so-called J-factor: where the first integral is performed over the solid angle of observation (∆Ω), the second one along the line of sight (l.o.s, λ), and ρ DM is the dark matter density profile of a given subhalo.
• Why high latitude? ΛCDM predicts the existence of DM subhalos at all Galactic latitudes distributed nearly isotropically from the Sun's position in the Galaxy. On the other hand, the majority of Galactic VHE astrophysical sources, such as pulsars, binaries, supernova remnants and pulsar wind nebulae, are expected to cluster heavily along the Galactic plane, where most of the stars and gas reside. Since many of these objects are expected to also be hidden among the pool of unIDs awaiting for a proper classification, we expect the distribution of unIDs to peak around zero Galactic latitude as well (as already discussed in Figure 1). For our purposes, these low-latitude sources only add contamination to our sample of potential DM subhalo candidates. Therefore, we apply a cut in our analysis at Galactic latitudes |b| ≤ 10 • . For consistency, this cut will also need to be done on our predicted subhalo distribution. On average, our Galactic cut removes only ∼11% of the simulated subhalos, while 87% of 2HWC unIDs are left out. Actually, in practice we cut a larger region |b| ≤ 30 • in the N-body simulation in order to match the 2/3 sky coverage of HAWC, this way having a totally fair one-to-one comparison among observed and simulated sky.
After this selection, we are left with only three candidates in the 2HWC catalog, namely 2HWC J1309−054, 2HWC J1040+308 and 2HWC J0819+157. All three candidates are reported as very faint, i.e. near the source detection threshold of the catalog. In fact, [35] recently analyzed 760 days of HAWC data, i.e., approximately 50% more integration time than the original 2HWC catalog, and found both 2HWC J0819+157 and 2HWC J1309−054 to be under the detection threshold. This behaviour of source flux versus integration time is not expected for actual sources but from spurious ones instead. Thus, in the following we only consider the remaining candidate for our analysis, 2HWC J1040+308. We note that a spectral analysis was already performed for this source in [35], no preference for DM being found. Yet, more data are probably needed before being able to definitely reject this unID as a DM subhalo. Thus, we conservatively keep this source as the only surviving target in our list of potential subhalo candidates and set DM limits by assuming the existence of one single DM subhalo candidate in HAWC data.

2HWC J1040+308 as a DM subhalo candidate
Before proceeding to set DM limits in the next section, we would like to highlight here some properties of 2HWC J1040+308 that make this unID particularly appealing in a DM context, as such properties may be compatible with a DM origin for the observed emission: • Spatial extension -Spatial extension has been hailed as a "smoking gun" for DM annihilation [24][25][26]47]. Indeed, low-mass yet sufficiently close DM subhalos may be expected to appear as extended unID sources. 2HWC J1040+308 is found to be spatially extended in HAWC data with a radius of 0.5 • [46]. This fact is even more notable for the case of very high latitude sources like this one (b = 61 • ), as sources at high latitudes are typically extragalactic and thus the majority of them are expected to be point-like. Indeed,it would be hard to explain all these features with a single astrophysical source -AGNs would appear as point-like sources, while Galactic sources could appear as extended, yet this source's high latitude would imply a small distance and thus it would be surprising not to have a detection in any other wavelengths. On the other hand, as these unIDs correspond to very faint sources near the detection threshold, it is currently unclear whether they would appear as extended, even if being actual DM subhalos, as the DM annihilation flux profile decreases very rapidly with distance to the subhalo center.
• Multi-wavelength search -As already mentioned, dark satellites are not massive enough to retain baryons and, as a result, they are not expected to shine at any wavelengths due to astrophysical processes. However, gamma-ray emission is expected to happen should these objects be composed of WIMP DM. 3 A dedicated search at different wavelengths was performed for 2HWC J1040+308, with null results. In particular, a combined search in less energetic gamma rays with Fermi-LAT and VERITAS was performed in [52], with no detection. An additional search can be performed with the SSDC online tool, 4 where no significant emission at lower energies is found.
• Heavy WIMP mass -It is interesting to note that a joint analysis of Fermi LAT [53], VERITAS [54] and HAWC data was recently done in [52] and no gamma-ray emission was reported for 2HWC J1040+308 in the comparatively lower energy range covered by the LAT and VERITAS. This is so despite the fact that 2HWC J1040+308 exhibits a hard spectrum with a photon spectral index of Γ = 2.08 ± 0.25 in the HAWC energy range [46], which would in principle make a detection at low energies easier to realize. If interpreted in a DM scenario, these results suggest heavy, TeV WIMP masses as we would not expect a detection by Fermi LAT or VERITAS only in case of significantly high values of the WIMP mass, for which a DM spectrum beyond the range of sensitivity of these two instruments would be generated. Also, interestingly, this unID slightly prefers a fit to a "Cutoff Power Law" instead of a "Simple Power Law" [35], which is a parametric form that better reproduces a typical DM annihilation spectrum [22,23]. Heavy WIMPs are well motivated [55][56][57] and, probably, favoured in the light of current DM constraints in the usual σv (velocity-averaged annihilation cross section) vs. m χ (WIMP mass) parameter space: IACTs are still far from being able to probe the thermal relic cross section values [58] for large, TeV WIMP masses, while for much lighter, O(GeV) masses there are already robust and strong constraints, e.g., [23,59].
All the above considerations make 2HWC J1040+308 an excellent candidate for being a DM subhalo. Yet, lacking a definitive answer for the moment we will proceed and will set constraints to DM annihilation using HAWC unIDs, just assuming that we cannot discard 2HWC J1040+308 as a DM subhalo.

DM constraints
The methodology we follow to set our DM constraints is explained in full detail in [23]. In short, these are computed according to i) the number of remaining HAWC unIDs as potential DM subhalo candidates, i.e. just one source (2HWC J1040+308), as discussed in Section 2; ii) the flux sensitivity of the instrument to DM annihilation; and iii) the J-factor predictions for the Galactic DM subhalo population as derived from N-body cosmological simulations. The annihilation cross section and WIMP mass σv and m χ are related by [23]: 3 DM-induced emission may also be expected in radio due to synchrotron radiation from secondary particles [48][49][50][51]; yet the high density of radio sources would make any potential association probably very challenging.
where F min is the minimum flux needed for a point-source detection, J is the astrophysical J-factor (see Eq. 1), and ∞ E th dN dE dE = N γ is the integrated DM spectrum. This DM spectrum is obtained from the PPPC4 DM ID tables [60], including electroweak corrections, and integrated from E th = 300 GeV for a variety of annihilation channels. Hereafter, though, and for the sake of clarity, we will focus on bb and τ + τ − channels, considering only pure annihilations (unity branching ratio).

Minimum detection flux
As mentioned previously, the minimum detection flux (F min ) is defined as the one required by the instrument to have a detection. More specifically, it is the source flux that is needed to reach TS = 25 over the background, where TS is the Test Statistic, defined as, L(H 1 ) and L(H 0 ) are, respectively, the likelihoods under the source and no source (null hypothesis) assumptions. Ideally, as done in [23], one would need to account for the variations of this quantity across different instrumental setups (integration time, energy threshold, instrument response functions, etc.) and adopting different DM spectra. 5 Finally, the source latitude should also be taken into account, however this is a second-order effect once the Galactic plane is removed from the analysis, as it is the case. Unfortunately, the VHE sky still lacks a proper model for the all-sky Galactic diffuse gamma-ray emission, which exists at lower energies [61]. This, and the fact that we do not possess the HAWC instrumental response functions (IRFs), private for this Collaboration's internal use at the moment, precludes a precise characterization of the HAWC F min . Instead, we will have to rely only on the differential sensitivity from figure 16 of [62]. 6 At each energy, the F min can be obtained by dividing the differential sensitivity (in TeV · cm −2 · s −1 ) value by the log-central energy value in each considered bin (to obtain ph · cm −2 · s −1 ). In our case, we need the F min corresponding to a declination of 30.87 • , which turns out to be very near the absolute minimum of the curve shown in figure 16 of [62]. 7 Note, however, that the reported flux in [62] assumes a power law spectrum with a photon spectral index Γ = 2.63 (Crab-like), while our DM subhalo candidate exhibits a considerably harder spectrum, Γ = 2.08, and therefore it should be comparatively easier to detect. Also, the sensitivity is expected to scale as ∼ √ t (assuming Poisson photon statistics), where t is the exposure time. Thus, as the reported flux in [46] was computed for 507 days of operation, the performance of HAWC would be underestimated for the 760 days of the current analysis.
With these limitations in mind, it is actually possible to compute a better estimate for F min than the one described above. We will define the improvement factor as the ratio between our estimation of F min and the official one. First, we take into account the improvement due to the photon spectral index -a source with smaller index (harder spectrum), will be easier to detect than a source with larger index (softer spectrum). We can use the scaling relation reported in [46], i.e., the variation of F min with the source spectral index, for an energy value of 7 TeV, used as the pivot energy. By calibrating the improvement at this energy, we can extend it to all the considered energy range (300 GeV -100 TeV), taking into account the spectral shapes of both power laws. Then, we rescale the obtained F min by differences in exposure times, i.e., a factor √ 507/760 = 0.82. 5 This is so because, by default, this quantity is computed for a source spectrum with a power law index Γ ∼ 2, while the DM is better parametrized by a "Power Law with SuperExponential Cutoff" [23], where the index and cutoff energy vary over the annihilation channel and WIMP mass. 6 This sensitivity is computed for the Crab, which is located at DEC∼ 22 • , while our source is at ∼ 30 • instead. Fortunately, the HAWC sensitivity for both locations is expected to be very similar [34].
Also, as our unID is an extended source, with a reported extension of 0.5 • , the HAWC sensitivity is expected to worsen with respect to that for a point-like source. Indeed, according to [34], a 0.5 • source would require a flux 1.93 times larger to be detected at the same significance, independently of the assumed spectral index. As we are tailoring the sensitivity to our candidate, we will take this correction factor into account. As a result, the disagreement between our results and those in [35] is reduced to a factor ∼ 25.
A comparison between both the originally reported F min in [62] (Γ = 2.63, t = 507 days, point-like) and the "improved" F min here described for our candidate (Γ = 2.08, t = 760 days, 0.5 • ) is shown in Figure 2 as the blue and red curves, respectively. The latter is the one that will be used to set our DM constraints in Section 3.3. 8  In blue, the F min as originally reported in [62] for the case of a point source described by a power law spectrum with index Γ = 2.63 and 507 days of exposure. The red curve shows the F min for our DM subhalo candidate instead (with a spatial extension of 0.5 • , Γ = 2.08 and 760 days of exposure). This one is computed taking into account i) the improvement factor of ∼ 5 derived for this unID at a pivot energy of 7 TeV [46], marked in the figure as a dashed grey vertical line; ii) the difference in exposure times, which makes the F min improve by a factor √ 507/760 = 0.82; and iii) a factor ∼2 worsening due to the extension of the source, according to [34]. See text for further details on these computations. The overall F min improvement is roughly an order of magnitude at low energies, while at the highest energies there is no improvement at all.

J-factor
The last ingredient needed to set the DM constraints is the J-factor. In the case of unID-based DM searches, we conservatively assume that all N DM subhalo candidates in our unID list (N = 1 in this work) are in fact DM subhalos [22,23]. Then, we require consistency of this supposedly observed number of DM subhalos with that obtained from N-body simulations, by assuming that the J-factor of the N th (1 st , in this work) most brilliant DM subhalo in our simulation sets the border between the population of detected and non-detectable subhalos. 8 The spectral index is reported to have an uncertainty Γ = 2.08 ± 0.25. We checked that this uncertainty translates into a factor ∼4 uncertainty in F min at most (factor 2 when compared to the benchmark Γ = 2.08).
As input, we use the N-body simulation work by [24], specifically designed to assess the relevance of low-mass subhalos for this kind of studies. The authors of [24] repopulate the original Via Lactea II DM-only N-body cosmological simulation of a Milky-Way-size halo [15] with subhalos well below its resolution limit by applying bootstrapping techniques and semi-analytical extrapolations of the relevant quantities. In particular, the subhalo mass function is extended down to 10 3 M , the subhalo distribution within the host halo is assumed to be similar to that of the resolved subhalos in the parent simulation, and a state-of-the-art concentration model [36] is used to model the subhalo structural properties and, ultimately, to compute the J-factors. Actually, hundreds of realizations of the parent simulation are performed that allow to derive statistical meaningful J-factor results. More precisely, as done in [23], in order to derive 95% C.L. upper limits on the DM annihilation cross section, we obtain the distribution of J-factor values for the N th subhalo under consideration and then as reference J-factor the one above which 95% of the J-factor distribution is contained; see Figure 3.
In this particular case, we cut out from the simulation results the Galactic region |b| ≤ 10 • in order to match the cut that we applied on the data to avoid the Galactic plane contamination (see Section 2.2). We also take into account that HAWC, unlike Fermi-LAT, does not achieve a uniform exposure for the whole sky but only for ∼ 2/3 of it. Specifically, to mimic this effect in the simulation side, we cut out the |b| ≤ 30 • region (i.e., we only keep 2/3 of the whole sky) in every realization. Finally, as we are interested in dark satellites alone, only M≤ 10 7 M subhalos (i.e., not massive enough to retain baryons) are considered. We checked that the resulting distribution looks isotropic as seen from the Earth. The result of applying and of not applying cuts on the J-factor distributions for the most brilliant subhalo in the simulations is shown in Figure 3. The value used to set the constraints, above which 95% of the distribution is contained, is log 10 (J) = 19.47, a factor 4.26 smaller than the value before cuts. 9 9 The J-factor of our brightest subhalo is comparable to the one typically quoted for dwarf spheroidal satellite galaxies (see e.g. [63]), which are expected to be well above the extragalactic isotropic DM-induced gamma ray background. On the other hand, the value of the DM-induced Galactic diffuse emission is comparable to the isotropic contribution at high latitudes [64], and therefore also negligible compared to the brightest subhalo J-factor. Finally, the boost due to Galactic unresolved substructure contributing to the diffuse emission in the line of sight of this unID is expected to be at the level of few percent and, thus, not relevant here. The boost due to substructures is only particularly important when integrating the subhalo signal for the whole host halo; see e.g. the discussions and results in [36].  [23] from the Via Lactea II N-body simulation. Red and blue histograms show, respectively, the values before and after applying our selection cuts on the simulated DM subhalos, namely a mass cut (M ≤ 10 7 M ) and a coordinate cut (|b| > 30 • ). The dashed, purple vertical line marks the value of the J-factor that will be used to set the 95% C.L. DM limits in Section 3, and it is defined as the one above which 95% of the J-factor distribution is contained.

Results
Once all involved quantities in Equation 2 have been properly characterized, we can set constraints on the DM parameter space at 95% confidence level. These are plotted in figure 4 for the bb and τ + τ − annihilation channels. We remind, once again, that these limits adopt a single HAWC unID as being a potential DM subhalo (2HWC J1040+308).
The obtained 95% upper limits to the DM annihilation cross section are most stringent for masses of ∼20 (2) TeV for the bb (τ + τ − ) annihilation channel, reaching cross section values around 5 · 10 −25 (8 · 10 −25 ) cm 3 s −1 . This is roughly one order of magnitude far from the thermal relic cross section value, yet they are competitive with today's IACTs best constraints (derived from H.E.S.S. observations of the Galactic center [65]) for very heavy WIMPs annihilating to bb above ∼50 TeV. It is worth emphasizing that our DM constraints are independent and complementary to the ones obtained for other astrophysical targets and by means of different analysis techniques.
We must note here that we observe a mismatch between the DM constraints obtained in this work and the ones shown in [35]. This disagreement, that reaches a factor ∼ 25 for τ + τ − , is hard to understand in detail but it is probably due to several factors. In particular, our F min is expected to be a fair approximation while [35] adopts the actual sensitivity of HAWC to a DM spectrum (only computable with the IRFs). In that work, the structural properties of subhalos were described according to [66], which is well-suited for field halos, not subhalos. Indeed, we use [36], where subhalos exhibit concentration values a factor 2-3 larger than field halos of the same mass. As the J-factor is roughly proportional to the third power of the concentration, this correction would thus yield J-factors about a factor 8-9 times larger, leading to better constraints. Other possible systematics, such as those coming from the Galactic diffuse emission model and the observation strategy adopted in [35], may turn out to be particularly relevant. Also, we note that they show cross-section values needed for detection,  Figure 4. 95% upper limits on the DM annihilation cross section as derived from HAWC unIDs and predictions from N-body cosmological simulations. Indeed, these constraints use 2HWC J1040+308 alone, i.e., the only HAWC unID located at |b| ≥ 10 • and surviving our selection criteria in Section 2.2. Left (right) panel shows the 95% C.L. upper limits for the bb (τ + τ − ) annihilation channel. For comparison, we also show as a blue, dot-dashed line the DM constraints obtained from Fermi-LAT unIDs using the same methodology [47]. Limits from the observation of the Galactic center region by H.E.S.S., i.e., the best DM constraints achieved from IACT observations at present, are included in both panels as well as a green, dashed line [65]. Finally, a dotted, orange line shows the constraints from Fermi-LAT observations of dSphs [59]. Note that we did not include in this Figure the results in [35], as a one-to-one comparison would be misleading given the fact that the methodology and underlying assumptions in both works are significantly different; see discussion in Section 3 for details.
while we here show 95% upper limits. All in all, a precise one-to-one comparison with the results in [35] is not possible, as the methodology and underlying assumptions are different.
For the sake of comparison, we also show in figure 4 the constraints found by [23] following a similar methodology, but using unIDs in Fermi-LAT catalogs instead. Both these low and high energy unID-based DM limits meet at roughly 2 (0.7) TeV for bb (τ + τ − ) annihilation channel.
Finally, it is important to note that our DM limits are subject to potentially large uncertainties. In addition to the already mentioned caveats in the computation of the minimum detection flux (see Section 3.1), 10 there exist important theoretical uncertainties in the N-body simulation predictions, such as the survival probability of the lightest subhalos near the Galactic center [67][68][69], 11 the precise subhalo structural properties, the impact of baryons on the subhalo population [22,70] or the value of the minimum subhalo mass that is adopted to separate between visible and dark satellites [71]. These and other issues that affect the Galactic subhalo population will be addressed elsewhere. Also, for WIMP masses above ∼ 10 TeV, the theoretical spectra of [60] should include higher-order EW corrections, although these are expected to be potentially relevant only for leptonic channels such as e + e − and µ + µ − .

Summary and conclusions
In this work, we have studied the HAWC unIDs reported in the 2HWC catalog to search for potential DM subhalos. HAWC proves to be the best VHE gamma-ray observatory to perform this kind of search at present, as it is the only one that surveys a significant portion of the sky in this energy window.
Starting from a total of 23 unIDs in the 2HWC, we apply several selection criteria in Section 2.2 based on expected DM signal and subhalo properties. More precisely, at first we only select those HAWC unIDs located at Galactic latitudes |b| ≥ 10 • to avoid contamination from astrophysical sources along the Galactic plane. 87% of the 2HWC unIDs do not pass this cut, i.e. we are left with only three unIDs. We then consider these three unIDs as DM subhalo candidates, yet two of them do not even reach the flux detection threshold when adding ∼250 more days of data [35]. Therefore, only one unID, labeled as 2HWC J1040+308, remains as a potential DM subhalo after our filtering procedure. Indeed, the measured flux of this unID grows as expected when considering more integration time, which reinforces its actual existence.
When scrutinized in further detail, 2HWC J1040+308 turns out to be specially interesting in a DM context, as it is reported in the 2HWC catalog as spatially extended, it exhibits a hard spectrum and does not possess visible counterparts neither at other wavelengths nor for other gamma-ray observatories operating at lower energies such as Fermi-LAT and VERITAS. The combination of the latter with the hard source spectrum might point towards particularly heavy WIMP masses, if the gamma-ray emission was actually produced by DM annihilation.
While waiting for more data that can help to shed further light on the true nature of 2HWC J1040+308, we proceed in Section 3 and set limits on the DM annihilation cross section by conservatively assuming that this source is actually a subhalo. To do so, we follow the methodology in [23] and compare the HAWC unID observations with predictions for the Galactic subhalo population as derived from N-body cosmological simulations of a galaxy like our own. In particular, we follow the simulation work in [23,24], where the VL-II DM-only simulation was repopulated to include subhalos with masses down to 10 3 M . This work provides the expected subhalo J-factors, which we combine with our own estimate of the HAWC instrumental sensitivity to DM subhalos in order to set 95% C.L. upper limits in the σv -m χ parameter space.
Our DM limits are most stringent for masses of ∼20 (2) TeV for the bb (τ + τ − ) annihilation channel, reaching cross section values around 5 · 10 −25 (8 · 10 −25 ) cm 3 s −1 . Although roughly one order of magnitude far from the thermal relic cross section value, our limits are competitive with today's IACTs best constraints [65]) above ∼ 50 TeV for bb annihilation channel. Also, these constraints are independent and complementary, yet competitive, to those obtained for other astrophysical targets adopting diverse analysis techniques.
The analysis presented in this work is subject to some potentially important uncertainties. In particular, a proper determination of the minimum HAWC detection flux taking into account the specific instrumental setup used in the observations, as well as the DM spectrum, annihilation channel and WIMP mass, would be required. The HAWC IRFs are nevertheless not public, thus this study is beyond the scope of the current paper. Instead, we used public information in [46,62] to come up with a fair estimate of the HAWC sensitivity to DM subhalos (Section 3.1). Another source of uncertainty is the lack of a proper Galactic diffuse gamma-ray background model, which prevents us from including it in our estimate of the HAWC sensitivity to DM. Yet, as in the end only one unID survives our cuts, this single unID is located at high Galactic latitudes (where the diffuse emission is expected to be subdominant) and exhibits a hard spectrum well fitted by a power law, we expect this uncertainty from the Galactic diffuse model to be a second order effect. Finally, there are also theoretical uncertainties hidden behind the adopted N-body simulation work. In particular, the survival probability of the smallest subhalos near the Galactic center -thus subject to strong tidal forces -is a matter of fierce debate in the community at present, with different works [67-69] providing diverse answers. In addition, the impact of baryons on the subhalo population was not considered here. Baryons may induce a reduction of the number of subhalos near the Sun's position, which would worsen our DM constraints probably by a factor of a few [22,70].
Looking into the future, more data, preferably not only in the VHE regime, are needed in order to elucidate the true nature of 2HWC J1040+308, i.e., the only HAWC unID that survives our selection cuts and, indeed, a promising DM subhalo candidate. Related to this, we note that the simplest and most robust way to improve our DM constraints is to also associate 2HWC J1040+308 to a known astrophysical source. In fact, a simple detection of variability or multi-wavelength emission would be enough to discard it as a DM subhalo (see several examples of unID rejections in [23]). As HAWC accumulates exposure time, also its minimum detectable flux will be lowered, thus allowing for better DM limits for a fixed number of unIDs. New HAWC observations and additional exposure time will surely bring up new unIDs though, thus in the end the resulting DM limits may improve or worsen depending on the number of unIDs that will be potential DM subhalos [22,23].
Finally, the upcoming Cherenkov Telescope Array (CTA) [72], with its superb instrumental capabilities, enhanced sensitivity to DM, improved angular resolution to pinpoint source extension, and already planned surveys of large portions of the VHE sky with unprecedented sensitivity, should be able to set the most competitive DM limits in the TeV energy range by means of the unID-based method used in this work. Alternatively, it is also possible that CTA may bring the first robust DM subhalo discovery.