Searching for TeV Dark Matter in Irregular Dwarf Galaxies with HAWC Observatory

While we have clear evidence of the existence of dark matter (DM) in the universe, we do not know the nature of DM, and various candidates have been proposed. One of the most common DM candidates is the weakly interacting massive particles (WIMPs), a family of particles with interaction at the weak scale that can be thermally produced in the early universe. More details can be found in Bertone et al. (2005), Bertone & Hooper (2018), and Arbey & Mahmoudi (2021). We assume that the generic WIMP (as the neutralino predicted by SUSY) is the only component of DM. From direct searches and collider experiments (see Mitsou 2015; Schumann 2019; Arbey & Mahmoudi 2021 and references therein), very restrictive constraints have been established in the mass range of 5 GeV to 1 TeV. Indirect searches explore a broad range of masses, from 10 GeV up to the unitarity limit (400 TeV for a Majorana particle; Enqvist & Kainulainen 1991; or 144 TeV considering the effect of bound states; Smirnov & Beacom 2019) by assuming that the DM particles can annihilate or decay to Standard Model (SM) particles producing stable final states as protons, neutrinos, and gamma rays. Then, indirect DM searches allow us to explore candidates with masses (above 10 TeV) that are not accessible to LHC or direct detection experiments. Indeed, by using observational data from wide-field gamma-ray observatories we can constrain the range of fundamental parameters of DM candidates. To do this, we use DM-dominated astrophysical objects where the probability to observe a DM signal is high. The target with the highest expected DM signal in the gamma-ray sky is the Galactic Center, but the analysis of this region involves the understanding of the distribution and amount of astrophysical gamma-ray sources in the region around the center. We can also observe regions where we expect negligible gamma-ray emission from astrophysical processes, which would lead to a clean DM signal. Examples of these objects are not only very well-known dwarf spheroidal (dSph), but also dwarf irregular (dIrr) galaxies. The dIrr galaxy's population is characterized by very low star formation rates (SFRs; McGaugh et al. 2017) and the lack of massive stars (Schulte-Ladbeck & Hopp 1998; Dunn 2007). Hence, there is little background gamma-ray emission at energies above 1 TeV and dIrr can be considered as essentially background-free objects in indirect DM searches (see, e.g., Gammaldi et al. 2018). Here, we use a sample of dIrr galaxies within the High Altitude Water Cherenkov (HAWC) field of view to show the HAWC sensitivity to look for DM signatures in these objects, compared with the limit obtained by benchmark targets (Albert et al. 2018). This paper is organized as follows: in Section 2 we briefly introduce the HAWC Observatory. In Section 3, we present the sample of dIrr galaxies used for this study and discuss the DM density profiles computed for dIrrs. In Section 4, we discuss the data set and the analysis for individual sources, Section 4.1, and the combined analysis for this population, Section 4.2. In Section 5, we show the exclusion limits computed for 31 dIrr galaxies and the comparison with limits from dSph galaxies. Finally, we conclude in Section 6.

Located in Sierra Negra, México at an altitude of 4100 m, the HAWC Observatory is an extended array of 300 Water Cherenkov Detectors (WCDs) to detect air showers produced by very high energy gamma rays and cosmic rays. Every WCD is 7.3 m in diameter and 4.5 m deep, filled with 200,000 l of purified water, and has instruments with four photo multiplier tubes (PMTs) anchored to the bottom. The PMTs collect the Cherenkov light produced by charged particles passing through the WCDs. One of the PMTs is located in the middle, and the other three are located at a distance of 1.8 m from the central PMT and equally spaced (120°). The HAWC Observatory is sensitive to gamma rays with energies in the range between 1 and 100 TeV, with a wide field of view, covering 2/3 of the sky each day and a duty cycle >95%. This allows the HAWC Observatory to study large populations of astrophysical sources and constrain parameters that are common to all targets in the sample.

When an air shower reaches the detector, hits recorded during the air shower event are selected within a time window of 250 ns around the trigger time. Hits that survive the time selection and quality cuts are used for the reconstruction of the air shower parameters. After the reconstruction, analysis of air shower events is based on cuts depending on the energy and size of events measured by the detector. Then, HAWC events are organized into bins according to the size of the area covered by each event recorded on the ground. The event size is defined as the ratio of the number of PMT hits used for the reconstruction to the total number of available PMTs for the reconstruction. A range of values of this ratio is called a fractional bin. Table 1 shows the definition of the nine fractional bins, f_hit, used in HAWC. This definition of fractional bins according to the size of the air shower recorded in the HAWC Observatory is weakly correlated to the energy of the primary particle (Abeysekara et al. 2017a, 2019). For more details about the reconstruction and estimation of physical parameters of the air showers in HAWC see Abeysekara et al. (2017a, 2017b) and Albert et al. (2020).

