Search for Ultraheavy Dark Matter from Observations of Dwarf Spheroidal Galaxies with VERITAS

Ultraheavy dark matter (UHDM) presents an alternative mass range for dark matter, and is partly motivated by the absence of a dark matter signature in the well-explored mass ranges suggested by the simplest dark matter models. Most dark matter searches have focused on the mass range of sub-eV (e.g., axion-like particles) or ∼GeV–TeV scales (e.g., weakly interacting massive particles). However, the dark matter particle mass is simply not limited to such ranges; indeed, there are many viable UHDM candidates (for a brief introduction, see Carney et al. 2022). If dark matter emerged as a thermal relic from the early universe, then as unitarity places an upper bound on its annihilation into Standard Model (SM) particles, this naively prohibits masses above (100) TeV if dark matter is a point-like particle (Griest & Kamionkowski 1990). Roughly, unitarity sets a maximal value for the dark matter annihilation cross section, and for dark matter heavier than ∼100 TeV, even the largest allowed cross section is insufficient to reduce the equilibrium abundance of dark matter to the observed value. However, if dark matter is made of composite states with geometrical cross sections (i.e., the cross section scales as π R², where R is the intrinsic size of the dark matter particle), the limit is easily evaded (e.g., Harigaya et al. 2016; Geller et al. 2018). One can also consider scenarios where the dark matter is not a simple thermal relic, with or without compositeness (e.g., Berlin et al. 2016; Contino et al. 2019). As discussed in Tak et al. (2022), the annihilation of UHDM particles can produce a γ-ray signal in the form of monoenergetic γ-ray lines in addition to a continuum contribution of photons with energy equal to and below the dark matter particle mass (E_γ ≲ M_χ ), with the exact spectrum determined by the particle physics underlying the annihilation. Given this, the authors demonstrated that current very-high-energy (VHE; ≥ 100 GeV) γ-ray observatories are sensitive to an annihilation signal from UHDM, for masses up to at least a few tens of PeV.

Among the best targets for indirect dark matter searches are dwarf spheroidal galaxies (dSphs) of the Local Group (located ∼20–200 kpc from Earth). Since they are dark matter–rich regions without known nearby VHE sources, ²⁷  they have been widely studied with current VHE observatories (e.g., Aleksić et al. 2014; Albert et al. 2018, 2020; Abdalla et al. 2018; Acciari et al. 2022). For instance, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) observed five dSphs and provided upper limits on the dark matter velocity-weighted annihilation cross section in the mass range from 100 GeV to 100 TeV (Archambault et al. 2017).

In this work, we revisit the VERITAS observations of four of the five dSphs (Segue 1, Ursa Minor, Boötes, and Draco; in total 216 hr of observations) and search for the indirect UHDM signal up to a mass of 30 PeV. The observation times for the targets are listed in Table 1. Note that we consider dSphs for which we have an estimate of the dark matter density profile; the Willman observation included in Archambault et al. (2017) is excluded. We derive upper limits on the UHDM velocity-weighted annihilation cross section from a joint-fit maximum likelihood estimation (MLE) analysis. We further interpret the derived limits in terms of the allowed radius of a composite UHDM particle.

Table 1. Table of the Four Dwarf Spheroidal Galaxies Considered by VERITAS in This Analysis, Showing the VERITAS Observational Results in the First Five Columns and the Assumed Properties of the Dwarf Galaxies

	Non	Noff	α	tobs	σ	ρs	rs	α	β	γ
				(hr)		(M⊙/pc3)	(pc)				(deg)	(GeV2 cm−5 · sr)
Segue 1	15895	120826	0.131	92.0	0.7	1.78	3.1 × 102	0.54	4.36	0.64	0.35	2.5 × 1019
Ursa Minor	4181	35790	0.119	60.4	−0.1	5.6 × 10−1	3.6 × 102	2.37	8.77	1.2 × 10−2	1.19	7.1 × 1018
Boötes	1206	10836	0.116	14.0	−1.0	6.7 × 10−4	1.2 × 104	2.81	4.87	1.08	0.47	1.7 × 1018
Draco	4297	39472	0.111	49.8	−1.0	8.2 × 10−3	2.6 × 103	1.96	6.09	0.95	1.41	1.3 × 1019

Note. Columns (1)–(5) show the counts recorded by VERITAS in the ON and OFF regions, ratio between the areas of the ON and OFF regions, and the exposure times and detection significances, respectively. The next five columns give the selected parameter set of the generalized NFW profile (Geringer-Sameth et al. 2015, described in Section 3) for the four dwarf spheroidal galaxies considered. The final two columns show the maximum angular distance considered in the J-factor calculation and the J-factor, respectively.