Figure 1 shows the sky map and position of our sample of 31 dIrrs in the HAWC field of view. All the dIrr galaxies in the sample have a distance to Earth smaller than 11 Mpc. All of these galaxies are extragalactic sources, with some of them belonging to the Local Group of galaxies having redshifts z < 0.001. The sample is taken from Karukes & Salucci (2016), where the authors computed the DM properties of 36 dIrr galaxies using kinematical data, and after selecting the galaxies within the HAWC field of view, the sample is reduced to 31 dIrr galaxies. They establish that these galaxies follow a low-mass version of the universal rotation curve (URC) function. The URC is a generic function of radial distance to the center of a spiral galaxy that describes the rotation curve of all galaxies in the local volume, and parameterizes the distributions of matter inside a galaxy (Persic et al. 1996). These parameters are found by fitting three contributions to the rotational motion in the galaxy: a stellar disk, a H i disk, and a spherical DM halo. In the case of a DM halo, the profile with the best fit is constrained to be a Burkert density profile (Burkert 1995):

$$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(r)={\rho }_{0}\displaystyle \frac{{r}_{0}^{3}}{({r}_{0}+r)({r}_{0}^{2}+{r}^{2})},\end{eqnarray} \tag{ 1 }$$

where r₀ is the core radius, and ρ₀ is the normalization density. These observational core profiles represent a modified, isothermal sphere that better fits the observations of the rotation curves that cannot be explained by DM-cusped profiles, resulting from N-body cosmological simulations (Navarro et al. 1997; Gammaldi et al. 2021). The ρ₀ and r₀ parameters used in this study are taken from Karukes & Salucci (2016). Table 2 shows the sample of galaxies, their coordinates, virial size, and the value of astrophysical factors computed for these galaxies.

Zoom In Zoom Out Reset image size

Table 2. Sample of dIrr Galaxies

Name	R.A.	Decl.	θvir	${\mathrm{log}}_{10}(\tfrac{{\rm{J}}}{{\mathrm{TeV}}^{2}\,{\mathrm{cm}}^{5}})$	${\mathrm{log}}_{10}(\tfrac{{\rm{D}}}{\mathrm{TeV}\,{\mathrm{cm}}^{2}})$
	(deg)	(deg)	(deg)
And IV	10.62	40.57	0.326	9.764	13.463
DDO 101	177.91	31.51	0.309	10.356	14.312
DDO 125	186.92	43.493	0.533	10.467	14.165
DDO 133	188.22	31.54	0.784	11.501	15.274
DDO 154	193.52	27.15	0.695	11.800	15.397
DDO 168	198.61	45.91	1.142	11.365	15.271
DDO 43	112.07	40.77	0.472	10.109	13.853
DDO 52	127.11	41.85	0.570	10.452	14.401
Haro 29	186.56	48.49	0.333	9.974	13.764
Haro 36	191.73	51.61	0.550	10.642	13.581
IC 10	5.10	59.29	3.857	11.857	15.619
IC 1613	16.19	2.13	2.361	11.632	15.325
NGC 3741	174.02	45.28	0.405	9.814	13.417
NGC 6822	296.23	−14.80	5.325	12.173	15.943
UGC 11583	307.56	60.44	0.853	10.676	14.605
UGC 1281	27.38	32.59	0.952	10.854	14.739
UGC 1501	30.31	28.84	1.032	10.937	14.843
UGC 2455	194.92	25.23	0.569	10.392	14.250
UGC 5272	147.59	31.48	0.765	10.721	14.731
UGC 5427	151.17	29.36	0.436	10.133	14.007
UGC 5918	162.40	65.53	0.612	10.512	14.420
UGC 7047	181.01	52.58	0.649	10.630	14.444
UGC 7232	183.43	36.63	0.654	10.845	14.581
UGC 7559	186.77	37.14	0.700	11.105	14.938
UGC 7603	187.18	22.82	0.652	11.368	15.251
UGC 7861	190.46	41.27	0.535	10.804	14.715
UGC 7866	190.56	38.50	0.496	10.672	14.462
UGC 7916	191.10	34.38	0.489	11.012	14.898
UGC 8508	202.68	54.91	0.584	10.362	14.071
UGC 8837	208.68	53.90	0.737	10.856	14.802
WLM	0.49	−15.46	2.609	12.062	15.777

Note. We show the 31 dIrr galaxies within the HAWC field of view used in this study. Columns: name of the galaxy (1), the R.A. (α) (2), and decl. (δ) (3) of the galaxy, the angular extension of the virial radius (4), the astrophysical factor for annihilation (5), and decay (6) computed with Clumpy (Hütten et al. 2019).

Download table as:  ASCIITypeset image

3.1. DM Photon Flux

In this section, we consider both the annihilation and decay of WIMPs. In the case of annihilation, the differential flux of the gamma-ray photons is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\rm{\Phi }}}_{\gamma }^{\mathrm{ann}}}{{dE}} & = & \displaystyle \frac{\langle {\sigma }_{\chi }v{\rangle }^{\mathrm{ref}}}{8\pi {m}_{\chi }^{2}}\displaystyle \sum _{f}{B}_{f}\displaystyle \frac{{{dN}}_{f}^{\mathrm{ann}}}{{dE}}\\ & & \times \,{\displaystyle \int }_{{\rm{\Delta }}{\rm{\Omega }}}{\displaystyle \int }_{{\rm{l}}.{\rm{o}}.{\rm{s}}.}\,{dld}{\rm{\Omega }}\,{\rho }_{\mathrm{DM}}{\left(r(l)\right)}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${{dN}}_{f}^{\mathrm{ann}}/{dE}$ is the differential spectrum of photons produced for annihilation channel f, m_χ is the mass of the WIMP, and 〈σ_χ v〉^ref is the reference thermal averaged annihilation cross section. The sum is over the total number of channels with branching ratios B_f that contribute to the production of photons (∑_f B_f = 1). We assume that branching ratios have a value of 1 for a specified channel, while the others have B_f = 0. The term in the integral is called the astrophysical factor J, or the J factor. The J factor is the double integral of the DM density profile squared along the line of sight l and over the solid angle ΔΩ around the line of sight.

For decaying DM, the differential flux of photons is computed by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\rm{\Phi }}}_{\gamma }^{\mathrm{dec}}}{{dE}} & = & \displaystyle \frac{1}{4\pi {m}_{\chi }{\tau }_{\chi }^{\mathrm{ref}}}\displaystyle \sum _{f}{B}_{f}\displaystyle \frac{{{dN}}_{f}^{\mathrm{dec}}}{{dE}}\\ & & \times \,{\displaystyle \int }_{{\rm{\Delta }}{\rm{\Omega }}}{\displaystyle \int }_{{\rm{l}}.{\rm{o}}.{\rm{s}}.}\,{dld}{\rm{\Omega }}\,{\rho }_{\mathrm{DM}}\left(r(l)\right).\end{array}\end{eqnarray} \tag{ 3 }$$

As in the case of annihilation, ${{dN}}_{f}^{\mathrm{dec}}/{dE}$ is the differential spectrum of photons produced for a decay channel f, and ${\tau }_{\chi }^{\mathrm{ref}}$ is the reference lifetime of the DM candidate. For DM decay, the astrophysical factor D, or D factor, is the double integral of the DM density profile ρ_DM along the line of sight l and over the solid angle ΔΩ around the line of sight. We show the sample of galaxies and their values of J and D factors in Table 2.

3.1.1. Photon Spectra

The production of photons from the annihilation (decay) of DM particles is due to the decay or hadronization processes of the unstable products. In both cases, the spectrum of photons is continuous and has an energy cutoff at the energy available in the process, the mass (half mass) of the DM particle for annihilation (decay). For this work, we considered WIMP masses in the range of 1–100 TeV and annihilation (decay) to five channels: b and t quarks, μ and τ leptons, and the W boson. The spectrum of photons is calculated with Pythia 8 (Sjöstrand et al. 2015). Figure 2 shows the spectrum of photons for the annihilation of a WIMP with mass m_χ = 60 TeV.

Zoom In Zoom Out Reset image size

3.1.2. Astrophysical Factor

For the set of dIrr Galaxies we studied here, the structural properties of luminous and DM contributions are constrained using kinematical data taken from Karukes & Salucci (2016). The DM density is described by a Burkert profile (Burkert 1995). The Burkert profile is a density distribution that resembles an isothermal profile in the inner regions (r < r₀) and a distribution with slope −3 in the outer regions. For this study, the virial radius, R_vir, and other DM-related parameters were computed assuming an overdensity Δ equal to 200, and an M_DM–c_vir taken from Sánchez-Conde & Prada (2014). We do not take into account the presence of DM subhalos as estimates of the boost factor may increase the total gamma-ray flux up to ∼10× the contribution of the main halo. As we assume a point-source model for the DM-induced gamma-ray emission, the total boost factor is used as a multiplicative factor to the exclusion limits. However, this boost factor does not change the conclusions of the analysis. Galaxies with θ_vir > 1° are left for a future study, and we assume that the spatial emission from the galaxies is coming from a point source as the angular resolution for bins 1 and 2 is 0.98° and 0.71°, respectively, see Table 2. In order to compute the integrals for the J and D factors, we use the Clumpy package (Hütten et al. 2019). We compute the astrophysical factors over the total extension of the DM halo.