Download table as:  ASCIITypeset image

VERITAS is an array of four Imaging Atmospheric Cherenkov Telescopes (IACTs). The instrument is located at the Fred Lawrence Whipple Observatory in southern Arizona (31° 40′ N 110°57′ W). The telescope optics utilize a Davies–Cotton design. The reflectors are 12 m in diameter and composed of 350 hexagonal mirrors. The VERITAS cameras are composed of 499 photomultiplier tubes (PMTs), and have a field of view of 35 (Holder et al. 2008). VERITAS precisely reconstructs γ-rays with energies between ∼100 GeV and ∼30 TeV and is sensitive to even higher energy γ-rays, up to ∼100 TeV. This is of particular relevance for this study, corresponding to sensitivity to an annihilation signal from a UHDM particle with mass up to a few tens of PeV. The angular resolution of VERITAS is ∼01 at 1 TeV (68% containment), while the energy resolution is 15%–20% at 1 TeV. VERITAS can detect a point source with a flux of 1% of the Crab Nebula flux in ∼25 hr of observation (Park et al. 2015).

Observations of the four dSphs considered here were made between 2007 and 2013. During this time period, VERITAS underwent two upgrades. The first took place in the summer of 2009, in which the position of one of the telescopes was altered to produce a more symmetric array. The second upgrade was made in the summer of 2012, in which the camera PMTs were exchanged for a model with a higher quantum efficiency and the trigger system was upgraded, yielding a 50% increase to the photon collection efficiency (Kieda et al. 2013). As the sensitivity of the instrument and the value of the energy threshold changed with each of these upgrades, dedicated Monte Carlo models and instrument response functions (IRFs; including effective areas, energy dispersion matrices, and point-spread functions) are available for each of the three array epochs. All data were collected in wobble mode (Fomin et al. 1994).

Data were reduced using one of the standard VERITAS calibration and event reconstruction pipelines (Cogan 2007). As described in Archambault et al. (2017), a novel crescent-background technique was used to define the OFF region for background estimation, while the ON region was centered on the target location. The number of ON and OFF counts, the ratio α between the size of the ON and OFF regions, and the detection significance (Li & Ma significance; Li & Ma 1983) is given in Table 1. No low-level data reanalysis was performed; the event lists and IRFs from Archambault et al. (2017) were used for this analysis.

In the previous VERITAS dark matter study using dSphs (Archambault et al. 2017), the so-called event-weighting method (Geringer-Sameth et al. 2015) was exploited to search for a dark matter signature in the observed data. In this work, we rather adopt a commonly used and extensively documented method, MLE, and perform an unbinned likelihood analysis. To perform the MLE analysis, we introduce a likelihood function, quantifying the consistency of the observed dSph data (D) with a given dark matter model,