In Table 2, we report the angular virial radius θ_vir (Column 4), and the J and D factors (Columns 5 and 6) for dIrrs. According to these values, we observe that the best targets in the sample are the galaxies NGC 6822, IC 10, IC 1613, WLM, and DDO 154 because of their large J and D factors. In particular, the galaxy DDO 154 is located in a decl. band that is favorable to the HAWC Observatory, see Figure 1. Furthermore, the values of D and J factors in our work are comparable to the values reported in Gammaldi et al. (2018, 2021). In particular, the slight difference with respect to the values in Gammaldi et al. (2018) is because the authors assume a value of overdensity, Δ, equal to 100. Differences with respect to the values reported in Gammaldi et al. (2021) are due to the reanalysis of available kinematic data to obtain the parameters of the DM profile. We also note that the best targets in the sample reported here are consistent with the results obtained in Gammaldi et al. (2018, 2021).

Additionally, we can compare the astrophysical factors to the values for dSph galaxies, in particular with galaxies within the field of view of the HAWC Observatory. We used the sample of galaxies studied in Albert et al. (2018). As we describe later (see Section 5.3), the population of dSph can be divided into two subclasses according to the number of stars hosted in the galaxies: classical and ultrafaint galaxies. The values of the J and D factors for dSph galaxies are shown in Table 3. We can observe that values of J factors for dIrr galaxies are roughly between 2 and 3 orders of magnitude smaller than the values for ultrafaint galaxies, while this difference is reduced for classical dSph galaxies. This may indicate that dIrr galaxies are not very suited to perform searches of annihilating DM. However, one should consider that values of J factors for ultrafaint dSphs, in the majority of cases, have larger uncertainties due to the lack of stellar data to constrain the DM profile. The case for decay is more interesting, as we observe that the value of D factors for the three different populations are very similar, so we will consider combined analysis between dIrr and dSph populations for a future study.

Table 3. Values of J and D Factors for dSph Galaxies Observed with HAWC

Name	${\mathrm{log}}_{10}(\tfrac{{\rm{J}}}{{\mathrm{TeV}}^{2}\,{\mathrm{cm}}^{5}})$	${\mathrm{log}}_{10}(\tfrac{{\rm{D}}}{\mathrm{TeV}\,{\mathrm{cm}}^{2}})$
Ultrafaint
Triangulum II	14.44	15.42
Segue I	13.66	15.64
Coma Berenices	13.32	15.71
Classical
Draco	13.37	16.15
Leo I	11.57	15.04
Leo II	12.11	14.33
Ursa Minor	13.24	14.92

Note. Values are taken from Albert et al. (2018). We convert the values to the same units for J and D factors for the dIrr galaxies reported in Table 2.

Download table as:  ASCIITypeset image

3.2. Astrophysical Gamma-Ray Emission

Galaxies that present star formation regions have been reported as gamma-ray emitters (Ohm et al. 2010). Because of this, it is necessary to consider this emission to compute the photon fluxes. However, the galaxies we use in this study are characterized by very low SFRs (Schulte-Ladbeck & Hopp 1998; Dunn 2007; McGaugh et al. 2017). Therefore, gamma-ray emission to teraelectronvolt energies is considered negligible. We use the method described in McGaugh et al. (2017) to estimate the SFR of the galaxies in our sample. Then, following the method described in Martin (2014), we estimate the gamma-ray flux of the dIrr galaxies to be ≲7 orders of magnitude below the Crab flux (see Figure 3). The energy-integrated sensitivity to a point source with a spectral index of α = −2.5 in the position of the Crab Nebula is 0.028 Crab Units (CU), see Figure 2 in Albert et al. (2020), so we do not consider the SFR-induced gamma-ray emission in our gamma-ray model.

Zoom In Zoom Out Reset image size

To derive the exclusion limits, we use data from the HAWC Observatory comprising 1017 days of data. We select this period of time to be able to compare it with the results for dSph galaxies in Albert et al. (2018).

To test our DM model, we use the maximum likelihood method to estimate best-fit values and exclusion limits for the parameters of interest. This method constrains the values of free parameters by maximizing the likelihood function given a data set D. For the HAWC Observatory, we used the binned log-likelihood function given by

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}L({S}_{i,k},{B}_{i,k}| D)=\displaystyle \sum _{i,k}\left[{N}_{i,k}\mathrm{log}({B}_{i,k}+{S}_{i,k})\right.\\ \quad -\,\left.\mathrm{log}({N}_{i,k}!)-({B}_{i,k}+{S}_{i,k})\right],\end{array}\end{eqnarray} \tag{ 4 }$$

where S_i,k is the number of expected events from our DM model, B_i,k is the number of observed background events in the region of interest (ROI) around the position of the galaxy, and N_i,k is the total number of observed counts. The sum is performed over i spatial bins and k f_hit fractional bins. The signal events S_i,k were obtained from the convolution of the photon flux of our DM model for five annihilation (decay) channels and the response matrix of the HAWC Observatory.

To measure how the DM model fits the data set D, we use the test statistic (TS) provided by the likelihood profile

$$\begin{eqnarray}&&\mathrm{TS}=-2\mathrm{log}\left(\displaystyle \frac{L({S}_{i,k}=0,{B}_{i,k})}{L({\hat{S}}_{i,k},{B}_{i,k})}\right).\end{eqnarray} \tag{ 5 }$$

The numerator is the maximum likelihood value for the null hypothesis, assuming S_i,k equal to zero in our DM model and the denominator is the maximum likelihood value when the DM model has a signal ${\hat{S}}_{i,k}$ different from zero. As usual, the statistical significance is obtained by $\sigma =\sqrt{\mathrm{TS}}$. We performed two kinds of analysis: for every dIrr Galaxy (individual source analysis) and combined analysis. We did not see any statistically significant excess from the data, and the significance is transformed into (one-sided) exclusion limits for the annihilation cross section and the decay lifetime of a generic DM candidate with masses in the range of [1, 100] TeV.

4.1. Individual Source Analysis

For the individual analysis, the values of the background and signal are estimated for spatial pixels within a circular ROI of 5° in diameter. The number of expected events for every source is computed by convolving the DM photon flux for the annihilation (decay) and the response matrix of HAWC. The DM photon flux is calculated assuming reference values for the annihilation cross section 〈σ_χ v〉^(ref) and the decay lifetime ${\tau }_{\chi }^{(\mathrm{ref})}$ of DM particles. Then, $-\mathrm{log}L({S}_{i,k},{B}_{i,k}| D)$ is computed and minimized using the Minuit package (Hatlo et al. 2005).

In all the cases the observations are consistent with the null hypothesis (no photons produced by annihilation or decay of DM particles), and the significances were converted into upper (lower) exclusion limits for the DM annihilation cross section (decay lifetime). To estimate the (one-sided) exclusion limit (at 95% C.L.), we found the value of the signal parameter ξ where the log-likelihood ratio changes by 2.71 with respect to the position of the maximum

$$\begin{eqnarray}&&\left[{\mathrm{TS}}^{(\max )}-\mathrm{TS}(\xi )\right]-2.71=0.\end{eqnarray} \tag{ 6 }$$

The signal parameter ξ is a global scale factor that only has an effect on the signal parameter. TS(ξ) is computed by

$$\begin{eqnarray}&&\mathrm{TS}(\xi )=2\displaystyle \sum _{i,k}\left[\mathrm{log}\left(\displaystyle \frac{{B}_{i,k}+\xi \times {S}_{i,k}}{{B}_{i,k}}\right)-\xi \times {S}_{i,k}\right].\end{eqnarray} \tag{ 7 }$$

When the signal parameter ξ is found, the value of the annihilation cross section and decay lifetime are obtained by the multiplication of the scale factor and the reference values used to compute the DM photon flux:

$$\begin{eqnarray}&&\langle {\sigma }_{\chi }v{\rangle }^{(95 \% )}=\langle {\sigma }_{\chi }v{\rangle }^{(\mathrm{ref})}\times \xi \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\tau }_{\chi }^{(95 \% )}={\tau }_{\chi }^{(\mathrm{ref})}\times {\xi }^{-1}.\end{eqnarray} \tag{ 9 }$$

The values obtained for 〈σ_χ v〉^(95%) and ${\tau }_{\chi }^{(95 \% )}$ do not depend on the value chosen for 〈σ_χ v〉^(ref) and ${\tau }_{\chi }^{(\mathrm{ref})}$, whose values are selected only for computational convenience. The methods described here are included in the analysis software for the HAWC Observatory as the Liff package (Younk et al. 2015).

4.2. Combined Analysis

We used the joint analysis technique to compute a combined limit for all the dIrr galaxies. The joint analysis allows us to estimate the value of parameters that are common to different sets of data; either observations of a source carried out by different experiments or observations of different sources performed by the same experiment. In our DM model, the common parameter for all the sources is the annihilation cross section (decay lifetime) used to compute the expected events for a specific source. The joint analysis is based on the joint likelihood function that results from the multiplication of the likelihoods for every data set. In the case of the joint analysis of a sample of astrophysical sources, the joint log-likelihood is

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}{L}_{\mathrm{joint}}(A\times {S}_{i,k,m},{B}_{i,k,m}| D)\\ \quad =\displaystyle \sum _{m=1}^{M}\mathrm{log}{L}_{d}(A\times {S}_{i,k,m},{B}_{i,k,m}| D),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{l}=\displaystyle \sum _{m=1}^{M}\displaystyle \sum _{i,k}\left[{N}_{i,k,m}\mathrm{log}({B}_{i,k,m}+A\times {S}_{i,k,m})\right.\\ \quad -\left.\,\mathrm{log}{N}_{i,k,m}!-({B}_{i,k,m}+A\times {S}_{i,k,m})\right],\end{array}\end{eqnarray} \tag{ 11 }$$