This likelihood is a product of the likelihoods modeling the total counts in the ON and OFF regions, as well as the predicted energy distribution of the counts in the ON region. In more detail, N_on and N_off are the number of observed ON- and OFF-region counts, respectively, and α is the relative exposure time between the ON and OFF regions. The nuisance parameter b is the expected number of background counts. Two probability density functions (PDFs) are required in this unbinned likelihood function: one for the dark matter signal (p_s ) and the other for the background (p_b ). The background PDF is obtained from the normalized OFF-region event distribution. The dark matter signal PDF and the dark matter signal counts (N_s ) expected to be observed by the instrument within the ON region, of size ΔΩ, are determined by the dark matter spectrum ( and J-factor (J(ΔΩ)), which is the square of the dark matter density integrated along the line of sight within the ON region. In detail,

Here 〈σ v〉 and M_χ are the velocity-averaged dark matter annihilation cross section and dark matter particle mass, respectively. Although not shown here, these results are convolved with the IRF of VERITAS to obtain p_s and N_s , which accounts for the finite angular and energy resolution of the instrument. For more details, see Archambault et al. (2017).

For the γ-ray spectrum from dark matter annihilation at production, , we use HDMSpectra (Bauer et al. 2021) ²⁸  instead of the widely used PPPC4DMID spectrum (Cirelli et al. 2011). This is because the former provides dark matter annihilation spectra for various channels in a broad mass range from 1 TeV up to the Planck energy, while PPPC4DMID extends only to a dark matter mass of 100 TeV. With HDMSpectra, we obtain a set of nine final-state photon spectra, assuming a 100% branching ratio of dark matter particles in nine different annihilation channels: e⁺ e⁻, μ⁺ μ⁻, τ⁺ τ⁻, , , W⁺ W⁻, ZZ, γ γ, and . In considering the differences between the production spectrum and the photon spectrum observable by VERITAS, it is important to note that the UHDM signature (from e.g., the annihilation of a 30 PeV dark matter particle) results in observed γ-rays below 100 TeV. Consequently, absorption on ambient photon fields can be ignored.

For the dark matter density profile, ρ(r), we adopt the generalized Navarro–Frenk–White (NFW) profile (Hernquist 1990; Zhao 1996; Geringer-Sameth et al. 2015), which is a function of five parameters,

The values of the free parameters used for each dSph and the resulting (unconvolved) J-factors are given in Table 1. The parameters are adopted from Geringer-Sameth et al. (2015).

For the joint-fit analysis, in which data from the four dSphs are combined to maximize statistical power, we combine the individual likelihood functions to form a joint one,

The significance of the dark matter signal over background can be obtained by comparing two likelihoods,

If the significance of the dark matter signal is below the threshold to claim a detection (λ ≳ 25), we compute an upper limit on the dark matter velocity-weighted annihilation cross section by using the likelihood profile. The one-sided 95% confidence level upper limit on the dark matter velocity-weighted annihilation cross section is the value of the cross section corresponding to of 1.35 compared to the likelihood maximum.

We do not detect a UHDM signal above background. From the individual and the joint-fit analyses, we obtain λ less than our threshold in all annihilation channels and for all masses from 1 TeV to about 30 PeV (see Appendix A).

4.1. Upper Limits on the UHDM Velocity-weighted Annihilation Cross Section

We compute upper limits on the dark matter velocity-weighted annihilation cross section for each channel. Figure 1 shows upper limits obtained in the joint-fit analysis, each with a systematic uncertainty band resulting from the limited understanding of the dark matter density distribution. The uncertainty band is obtained from 300 realizations with different dark matter density profile parameter sets from Geringer-Sameth et al. (2015); each parameter set can sufficiently describe the stellar-kinematic data observed from the selected target. In the case of Segue 1, the ambiguity of selecting member stars significantly affects the dark matter profile, such that the total density can differ by 2 orders of magnitude (Bonnivard et al. 2016). For this reason, we additionally present the combined upper limits excluding the Segue 1 data.

Zoom In Zoom Out Reset image size

We note that the discontinuity in the γγ channel at around 100 TeV, the maximum value for which we consider γ-ray events, is expected. Above 100 TeV, the dominant contribution from the delta function/line annihilation signal at M_χ = E_γ results in final-state γ-rays whose energies are above the VERITAS sensitive energy range, leaving only the continuum spectrum. The continuum spectrum is more challenging to detect in comparison to a line signature, resulting in less-sensitive limits when the line component is no longer detectable.

4.2. Comparison with the Background Fluctuation

We test whether the distribution of ON-region events is consistent with the Poisson fluctuation of the background. To do this, we estimate an expected upper limit from a simulated ON region for which events and their energy are randomly selected from the observed OFF-region events. The number of simulated ON events is selected from a Poisson distribution with mean equal to the observed number of OFF-region events, scaled by the ratio of the areas of the ON and OFF regions, . For each channel, we repeat this process 300 times and obtain an expected upper-limit band with the width determined by the magnitude of the background fluctuation.

Figure 2 shows the comparison of the observed upper limits with the expected upper-limit bands. Each solid line (blue) is an upper-limit curve from the parameter set listed in Table 1, and the expected upper-limit band is depicted in orange with 68% (solid) and 95% (dotted line) containment. For all annihilation channels, the observed upper limits are consistent with the expected upper limits within the 95% confidence level. This result supports the nondetection of the UHDM annihilation signal, as well as quantifying the impact of statistical uncertainty on the derived limits.

Zoom In Zoom Out Reset image size

As mentioned at the outset, a requirement for UHDM is to evade the so-called unitarity limit. The (100 TeV) bound assumes that the dark matter is a point-like particle, which is in thermal equilibrium with SM particles in the early universe. However, one straightforward way to evade this limit takes point-like dark matter that captures into bound states. These additional channels can achieve a larger annihilation cross section while respecting unitarity. Individual partial-wave contributions must respect their associated unitarity bound, but the total cross section is given by the sum of all partial-wave contributions. This effect is even seen in medium-sized representations of electroweak SU(2), allowing them to be simple thermal-relic UHDM candidates (Bottaro et al. 2022).

One class of UHDM models that further relaxes the unitarity bounds on mass are composite dark matter models, where UHDM is not a point-like particle and thus possibly has a geometrical cross section. In the case of an interaction with a geometrical cross section, the unitarity bound becomes

where v_rel is the average velocity between dark matter particles (in our case, in dSph halos), and R is the size of the particle. Note that the unitarity limit for a point-like particle can be reproduced with R = 0, whereas if the particle mass is large enough, we can reproduce the classical cross section of 〈σ v〉 = 4π R² v_rel. For a detailed discussion, see Tak et al. (2022). Figure 3 shows our upper limits for two annihilation channels (blue solid lines), τ⁺ τ⁻ and , as well as the theoretical bounds: the standard thermal-relic limit (red solid), the unitarity limit for a point-like particle respecting the partial-wave unitarity bound (purple solid), and the unitarity limits for a composite particle (purple dashed lines). Note that we assume v_rel/c = 2 × 10⁻⁵ for the relative velocity between dark matter particles in dSph galaxies (Martinez et al. 2011; McGaugh et al. 2021). Our results not only constrain part of the allowed region of a point-like dark matter cross section, but also limit the radius of UHDM in a mass range from about 100 TeV to 30 PeV. This is visible from Figure 3 and depicted explicitly in Figure 4. For example, below a dark matter mass of approximately 1 PeV, a UHDM model with the UHDM particle size of 0.6 fm or larger can be rejected at the 95% confidence level in all annihilation channels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows our upper limits compared with results from the Fermi-LAT, MAGIC, VERITAS, H.E.S.S., and HAWC collaborations. Since we use the previously published VERITAS observations, our results are similar to the published ones below 100 TeV, with the differences coming entirely from the method of extracting upper limits. Our results extend limits on the dark matter velocity-weighted annihilation cross section into a mass range which has not previously been explored.

Zoom In Zoom Out Reset image size

In this paper, we have presented an indirect search for a UHDM annihilation signal, using previously published VERITAS observations to access a novel dark matter parameter space. We search for final-state γ-rays from nine annihilation channels, using 216 hr of observations of four dwarf spheroidal galaxies: Segue 1, Ursa Minor, Boötes, and Draco. In the absence of a detection, we have shown upper limits on the dark matter velocity-weighted annihilation cross section for dark matter particle masses from 1 TeV to 30 PeV with a joint-fit MLE analysis. This work has reported a new UHDM search with IACT observations, detailed a robust method for such searches, and should provide insight for future UHDM studies with the deep observations from the current IACTs and/or the future sensitive observatories such as the Cherenkov Telescope Array (Acharya et al. 2019).

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation, and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science, and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument. D.T. and E.P. acknowledge the Young Investigators Program of the Helmholtz Association, and thank Alex Geringer-Sameth, Savvas M. Koushiappas, and Matthew Walker for providing the parameter sets for the J-factors. D.T. is supported by the National Research Foundation of Korea (NRF) grant, No. 2021M3F7A1084525, funded by the Korea government (MSIT). M.B. is supported by the DOE (HEP) Award DE-SC0019470.

Figure 6 shows the signal significance (given by ) as a function of dark matter particle mass in the nine annihilation channels. The significance curves for the individual dSphs are shown, as well as the combined results. For no dark matter particle mass, dSph, or annihilation channel does the signal significance reach 2σ.

Zoom In Zoom Out Reset image size

We note that the significance in Figure 6 is calculated from the likelihood analysis with observed ON and OFF regions. This result shows the nondetection of a dark matter signal. Figure 2, in contrast, compares observed upper limits with expected upper-limit bands assuming a simulated ON region made up of randomly sampled observed OFF-region events. The observed agreement between the observed limits and expected limit band implies that observed ON region is consistent with the Poisson fluctuation of observed OFF regions. These two approaches lead to the same conclusion that we do not observe any excess (a possible dark matter signal) from our observations.

Figure 7 shows the upper limits on the UHDM velocity-weighted annihilation cross section as a function of particle mass for each dSph considered, as well as the combined limit. As in the main text, nine annihilation channels are considered. As expected based on the J-factors listed in Table 1, Segue 1 generally provides the most constraining limits, followed by Ursa Minor and Draco, with the weakest limits coming from Boötes.

Zoom In Zoom Out Reset image size