where S_i,k,m is the number of expected events for the mth source, the ith spatial bin, and the kth f_hit fractional bin. The common parameter for all the dIrrs is represented by a global normalization factor A. N_i,k,m and B_i,k,m are the observed counts and the background counts for the mth source, the ith spatial bin, and the kth nHit fractional bin, respectively. In this way, the joint likelihood is the usual method of taking into account the fluctuations of the background of every galaxy in the estimation of the signal parameter. Now, we can define the TS as usual:

$$\begin{eqnarray}&&{\mathrm{TS}}_{\mathrm{joint}}=-2\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{joint}}(S=0,{B}_{i,k,m})}{{L}_{\mathrm{joint}}(\hat{S},{B}_{i,k,m})}\right).\end{eqnarray} \tag{ 12 }$$

In the case of no detection, the statistical significances are converted into (one-sided) exclusion limits following the same method as in the individual source analysis. We use the Minuit package (Hatlo et al. 2005) to find the common signal that minimizes the negative of the joint likelihood function.

4.3. Expected Limits

Here, we show the exclusion limits obtained for 31 dIrr galaxies within the field of view of the HAWC Observatory. We present the limits for DM candidates with masses between 1 and 100 TeV. As explained in previous sections, the parameters of interest are the annihilation cross section, 〈σ_χ v〉, and the decay lifetime, τ_χ , of the DM particle. Here, we show and focus the discussion on the combined limits. The results for the individual limits are given in Appendices A and B.

5.1. Annihilation

The calculation of the combined upper limits on the annihilation cross section used for all galaxies in the sample. For simplicity, we assumed a point-source model for all galaxies. We show results for five annihilation channels: b, μ, τ, t, and W, see Figure 4. We observe that the most constraining limits are for the annihilation channel to τ lepton.

Zoom In Zoom Out Reset image size

We do not observe any difference between the combined limit (solid blue line) and the observed DDO 154 constraint (solid red line) in the lower mass range, see Figure 5. This effect is due to the large dispersion in the characteristics of the galaxies. In particular, galaxies with a small expected signal-to-background ratio contribute little to the combined limit.

Zoom In Zoom Out Reset image size

The combined limit clearly shows a separation from the galazy DDO 154 observed, though the magnitude of this effect depends on the channel. This separation starts at masses above 10 TeV for annihilation channels to μ and τ leptons; and above 30 TeV for the t, b, and W channels. Further studies are needed to fully address this effect, especially when including more galaxies in the combined analysis. However, the improvement for massive DM candidates shows the importance of combined analysis to obtain stronger constraining limits on DM parameters.

In addition, we calculate the expected limit for the DDO 154 galaxy. We obtain the average value, and 68% and 95% CIs after fluctuating the background events in the HAWC maps, the solid-black line in Figure 5. We chose the DDO 154 galaxy because it has the most constraining limit of all the galaxies in our sample and contributes most strongly to the combined limit.

We observe that, in general, the observed limit for the galaxy DDO 154 (red line) is more restrictive than the expected DDO 154 limit (black line) for DM candidates with mass below 10 TeV (30 TeV) for μ and τ (b, t, and W) channels. For more massive DM particles, the observed DDO 154 limit is closer to the expected DDO 154 exclusion limit. This behavior is related to the fact that the observed and expected counts of the background are approximately equal to the high f_hit.

For lower masses, the discrepancy between the observed DDO 154 and the expected DDO 154 limits may be explained by a 2σ under-fluctuation observed in the TS map for the region around the galaxy DDO 154. This under-fluctuation tells us that a DM model is likely less probable to describe the data observed in the galaxy DDO 154, and more consistent with the background-only hypothesis.

Finally, Figure 5 also shows the limits obtained from observations of other experiments. The data sets correspond to the exclusion limits obtained from observations of the galaxy WLM with H.E.S.S. (orange line), Abdallah et al. (2021); and the combined limit for seven dIrr galaxies within the field of view of the Fermi-LAT experiment (yellow line), Gammaldi et al. (2021). We can observe that there is a nice complementarity between the three exclusion limits. Also, we observe that the HAWC combined limit is the most restrictive for dIrr galaxies for masses above 10 TeV.

5.2. Decay

While it is often supposed that WIMPs are stable particles, in some supersymmetric models, the DM candidate can decay to SM particles. Adding a term that breaks down R parity allows this process. To match the DM relic density at the present epoch, the candidate must live much longer than the age of the universe (Ando & Ishiwata 2015). In supersymmetric models, the estimated lifetime is usually τ_χ ∼ 10²⁷ s (Buchmuller et al. 2007; Choi et al. 2014). Figure 6 shows the lower combined limits obtained for our galaxy sample, and Figure 7 also shows the expected limits for the galaxy DDO 154. As in the annihilation case, we observe that the most restrictive combined limit is for the decay channel to τ leptons, where we obtain a lifetime >10²⁶ s for masses above ∼20 TeV. We note the same behavior for the combined and DDO 154 limits, with an improvement in the high-mass regime showing the importance of combined analysis for these targets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.3. Comparison with dSph Galaxies

We have obtained the exclusion limits under the assumption that dIrr galaxies do not show evidence of astrophysical processes contributing to the gamma-ray flux. The same hypothesis is valid for dSph galaxies, and we can compare the results between both target classes. For the comparison, we used the dSph limits reported in Albert et al. (2018). There exists a classification of the dSph galaxy population based on the number of star members in the galaxy: classical and ultrafaint dSph galaxies. The designation also refers to the quality of available kinematical data, which impacts the level of uncertainty on the astrophysical factors: better data, smaller uncertainty. We used the DDO 154 and combined limits from the dIrr sample to observe possible similarities between the populations. The objective is to encourage further studies combining these galaxy populations and put more constraining limits on the DM parameter space.

Figures 8 and 9 show a comparison between dIrr and dSph galaxies for DM particles annihilating or decaying to five different channels. The black and red lines represent the dIrr combined and DDO 154 limits, respectively. The dotted lines show the exclusion limits for classical (Draco, Leo I, Leo II, and Ursa Minor) dSphs and the solid lines for ultrafaint (Triangulum II, Segue I, and Coma Berenices) dSph galaxies. We observe that, in general, the limits from dIrr galaxies are less restrictive than the ultrafaint dSphs up to 2 orders of magnitude for annihilation and 1 order of magnitude for decay. The reason is that ultrafaint galaxies have astrophysical factors larger than our sample. Also, one must remember that the associated uncertainty is considerably higher for ultrafaint galaxies. For example, in recent studies (Pace & Strigari 2018), if the astrophysical factor obtained for Triangulum II is an upper limit, then the actual constraint could be less restrictive. We also observe that the combined limit for dIrr galaxies is comparable (has the same order of magnitude) to limits from classical dSph galaxies. The reason is that the astrophysical factors between the two populations are very similar. Note that this does not imply that both galaxy populations have or share the same properties (dIrr galaxies are farther than spheroidals, but dIrr are more massive, for example).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also observe that the combined limits of the dIrr and dSph galaxies have similar behavior. The improvement is relatively small in both populations, probably caused by background fluctuations in galaxies in decl. bands where HAWC has lower sensitivity, or by an effect due to the extension of the target. They both impact the calculation and minimization of the likelihood function. As explained in Section 4.2, the first scenario can lead to a decrease in the best limit when different galaxies in the sample have background counts differing from each other. In our case, one possible solution could be only to use galaxies within decl. bands around the target that contributes the most to the combined limit (DDO 154). Also, including only targets with similar expected signal intensity (or J and D factors) may help increase the combined limit. We will explore this possibility in a future analysis.

We have shown that dIrr galaxies are a new suitable population of objects to perform indirect DM searches by extended-array experiments, like the HAWC Observatory. This is possible because of the environmental conditions of these galaxies, with low SFRs and a low population of massive stars. We computed the exclusion limits at 95% C.L. for the annihilation and decay of generic WIMPs with masses between 1 and 100 TeV. In both cases, the best limits are obtained for the galaxy DDO 154, located in a decl. band where HAWC has good sensitivity. The computed limits were obtained under the assumption that the DM density is described by a Burkert profile and no enhancement by substructure was considered, but it is possible that the presence of a central black hole in these galaxies could modify the slope of the DM distribution in the central region of dIrrs. This would lead to an enhancement in the annihilation exclusion limit, at least for 1 or 2 orders of magnitude, according to the dynamic history of the host galaxy (Gonzalez-Morales et al. 2014), and the fact that the profile is modified via the rotational motion of the galaxy (Oman et al. 2018), leading to more cusped profiles, like the Navarro–Frenk–White profile (Navarro et al. 1997). More studies on this approach are planned for an upcoming publication. Also, we observe more dIrr galaxies in the local universe (∼11 Mpc), and we can obtain their mass distribution from rotation curve data. Therefore, combined analysis as described in this work will lead to more stringent constraints on DM candidates. Moreover, given the established universality of their mass distribution, useful constraints on the DM cross sections can be obtained even for a limited amount of observations for each of them. We also show that our exclusion limits are comparable to the limits obtained for classical dSph galaxies, so future studies may consider a possible joint analysis between the two populations to be able to constrain even more the properties of DM candidates.

The goal of this paper is to show the capabilities of the HAWC Observatory to use large populations of galaxies to constrain parameters to DM models through the search of gamma-ray signatures of annihilation or decay of DM particles. Here, we used a sample of dIrr galaxies. However, the present analysis is the starting point for more detailed studies. In this analysis, we assumed that the DM-induced gamma-ray emission is described by a point source, but this is not the case for galaxies like IC 10, IC 1613, and NGC 6822, with an extension across the sky of ≳2° (see Table 2). Using an extended emission model should lead to more realistic constraints. Furthermore, with an analysis comprising a longer period of time and improved HAWC energy estimators, we can improve the constraints, both on annihilation and decay. We may also obtain limits to DM candidates' masses in the high-end range of masses expected for WIMPs (Enqvist & Kainulainen 1991; Smirnov & Beacom 2019). Therefore, we will be probing the hypothesis of the thermal production of WIMPs in the early universe. Note that this can be an indication of other interesting scenarios as multicomponent or light DM.

Second, as we pointed out at the beginning of this section, dIrr galaxies are abundant in the local universe and combined analysis using more targets should be addressed. Moreover, as we discussed in Section 5.3, this also gives us the opportunity to test different strategies for combining the data from dIrr galaxies. For example, in order to reduce the background fluctuations we may use targets with a similar expected signal-to-background ratio in the same decl. band to check if the combined limit can be improved.

Finally, combined limits using not only the dIrr galaxy populations, but also classical dSph galaxies should help increase the constraints on the different DM parameters.

We acknowledge the support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grant Nos. 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grant Nos. IG101320, IN111716-3, IN111419, IA102019, IN110621, IN110521, IN102223; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant No. DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society-Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, NXPO (grant No. B16F630069); Coordinación General Académica e Innovación (CGAI-UdeG), PRODEP-SEP UDG-CA-499; Institute of Cosmic Ray Research (ICRR), University of Tokyo, H.F. acknowledges support by NASA under award No. 80GSFC21M0002. We also acknowledge the significant contributions over many years of Stefan Westerhoff, Gaurang Yodh, and Arnulfo Zepeda Dominguez, all deceased members of the HAWC Collaboration. We thank Scott Delay, Luciano Díaz, and Eduardo Murrieta for their technical support. Viviana Gammaldi's contribution to this work has been supported by Juan de la Cierva-Incorporación IJC2019-040315-I grants, by grant Nos. PGC2018-095161-B-I00, CEX2020-001007-S, and PID2021-125331NB-I00 all funded by MCIN/AEI/10.13039/501100011033 and by ERDF "A way of making Europe."

Here, we present the individual exclusion limits for all the targets in our sample. For completeness, we also show the combined limit as described in Section 4.2.

In this section, we show the exclusion limits for the 31 dIrr galaxies in our sample, see Figure 10. We observe that the best individual constraining limit is for the galaxy DDO 154. While its astrophysical factor is not the largest in our sample, the galaxy DDO 154 is located in a decl. band where HAWC has good sensitivity. For completeness, we also show the combined limit in Figure 10. We can also observe that the DDO 154 and combined limits are similar for energies below 10 TeV. For channels, such as μ and τ, the combined limit shows an improvement above 10 TeV. For the other annihilation channels, the improvement in the combined limit only appears at the high-energy regime. An explanation could be the high background rejection at high energies by HAWC, providing cleaner sample data for all the dIrr galaxies.

Zoom In Zoom Out Reset image size

As we pointed out in Section 5, another explanation for the null difference between the combined and the DDO 154 limits is the large scatter in the properties of all galaxies in the sample. In particular, the individual backgrounds should contribute to this effect.

In this section, we show the exclusion limits for the 31 dIrr galaxies in our sample (see Figure 11). Again, we observe that the best individual limit is for the galaxy DDO 154. As in the annihilation case, the probable reason is that the galaxy DDO 154 is in a favorable decl. band. For completeness, we also show the combined limit. We observe the same features between the combined and the DDO 154 limits. As we discuss in Section 5, the decay limits (individual and combined) for dIrr galaxies are competitive with those obtained for other targets such as dSph galaxies. The main reason is the high mass of dIrr galaxies.

Zoom In Zoom Out Reset image size