VERITAS Collaboration

An indirect search for dark matter with a combined analysis of dwarf spheroidal galaxies from VERITAS

A. Acharyya CP3-Origins, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark C. B. Adams Physics Department, Columbia University, New York, NY 10027, USA P. Bangale Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA J. T. Bartkoske Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA P. Batista DESY, Platanenallee 6, 15738 Zeuthen, Germany W. Benbow Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA J. L. Christiansen Physics Department, California Polytechnic State University, San Luis Obispo, CA 94307, USA A. J. Chromey Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA A. Duerr Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA M. Errando Department of Physics, Washington University, St. Louis, MO 63130, USA A. Falcone Department of Astronomy and Astrophysics, 525 Davey Lab, Pennsylvania State University, University Park, PA 16802, USA Q. Feng Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA G. M. Foote Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA L. Fortson School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA A. Furniss Department of Physics, California State University - East Bay, Hayward, CA 94542, USA W. Hanlon Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA D. Hanna Physics Department, McGill University, Montreal, QC H3A 2T8, Canada O. Hervet Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA C. E. Hinrichs Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA and Department of Physics and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH 03755 USA J. Holder Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA T. B. Humensky Department of Physics, University of Maryland, College Park, MD, USA and NASA GSFC, Greenbelt, MD 20771, USA W. Jin Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA M. N. Johnson Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA P. Kaaret Department of Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA M. Kertzman Department of Physics and Astronomy, DePauw University, Greencastle, IN 46135-0037, USA D. Kieda Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA T. K. Kleiner DESY, Platanenallee 6, 15738 Zeuthen, Germany N. Korzoun Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA S. Kumar Department of Physics, University of Maryland, College Park, MD, USA M. J. Lang School of Natural Sciences, University of Galway, University Road, Galway, H91 TK33, Ireland M. Lundy Physics Department, McGill University, Montreal, QC H3A 2T8, Canada G. Maier DESY, Platanenallee 6, 15738 Zeuthen, Germany C. E McGrath Conor McGrath, mcgrathconor94@gmail.com School of Physics, University College Dublin, Belfield, Dublin 4, Ireland M. J. Millard Department of Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA C. L. Mooney Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA P. Moriarty School of Natural Sciences, University of Galway, University Road, Galway, H91 TK33, Ireland R. Mukherjee Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA W. Ning Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA S. O’Brien Physics Department, McGill University, Montreal, QC H3A 2T8, Canada and Arthur B. McDonald Canadian Astroparticle Physics Research Institute, 64 Bader Lane, Queen’s University, Kingston, ON Canada, K7L 3N6 R. A. Ong Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA N. Park Department of Physics, Engineering Physics and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada M. Pohl Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam-Golm, Germany and DESY, Platanenallee 6, 15738 Zeuthen, Germany E. Pueschel Fakultät für Physik & Astronomie, Ruhr-Universität Bochum, D-44780 Bochum, Germany J. Quinn School of Physics, University College Dublin, Belfield, Dublin 4, Ireland P. L. Rabinowitz Department of Physics, Washington University, St. Louis, MO 63130, USA K. Ragan Physics Department, McGill University, Montreal, QC H3A 2T8, Canada P. T. Reynolds Department of Physical Sciences, Munster Technological University, Bishopstown, Cork, T12 P928, Ireland D. Ribeiro School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA E. Roache Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA J. L. Ryan Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA I. Sadeh DESY, Platanenallee 6, 15738 Zeuthen, Germany L. Saha Center for Astrophysics

|

Harvard & Smithsonian, Cambridge, MA 02138, USA G. H. Sembroski Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA R. Shang Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA M. Splettstoesser Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA D. Tak Donggeun Tak, donggeun.tak@gmail.com SNU Astronomy Research Center, Seoul National University, Seoul 08826, Republic of Korea. A. K. Talluri School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA J. V. Tucci Department of Physics, Indiana University-Purdue University Indianapolis, Indianapolis, IN 46202, USA V. V. Vassiliev Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA A. Weinstein Department of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA D. A. Williams Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA S. L. Wong Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

(July 25, 2025)

Abstract

Understanding the nature and identity of dark matter is a key goal in the physics community. In the case that TeV-scale dark matter particles decay or annihilate into standard model particles, very-high-energy (VHE) gamma rays (greater than 100 GeV) will be present in the final state. The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an imaging atmospheric Cherenkov telescope array that can indirectly detect VHE gamma rays in an energy range of 100 GeV to $>$ 30 TeV. Dwarf spheroidal galaxies (dSphs) are ideal candidates in the search for dark matter due to their high dark matter content, high mass-to-light ratios, and their low gamma-ray fluxes from astrophysical processes. This study uses a legacy data set of 638 hours collected on 17 dSphs, built over 11 years with an observing strategy optimized according to the dark matter content of the targets. The study addresses a broad dark matter particle mass range, extending from 200 GeV to 30 PeV. In the absence of a detection, we set the upper limits on the dark matter velocity-weighted annihilation cross section.

Dark Matter

I Introduction

Cold dark matter is a key component of the current cosmological picture, comprising $\sim$ 85% of the matter of the universe [1]. Indirect astrophysical dark matter searches provide a means of probing the nature of dark matter, complementary to direct detection and collider searches. Assuming that dark matter is made up of particles, and that these particles interact with standard baryonic matter, it is possible to search for annihilation or decay of dark matter particles to Standard Model particles. The focus of this paper is searching for dark matter self-annihilation signatures, assuming self-conjugate dark matter.

Gamma rays are typically present in the final state following dark matter annihilation, whether as mono-energetic “line” signatures from promptly produced photons, or as a continuum gamma-ray spectrum with a cut-off at the dark matter particle mass. For the weakly interacting massive particle (WIMP) class of dark matter candidates, the preferred dark matter mass lies in the GeV to TeV range. The annihilation of such particles leads to gamma rays in the sensitive energy range of ground- and space-based gamma-ray instruments. Indirect searches for dark matter annihilation with gamma rays have been conducted and limits set on the dark matter velocity-weighted annihilation cross section with the space telescope $Fermi$ -LAT [2, 3, 4], the imaging atmospheric Cherenkov telescopes (IACTs) H.E.S.S. [5, 6, 7, 8], MAGIC [9, 10] and VERITAS [11], and the water Cherenkov observatory HAWC [12, 13, 14]; the cited publications highlight the most recent results.

The Galactic Center and the Milky Way satellite dwarf spheroidal galaxies (dSphs) are canonical targets for indirect searches with gamma rays from dark matter annihilation. While the Galactic Center is the nearest large repository of dark matter and hence has the largest predicted annihilation signal, dSphs offer a set of targets with low astrophysical backgrounds and modest angular extensions compared to the point-spread functions of the gamma-ray instruments. The dSphs can be divided into classical and ultra-faint objects. The former contain on the order of hundreds of stars, while the latter contain on the order of tens of stars. As stellar motions are used to constrain the dark matter content of these systems, estimates of the dark matter content of classical dSphs tend to be more robust, as their stellar populations are well-measured. This paper analyzes the complete set of observations of dSphs, including classical and ultra-faint targets, collected by VERITAS over an 11 year period, and it utilizes the lack of a detected signal to derive upper limits on the dark matter velocity-weighted annihilation cross section in the standard WIMP mass range below $\sim$ 100 TeV.

A recent work by [15] has highlighted the theoretical motivation for dark matter annihilation searches outside of the standard WIMP mass range. The unitarity limit at $\sim$ 194 TeV [16, 17] is violated by point-like thermal-relic dark matter. However, dark matter composed of composite states or non-thermal relics can evade this limit [e.g., 18, 19]. Also, [15] showed that the annihilation of ultra-heavy dark matter particles (UHDM) up to 30 PeV can produce a sufficient number of photons below 100 TeV to enable VERITAS to constrain the velocity-averaged cross section for PeV dark matter. This UHDM mass range has been probed with gamma-ray and neutrino searches [20, 21], but a broad range of dark matter annihilation cross sections remains unconstrained. Consequently, we also present limits on annihilation of UHDM, with particle masses between 194 TeV and 30 PeV.

II Predicted gamma-ray flux and target selection

The predicted observable gamma-ray flux from a dSph, $\phi_{s}$ , is given by the product of two terms; the “particle physics factor” ( $\Phi_{PP}$ ) which is based on theoretical models of the annihilation of dark matter particles producing gamma-ray products, and the “astrophysical factor” or “J-factor”, which is determined by the dark matter content of the region observed:

\phi_{s}(\Delta\Omega)=\underbrace{\frac{1}{4\pi}\frac{\left<\sigma v\right>}{2M^{2}_{\chi}}\int^{E_{max}}_{E_{min}}\frac{dN}{dE}dE}_{\Phi_{PP}}\times\underbrace{\int_{\Delta\Omega}\frac{dJ}{d\Omega}d\Omega}_{\text{J-factor}},

(1)

where $M_{\chi}$ is the dark matter particle mass, $\left<\sigma v\right>$ is the thermally averaged, velocity-weighted annihilation cross-section of the dark matter particles, $dN/dE$ is the final-state photon spectrum, and $dJ/d\Omega$ is the differential annihilation J-factor with respect to solid angle $\Omega$ (Equation 2).

The photon spectrum from dark matter annihilation depends on the annihilation channel of the dark matter pair into Standard Model particles. Here we assume a 100% branching ratio of dark matter particles into the following channels: $e^{+}e^{-}$ , $\mu^{+}\mu^{-}$ , $\tau^{+}\tau^{-}$ , $t\bar{t}$ , $b\bar{b}$ , $W^{+}W^{-}$ , $ZZ$ , $\gamma\gamma$ , and $\nu_{e}\bar{\nu}_{e}$ . For each of the annihilation channels, we exploit the spectrum provided by PPPC4DMID, which offers trustworthy spectra in the 5 GeV to 100 TeV energy range [22], and HDMSpectra, which provides reliable spectra from 1 TeV to the Planck energy [23]. The latter includes many of the interactions in the full unbroken Standard Model that are relevant for estimating accurate UHDM annihilation spectra. Thus, we adopt the spectra from PPPC4DMID for analyzing masses from 200 GeV to $<$ 100 TeV, and HDMSpectra for analyzing masses from 100 TeV to 30 PeV. We stress that the spectra from the two packages are consistent within the mass range from 1 TeV to 100 TeV [15], implying a low systematic uncertainty in the dark matter photon spectra. We also note that another dark matter photon-spectrum calculation method, CosmiXs (energy range: 5 GeV to 100 TeV) [24], also shows differences for indirect dark matter searches with IACTs that are smaller than the uncertainties due to the dwarf geometrical factors and the Poissonian fluctuations [e.g., 25].

The differential annihilation J-factor with respect to solid angle is obtained by integrating the dark matter density profile squared ( $\rho^{2}$ ) along the line-of-sight ( $l$ ) through the dSph for a given observation angle:

\frac{dJ}{d\Omega}=\int_{los}\rho^{2}(l,\Omega)\,dl

(2)

Using the commonly employed simplifying assumption that the dark matter density profiles are spherically symmetric, we can write the angular J-profile as a function of the observation angle ( $\theta$ ) from the center of the dSph as:

\frac{dJ}{d\theta}=2\pi\sin\theta\,\frac{dJ}{d\Omega}

(3)

The above J-profile can be integrated over a given angular range to determine the J-factor for a given dSph analysis.

For the dark matter density profile we utilize the Navarro-Frenk-White (NFW) model:

\rho(r)=\frac{\rho_{s}}{(r/r_{s})(1+(r/r_{s}))^{2}},\quad\textrm{if {\hbox{r<r_t}}, otherwise {\hbox{\rho}} = 0,}

(4)

where $\rho_{s}$ is a characteristic density, and $r_{s}$ and $r_{t}$ are a scale radius and a tidal truncation radius, respectively.

The set of NFW parameters ( $\rho_{s}$ , $r_{s}$ , and $r_{t}$ ) may be estimated from dynamical analyses of kinematic data of member stars. For this work we make use of the publication by [26] (Ando+20) and the associated Python package¹¹1https://github.com/shinichiroando/dwarf_params that provides the probability density functions (PDFs) for those parameters. With their Bayesian approach, they estimated the posterior distributions of the parameters by adopting either a physically motivated prior on a circular velocity of a satellite in a host subhalo [28] for the ultra-faint galaxies, or a well-known J-factor prior constrained by the large number of member stars [29] for the classical dSphs. With the use of physically informative priors, Ando+20 takes into account the current understanding of the role of structure formation on dark matter subhalos in dSphs, such as the extended Press-Schechter formalism combined with tidal effects on subhalo evolution [e.g., 30, 31, 28, 29]. Ando+20 cites the use of informative priors as important for improving the accuracy of the J-factor calculation for ultra-faint dwarfs with scarce kinematic data. As the majority of the targets included in our analysis are ultra-faint, we consider the use of Ando+20 J-factors appropriate. The impact of this choice is elaborated on when assessing the results.

This work presents the analysis of the complete VERITAS observational data set obtained on dSphs and covers the time period of 2007 - 2018. The VERITAS observation strategy evolved with time to focus on deep observations of those objects with the largest predicted J-factors, but to also obtain exposures on additional dSphs, while taking into account uncertainties on the published J-factors and the possibility that J-factors could have been underestimated. VERITAS devoted deep exposures to three dark-matter-dominated dSphs²²2We note that Draco II is one of the most dark-matter rich dSphs, but is not included in our deep-exposure sample because, during the data collection phase, Ando+20 had not yet been published.: Segue 1, Ursa Minor, and Ursa Major II, with more than 100 hours obtained on each target. To supplement these deep observations, survey observations were made on an additional 14 dSphs.

Dwarf	$\rho_{s,\text{ref}}$	$r_{s,\text{ref}}$	$r_{t,\text{ref}}$	$\rho_{s}$	$r_{s}$	$r_{t}$	log ${}_{10}J(\theta_{\rm cut})$	log ${}_{10}J(0.5^{o})$
	[M ${}_{\odot}/{\rm pc}^{3}$ ]	[pc]	[kpc]	[M ${}_{\odot}/{\rm pc}^{3}$ ]	[pc]	[kpc]	[GeV²/cm⁵]	[GeV²/cm⁵]
Boötes	0.045	425	3.50	7.0 ${}^{+11.2}_{-4.3}\times 10^{-2}$	3.1 ${}^{+2.9}_{-1.4}\times 10^{2}$	2.1 ${}^{+4.0}_{-1.2}$	$17.67_{-0.25}^{+0.26}$	$17.77_{-0.24}^{+0.23}$
Coma Berenices	0.067	397	2.15	7.9 ${}^{+12.6}_{-4.9}\times 10^{-2}$	3.5 ${}^{+3.9}_{-1.8}\times 10^{2}$	2.7 ${}^{+5.4}_{-1.7}$	$18.16_{-0.32}^{+0.29}$	$18.37_{-0.33}^{+0.30}$
CVn I	0.042	702	11.92	4.7 ${}^{+6.5}_{-2.7}\times 10^{-2}$	6.4 ${}^{+4.7}_{-2.6}\times 10^{2}$	6.0 ${}^{+6.9}_{-3.1}$	$17.31_{-0.12}^{+0.15}$	$17.38_{-0.11}^{+0.11}$
CVn II	0.062	381	2.28	8.4 ${}^{+14.0}_{-5.2}\times 10^{-2}$	3.1 ${}^{+3.7}_{-1.7}\times 10^{2}$	2.3 ${}^{+5.1}_{-1.5}$	$17.16_{-0.44}^{+0.38}$	$17.19_{-0.47}^{+0.37}$
Draco II	0.054	741	8.19	5.9 ${}^{+7.4}_{-3.2}\times 10^{-2}$	7.0 ${}^{+5.2}_{-3.0}\times 10^{2}$	7.7 ${}^{+8.5}_{-4.2}$	$19.14_{-0.25}^{+0.25}$	$19.49_{-0.25}^{+0.20}$
Hercules I	0.055	300	1.15	9.2 ${}^{+15.7}_{-6.0}\times 10^{-2}$	2.1 ${}^{+2.5}_{-1.1}\times 10^{2}$	1.2 ${}^{+2.8}_{-0.7}$	$16.92_{-0.36}^{+0.31}$	$16.93_{-0.39}^{+0.34}$
Leo I	0.046	946	11.40	5.7 ${}^{+6.7}_{-2.9}\times 10^{-2}$	8.2 ${}^{+5.2}_{-3.3}\times 10^{2}$	10.0 ${}^{+8.8}_{-4.8}$	$17.63_{-0.08}^{+0.09}$	$17.70_{-0.08}^{+0.07}$
Leo II	0.044	808	11.42	5.5 ${}^{+7.5}_{-3.1}\times 10^{-2}$	6.7 ${}^{+5.5}_{-2.9}\times 10^{2}$	6.8 ${}^{+8.1}_{-3.6}$	$17.48_{-0.09}^{+0.09}$	$17.54_{-0.10}^{+0.10}$
Leo IV	0.093	181	0.53	1.2 ${}^{+2.3}_{-0.8}\times 10^{-1}$	1.5 ${}^{+2.4}_{-0.9}\times 10^{2}$	0.8 ${}^{+2.4}_{-0.6}$	$16.57_{-0.64}^{+0.52}$	$16.56_{-0.66}^{+0.57}$
Leo V	0.102	181	0.87	1.1 ${}^{+2.3}_{-0.8}\times 10^{-1}$	1.7 ${}^{+3.0}_{-1.1}\times 10^{2}$	0.9 ${}^{+3.4}_{-0.7}$	$16.57_{-0.69}^{+0.60}$	$16.58_{-0.69}^{+0.60}$
Segue 1	0.086	344	3.18	1.1 ${}^{+1.7}_{-0.7}\times 10^{-1}$	2.9 ${}^{+3.7}_{-1.7}\times 10^{2}$	2.2 ${}^{+5.2}_{-1.6}$	$18.61_{-0.41}^{+0.34}$	$18.91_{-0.48}^{+0.39}$
Segue 2	0.186	65	0.25	1.9 ${}^{+3.6}_{-1.3}\times 10^{-1}$	0.6 ${}^{+1.3}_{-0.4}\times 10^{2}$	0.2 ${}^{+0.9}_{-0.2}$	$17.15_{-0.79}^{+0.48}$	$17.23_{-0.99}^{+0.58}$
Sextans I	0.043	645	11.07	5.0 ${}^{+7.0}_{-2.9}\times 10^{-2}$	5.6 ${}^{+4.6}_{-2.4}\times 10^{2}$	5.1 ${}^{+6.5}_{-2.9}$	$17.87_{-0.28}^{+0.28}$	$18.05_{-0.29}^{+0.25}$
Triangulum II	0.100	134	1.64	1.4 ${}^{+3.0}_{-1.0}\times 10^{-1}$	1.0 ${}^{+2.4}_{-0.7}\times 10^{2}$	0.5 ${}^{+2.1}_{-0.3}$	$17.56_{-0.77}^{+0.60}$	$17.65_{-0.90}^{+0.68}$
Ursa Major I	0.056	638	3.43	7.4 ${}^{+11.0}_{-4.4}\times 10^{-2}$	4.2 ${}^{+4.4}_{-2.2}\times 10^{2}$	3.6 ${}^{+6.5}_{-2.3}$	$18.05_{-0.24}^{+0.21}$	$18.19_{-0.25}^{+0.22}$
Ursa Major II	0.065	452	2.41	4.9 ${}^{+6.5}_{-2.7}\times 10^{-2}$	7.5 ${}^{+5.5}_{-3.2}\times 10^{2}$	8.1 ${}^{+8.5}_{-4.3}$	$18.43_{-0.38}^{+0.35}$	$18.79_{-0.48}^{+0.36}$
Ursa Minor	0.049	765	9.48	6.0 ${}^{+8.0}_{-3.4}\times 10^{-2}$	6.0 ${}^{+5.1}_{-2.7}\times 10^{2}$	5.9 ${}^{+7.7}_{-3.4}$	$18.22_{-0.24}^{+0.25}$	$18.47_{-0.22}^{+0.20}$

Table 1: Table of the NFW parameters of the 17 dSphs used in this work. The first three columns show the reference NFW profile parameters, which approximate the median of the distribution of viable J-profiles [described in Section II; adopted from 26]. The next three columns give the medians with 68% confidence intervals of the NFW profile parameters. The last two columns show the logarithm of the median J-factors with 68% confidence intervals integrated to

\theta_{\rm cut}

and 0.5^∘ from the centers of the galaxies, respectively. Each dSph was analyzed with an appropriate value of

\theta_{\rm cut}

: 0.089, 0.110, or 0.141 degrees (see Section III and Table 2).

In Table 1 we present the reference NFW parameters $\rho_{s}$ , $r_{s}$ , and $r_{t}$ , as well as the median and 68% confidence interval for the parameters, as provided by Ando+20. The median and interval are calculated from a random sample of 1,000 viable parameter sets for each dSph. The parameter sets constitute the 1,000 randomly-drawn profiles, providing the value and errors of the J-factor integrated out to an observation angle of 0.5^∘, as well as $\theta_{\rm cut}$ (see Section III). Since a strong correlation exists among three NFW parameters [ $\rho_{s}$ , $r_{s}$ , and $r_{t}$ ; see 33], the use of the median parameters does not guarantee the median value of the J-profile. To provide a reference J-profile for each dSph, we select the individual randomly-drawn profile that most closely matches the median J-profile calculated from 1,000 randomly drawn profiles for each object.

III VERITAS Instrument, Observations and Analysis

III.1 Instrument

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is a four 12 m-diameter IACT array, located at the Fred Lawrence Whipple Observatory in southern Arizona. It has been in full operation since 2007. VERITAS indirectly detects gamma rays in the energy range 85 GeV to $>30$ TeV through imaging the fast (ns) flashes of Cherenkov light produced by extensive air showers in the atmosphere. Each of the four VERITAS telescopes has a reflector comprised of 350 hexagonal mirrors, mounted on a Davies-Cotton optical support structure, which images the Cherenkov light onto a 499-pixel photomultiplier-tube (PMT) camera mounted in the focal plane [34]. Stereoscopic analysis of the images allows showers initiated by gamma rays to be preferentially selected over showers initiated by hadronic cosmic ray particles and the energies and directions of the primary gamma rays to be estimated.

The performance of VERITAS varies with energy; for 1 TeV photons VERITAS has an energy resolution of 17 $\%$ , an angular resolution of 0.08 deg (68% containment radius) and an effective gamma-ray collection area on the order of 10⁵ m². The VERITAS source location accuracy is $<$ 50 arcseconds. The sensitivity of VERITAS allows a statistically significant detection (5 standard deviations, 5 $\sigma$ , above background) of a point source with flux equivalent to 1 $\%$ that of the Crab Nebula in $\sim$ 25 hours [35].

VERITAS has undergone two major configuration changes since it began full four-telescope operations in 2007; in 2009 Telescope 1 was moved to a new location in order to improve sensitivity, while in 2012 the cameras were upgraded with higher quantum-efficiency PMTs and a revised trigger system, resulting in an improved low-energy response [36]. In addition to the upgrades discussed above, the VERITAS instrument response has changed with time due to mirror reflectivity changes and PMT gain changes. We perform regular calibration measurements [37] that are used to produce time-dependent instrument response functions (IRFs), which are used in this study since the data were taken over a time period of roughly 11 years (2007-2018).

III.2 Observations

All observations were made in the “wobble” observing mode, whereby the target being observed is offset in the field of view by 0.5^∘ to allow for simultaneous measurement of background at the same radial distance as the target from the center of the field of view [38]. Observing runs are typically 20-30 minutes in duration, and the direction of the offset is cycled through the four cardinal directions to reduce systematic effects.

Dwarf	$N_{ON}$	$N_{OFF}$	Exposure	$\theta_{\rm cut}$	Significance	$\Phi^{95\%}$
	[counts]	[total counts]	[hours]	[deg]	[ $\sigma$ ]	[10^-12cm^-2s^-1]
Boötes	446	2569	13.98	0.141	0.8	0.91
Coma Berenices	1122	6770	39.76	0.110	-0.2	0.66
CVn I	411	2430	9.72	0.141	0.3	0.98
CVn II	335	1822	8.14	0.141	1.6	1.69
Draco II	223	1335	8.02	0.141	0.0	1.62
Hercules I	369	2187	9.46	0.141	0.2	1.18
Leo I	196	1182	5.66	0.141	-0.1	1.98
Leo II	550	3275	11.31	0.141	0.2	0.97
Leo IV	7	65	0.48	0.141	-1.2	2.58
Leo V	33	218	1.38	0.141	-0.5	1.55
Segue 1	3070	18336	126.29	0.110	0.2	0.30
Segue 2	487	3000	12.51	0.110	-0.5	0.80
Sextans I	213	1262	7.45	0.141	0.2	0.95
Triangulum II	751	4870	29.51	0.110	-2.0	0.36
Ursa Major I	358	2073	6.63	0.141	0.6	0.99
Ursa Major II	2266	13855	212.32	0.089	-0.8	0.20
Ursa Minor	2253	13608	135.3	0.110	-0.3	0.28

Table 2: Table of the VERITAS observations of the 17 dwarf spheroidal galaxies. The first and second columns give the counts of gamma-like events in the ON and OFF regions, respectively. The next two columns show the exposure times and radial cuts defining the ON-region, respectively. The detection significance is given in the following column in terms of standard deviations above the background. The final column gives the 95% confidence level upper limit on the flux above 300 GeV, assuming a spectral index of -2.4. In all 17 cases, the background normalisation value

\alpha

is 0.167.

Overall, VERITAS obtained a total quality-selected exposure of 638 hours on the 17 dSphs. The targets and exposures are shown in Table 2. Note that four dSphs—namely Leo I, Leo II, Sextans, and Ursa Minor—are classical dSphs, and the rest are classified as ultra-faint dSphs. This study uses approximately 3–4 times more data than in the previous VERITAS study [total 230 hours on four dSphs; 11], and the previous data set is included in this work.

III.3 Analysis

VERITAS scientific results are produced and validated with two independent software packages: VEGAS [39] and EventDisplay [40]. Since the previous VERITAS dark matter publication, both packages have been improved by the addition of updated analysis techniques; specifically, boosted decision trees [BDTs; 41] in EventDisplay and an image template model [ITM; 42] in VEGAS. Studies have shown a $\sim$ 20 $\%$ - 25 $\%$ increase in sensitivity from the addition of BDTs and a $\sim$ 30 $\%$ increase in sensitivity from the addition of ITM analysis. In this work, we present the results from the EventDisplay package, although they are verified using VEGAS and an independent analysis pipeline for performing the dark matter likelihood analysis.

The gamma-ray-selection/hadron-rejection parameters chosen for this analysis are optimized for a moderate energy threshold, which differs from the previous VERITAS dark matter search where the parameters used were optimized for the lowest possible energy threshold. The decision to raise the energy threshold for this analysis is made in order to focus the dark matter search on the energy range where VERITAS is most sensitive, while avoiding systematic effects associated with deep exposures that affect data with less restrictive gamma/hadron selection.

While VERITAS has the capability to detect events with energies as low as 85 GeV, we apply a low-energy threshold to avoid events with poorly reconstructed energy and angular information. The energy threshold for each run is determined by taking into account the VERITAS IRFs.³³3Specifically, a run’s low-energy threshold is established as the higher of two values: the energy corresponding to 15% of the maximum point in the effective area curve or the energy at which the deviation between the true energy and the reconstructed energy is less than 20%. This threshold typically ranges from 200 GeV to 300 GeV, depending on the observing conditions. From a pilot study comparing results with and without applying the low-energy threshold, we found that considering a threshold effectively mitigates systematic effects without compromising sensitivity at energies where VERITAS is most effective in the dark matter search.

The data were analyzed to select candidate gamma-ray events, and a reflected-region background estimation was performed [44] to search for an excess of events from the direction of each target. The ON-source region was defined to be a circular region centered on the target. Multiple OFF regions are defined as circular regions of identical size to the ON region again offset by 0.5^∘ from the center of the field-of-view. Note that for EventDisplay, the number of OFF regions is set to six for each run, giving the relative exposure time between the ON and OFF regions to be 1/6 ( $\alpha\simeq 0.167$ ) for all dSph observations.

The radius of the ON-source region is normally chosen be comparable to the VERITAS point-spread-function (PSF) for a point source, or larger for a spatially extended target. All of the dSph objects in this study have angular extensions greater than the VERITAS PSF. However, when one expands the ON-source region size, there are fewer available regions for background estimation. Furthermore, for large data sets the analysis is more susceptible to systematic effects such as gradients in the number of events recorded across the cameras due to the varying zenith angle across the field of view. This effect is exacerbated for a larger ON-source region. The expected field-of-view significance distributions for an empty field is a Gaussian distribution of unit width centered at zero. Gradients across the camera can dramatically broaden the distribution, resulting in an unreliable assessment of the signal significance at the target position.

Prior to this analysis being conducted, a study was performed to optimize the radial cut used to define the ON-source region for dSphs with deep and shallow exposures to maximize the J-factor enclosed while preserving an acceptable field-of-view significance distribution on a control data set. Due to the large number of targets, all with differing exposures, we considered three possible exposure times, corresponding to targets with 16 hours of exposure, 16-150 hours of exposure, and more than 150 hours of exposure. We iteratively increased the radial cut defining the ON-source region, and selected the maximum radial cut for which no systematic effects in the field-of-view significance distributions were visible on a test data set.

Based on the results of the study on the control sample, we use radial cuts to define the ON and OFF regions of 0.089^∘ ( $\theta^{2}=0.008~\text{deg}^{2}$ ) for the deepest exposure target (Ursa Major II), 0.110^∘ ( $\theta^{2}=0.012~\text{deg}^{2}$ ) for the remaining targets with more than 16 hours of exposure, and 0.141^∘ ( $\theta^{2}=0.02~\text{deg}^{2}$ ) for targets with less than 16 hours of exposure (see Table 2). Note that the data for Segue 2 were analyzed with a radial cut appropriate for a longer exposure, as its raw exposure of more than 16 hours becomes less than 16 hours after application of time cuts to remove periods with weather or hardware issues.

IV Maximum Likelihood Estimation

IV.1 Predicted signal

For a given VERITAS analysis we can calculate the number of gamma-ray events expected ( $\mathcal{S}$ ) for a given data set by taking the predicted flux (Equation 1) and folding it with the VERITAS IRFs:

\mathcal{S}=\frac{\left<\sigma v\right>T_{obs}}{8\pi M^{2}_{\chi}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{dN}{dE^{\prime}}J(E^{\prime})A(E^{\prime})D(E|E^{\prime})dE^{\prime}dE,

(5)

where $T_{obs}$ is the total exposure time, $J(E^{\prime})$ is the J-profile convolved with the VERITAS energy-dependent PSF and integrated to the appropriate value of $\theta$ , $A(E^{\prime})$ is the VERITAS effective area, and $D(E|E^{\prime})$ is the energy migration matrix, in which $E$ is the reconstructed photon energy and $E^{\prime}$ is the true photon energy. Time-averaged IRFs are used in this analysis for each dSph. Using a maximum likelihood analysis we can test for a significant excess of gamma rays and place constraints on the value of $\left<\sigma v\right>$ as a function of dark matter particle mass.

IV.2 Likelihood

The use of the maximum likelihood estimation (MLE) method has been shown to maximize the sensitivity to gamma-rays produced from dark matter interactions [e.g. 45]. We adopt the MLE method as described in [46]. Source analysis with an IACT such as VERITAS is performed by comparing the observed number of ON-region (source region) events with the observed number of OFF-region (background region) events. As the ON and OFF counts are Poisson-distributed, we construct the full likelihood function by multiplying two Poisson likelihood functions with the model-dependent likelihood function $P_{i}(E_{i}|M_{\chi},\langle\sigma v\rangle)$ ,

\mathcal{L}=\frac{(\mathcal{S}+\alpha\mathcal{B})^{N_{on}}e^{-(\mathcal{S}+\alpha\mathcal{B})}}{N_{on}!}\frac{\mathcal{B}^{N_{off}}e^{-\mathcal{B}}}{N_{off}!}\prod^{N_{on}}_{i=1}P_{i}(E_{i}|M_{\chi},\langle\sigma v\rangle),

(6)

where $\alpha$ is a background normalisation factor that accounts for the ratio of the number of ON regions to OFF regions, $\mathcal{B}$ is a nuisance parameter that describes the total expected number of background counts from multiple OFF regions, and $\mathcal{S}$ is the expected number of events from dark matter annihilation at a given dark matter mass and velocity-weighted cross section within the ON region (Equation 5). Finally, $P_{i}(E_{i}|M_{\chi},\langle\sigma v\rangle)$ is the predicted energy-dependent counts distribution in the ON region, where these counts can be from the dark matter self-annihilation or cosmic-ray background. The distribution is given by

P_{i}(E_{i}|M_{\chi},\left<\sigma v\right>)=\frac{\mathcal{S}p_{s}(E_{i})+\alpha\mathcal{B}p_{b}(E_{i})}{\mathcal{S}+\alpha\mathcal{B}},

(7)

where $p_{s}$ is from the normalised signal spectrum at a given energy, and $p_{b}$ is generated from a normalised histogram of the energies of all background events. In other words, they are the probability density functions for the dark matter signal and the background, respectively.

We can then maximize the log likelihood function, neglecting constant terms,

\log{\mathcal{L}(\langle\sigma v\rangle;\mathbf{\mathcal{B}})}=N_{off}\log{\mathcal{B}}-\mathcal{S}-(1+\alpha)\mathcal{B}+\sum_{i=1}^{N_{on}}{\log{(\alpha\mathcal{B}\,p_{b}(E_{i})+\mathcal{S}\,p_{s}(E_{i}))}},

(8)

with respect to $\left<\sigma v\right>$ and $\mathcal{B}$ , giving us a constraint on the dark matter velocity-weighted annihilation cross-section.

For the joint-fit MLE analysis, we combine data from the 17 dSphs to maximize statistical power using the joint likelihood function

\log{\mathcal{L}_{\rm joint}(\langle\sigma v\rangle;\mathbf{\mathcal{B}})}=\sum_{i=1}^{N_{\text{target}}}\log{\mathcal{L}(\langle\sigma v\rangle;\mathcal{B}_{i})}.

(9)

To determine the significance of the dark matter signal over the background, we compare two likelihoods (null and alternative hypotheses) using the following equation,

\lambda=-2\log\left(\frac{\mathcal{L}_{\mathcal{S}\equiv 0}}{\mathcal{L}_{\mathcal{S}\neq 0}}\right)\!.

(10)

When the resulting value of $\lambda$ is below the detection threshold of 25 (equivalent to $5\sigma$ detection), we calculate an upper limit (UL) on the velocity-weighted annihilation cross section by utilizing the likelihood profile. The UL is obtained by computing $\Delta\log\mathcal{L}$ of 1.35 compared to the likelihood maximum, corresponding to the one-sided 95% confidence level.

Refer to caption — Figure 1: Profile likelihood for a joint-fit analysis with 17 dSphs for a dark matter particle mass of 1 TeV (left) and the VERITAS 95% upper limit on the dark matter cross section (right), considering the $\tau^{+}\tau^{-}$ annihilation channel. The solid black line shows the joint-fit analysis result, and the dashed gray lines are from the individual dSph analyses.

V Results

V.1 Non-detection with Li & Ma analysis

The significance of a signal above background in the ON region is estimated based on the Li & Ma method [47], using the counts in the ON and OFF regions and the ratio of areas between the ON and OFF regions ( $\alpha$ =0.167 for this analysis). A summary of the counts, the detection significance, and the 95% confidence level upper limits on the flux is shown in Table 2. No dSph shows a significant signal.

V.2 Non-detection with maximum likelihood estimation

We test nine annihilation channels ( $e^{+}e^{-}$ , $\mu^{+}\mu^{-}$ , $\tau^{+}\tau^{-}$ , $t\bar{t}$ , $b\bar{b}$ , $W^{+}W^{-}$ , $ZZ$ , $\gamma\gamma$ , and $\nu_{e}\bar{\nu}_{e}$ ) using the MLE analysis and find that there is no evidence of dark matter annihilation signals from the 17 dSph observations; i.e., the flux and energy spectrum of observed events from the source region is consistent with that of the background fluctuations. The following sections describe our constraint on the velocity-weighted annihilation cross section of dark matter in various aspects. Note that in the case of the $\nu_{e}\bar{\nu}_{e}$ annihilation channel, production of final-state gamma rays is enabled by radiation of a $Z$ boson by an off-shell neutrino at sufficiently high energies, or decay of an off-shell neutrino to a $W$ boson and an electron or positron.

V.3 Upper limits on the dark matter velocity-weighted annihilation cross section

Fig. 1 describes the method of obtaining the dark matter annihilation cross section from the joint-fit MLE analysis and the results. The left panel shows an example of the profile likelihood from individual MLE analyses (dotted or dashed lines) and the joint-fit MLE analysis (solid black line) for the dark matter mass of 1 TeV. In this case, the joint-fit profile likelihood has $\Delta\log(\mathcal{L})$ of 1.35 at $\langle\sigma v\rangle=1.76\times 10^{-24}$ , corresponding to the one-sided 95% confidence UL. This process is repeated for other dark matter masses so that we can get the UL curve as seen in the right panel of Fig. 1. Generally, the joint-fit result provides a stronger constraint on the velocity-weighted dark matter annihilation cross section compared to those from the individual fits.

Since the J-factors have significant uncertainties, we compute for every dSph a set of ULs for each annihilation channel using 300 possible J-factors. Each J-factor is randomly generated by drawing the three NFW parameters from the posterior distributions provided by Ando+20 (see Table 2) and integrating over the line of sight and solid angle. From the 300 realizations, we obtain the median as well as 68% and 95% containments for ULs on the dark matter annihilation cross section, which reflects the systematic uncertainty due to the limited understanding of the dark matter distribution.

Fig. 2 shows the UL band from the joint-fit analysis, overplotted with UL curves from the three dSphs with the deepest exposures (Segue 1, Ursa Major II, and Ursa Minor) as a reference (see also Appendix A). The red dotted-dashed line corresponds to the theoretical expectation of the velocity-weighted annihilation cross section of thermal-relic dark matter, $\langle\sigma v\rangle\sim 2.4\times 10^{-26}$ cm³/s, extending to the unitarity bound at $\sim$ 100 TeV [17], except for the $\gamma\gamma$ channel. For the loop-suppressed $\gamma\gamma$ channel, we use $\langle\sigma v\rangle\sim 1\times 10^{-28}$ cm³/s [48]. Note that for the individual dSph analyses, we take the median of the J-factor distribution and compute a single UL. In general, Segue 1 primarily influences the joint-fit result in the low mass range (up to about 10 TeV), while the joint-fit result in the high mass range is predominantly determined by the observations of Ursa Major II. Ursa Minor and other dSphs have minimal impact on the joint-fit result.

The discontinuity observed in the limit for $M_{\chi}$ at $\sim 100$ TeV for the $\gamma\gamma$ channel is expected. We do not analyze photon events with energies above $100$ TeV, as they are beyond the energy range in which events can be reliably reconstructed. Hence, for $M_{\chi}>100$ TeV we only consider the continuum $\gamma\gamma$ signal that produces $\lesssim$ 100 TeV gamma rays, rather than the more easily detectable line signature located above 100 TeV.

V.4 Comparison with the null hypothesis

In addition to calculating “observed” ULs, we calculate “expected” ULs assuming the background-only (null) hypothesis. Simulated ON regions are generated through a procedure where events (and their corresponding energies) are randomly selected (with replacement) from the observed OFF-region events, which are assumed to be pure background. The events are sampled according to a Poisson distribution, with the mean equal to the observed number of OFF-region events scaled by the ratio of the areas of the ON and OFF regions; i.e., $N_{\rm on,sim}={\rm Pois}(\alpha N_{\rm off,obs})$ . Simulated OFF-regions are obtained in the same manner, considering Poisson fluctuations through random sampling. We repeat this process 300 times for each channel, resulting in an expected upper limit band. The width of this band is determined by the magnitude of the Poisson fluctuations of ON and OFF regions. We use the same J-profile for both the expected and observed limits (Section V.3)

Fig. 3 shows the comparison between the expected UL band (orange) and the observed UL (blue solid line). For each annihilation channel, the expected upper limit band shows the 68% containment (orange shaded region) and 95% containment (red hatched region). The observed upper limits remain consistent with the expected upper limits within the 95% confidence level across all annihilation channels. This result indicates that a dark matter annihilation signal has not been observed, while also quantifying the influence of statistical uncertainty on the derived ULs.

V.5 Ultra-heavy dark matter search

For a dark matter mass exceeding 194 TeV (UHDM regime), the observed data did not show significant deviation between the observed and expected limits (Fig. 3). In Fig. 4, we interpret this non-detection in terms of theoretical benchmark models [for details on the benchmark models, see 15]. The left panel of Fig. 4 shows both the unitarity limit for a point-like particle with angular momentum J = 0 (partial-wave unitarity limit; $\langle\sigma v\rangle_{unitarity}\propto 1/v_{rel}$ ) and a set of unitarity bounds for composite dark matter particles of different radii, $\langle\sigma v\rangle_{unitarity}(v_{rel},R)$ . Note that, when computing the partial-wave and composite unitarity bounds, we adopt $v_{rel}/c=2\times 10^{-5}$ , where $v_{rel}$ represents the relative velocity between dark matter particles in dSph galaxies [50, 51]. This velocity $v_{rel}$ is much slower than that of the thermal-relic dark matter particles in the early universe. While the limits from this dataset are unable to probe below the limit for a point-like particle, it is possible to rule out models with large dark matter particle radii. The right panel shows that we are able to exclude a certain range of values for the radius of the composite particle. If the mass of the dark matter particle is less than 1 PeV, its radius should be smaller than the proton charge radius, 0.84 fm [49].

V.6 Comparison with previous works

We compare the derived ULs with those from the previous VERITAS work [VERITAS 17; 11]. As described in Section III, in this work, we analyzed a larger dataset with improved analysis techniques. Another significant difference in this study is the utilization of updated J-factors from Ando+20 (Section II), whereas the previous VERITAS work relied on J-factors estimated using uniform (or uninformative) priors [GS+15; 52]. We note that the J-factors calculated using the NFW parameters of Ando+20 result in a lower J-factor on average than the parameters used in the previous VERITAS study using the NFW parameters from GS+15. This is discussed in more detail in Appendix B. To quantify the improvement, we compute ULs for Segue 1 for the published dataset and current analysis tools, and compare against the published Segue 1 ULs. We consider the $b^{+}b^{-}$ and $\tau^{+}\tau^{-}$ annihilation channels.

Fig. 5 shows UL curves from the extended Segue 1 observation for the two J-factors, GS+15 (orange dotted line) and Ando+20 (blue solid line), and the UL curve from VERITAS 17 (black dashed line). The limits indicate that our enhanced techniques, such as BDTs, ITM, and an optimized $\theta^{2}$ cut, along with the increased exposure time, contribute to constraints on the dark matter annihilation cross-section that are more stringent by a factor of 1.7 to 5.1. This is evident when we compare the limits obtained from the same J factor, represented by the orange dotted line (this work) and the black dashed line (the previous work) in Fig. 5. Note that while extending the exposure time on Segue 1 by 37% (126 h/92 h) lowers the upper limit slightly, this enhancement is negligible compared to the significant improvement achieved by the enhanced techniques, which lower the limits by an average factor of 3. In contrast, utilizing Ando+20 leads to a less restrictive limit on the dark matter annihilation cross section, as seen by the upper limit (blue solid line) increasing by approximately three times compared to the UL curve with GS+15 (orange dotted line). Note that the J-factor of Segue 1 from Ando+20 is smaller than that of GS+15 by a factor of about 3 (see Appendix B). Overall, taking into account both the positive and negative effects, we arrive at a UL curve similar to that of VERITAS 17.

VI Discussion and Conclusions

Fig. 6 shows our upper limits for two annihilation channels (black lines) in comparison to results utilizing dSph observations from the Fermi-LAT [blue dashed; 2], MAGIC [orange dashed; 9], H.E.S.S. [greeen dashed; 7], and HAWC [red dashed; 14] collaborations. Using the Ando+20 J-factors (black solid), this work demonstrates results on par with the most stringent results from other collaborations. Note that the results from the other collaborations were derived from a range of sources for the J-factor, including but not limited to GS+15.

Due to the significant influence of the J-factor on the constraint of the dark matter velocity-weighted annihilation cross-section, it is crucial to utilize an accurate J-factor when conducting indirect dark matter searches with gamma rays. In this work, we adopt the recent J-factor study by Ando+20, which employs physically motivated informative priors (Section II). We stress that this is in contrast to previous J-factor estimates (e.g., GS+15, used in the previous VERITAS work) which rely on uniform priors or neglect relevant physical information. If we account for the differences in the J factor between GS+15 and Ando+20, our upper limits are expected to be lowered by a factor of about three (see Appendix B). This effect arises from the inverse proportionality between the upper limit and the J-factor, $\langle\sigma v\rangle_{95}\propto 1/J(\Delta\Omega)$ .

In summary, this paper conducts an indirect search for a dark matter annihilation signal in gamma rays. We analyzed 638 hours of observations taken by the VERITAS array, spanning from 2007 to 2018. Seventeen dSphs were observed, employing two observational strategies: deep exposures on dark matter-dense dSphs and survey observations of dark matter-sparse dSphs. Our search targeted the detection of the final state gamma rays resulting from nine annihilation channels. Although no significant signals were detected, we derived upper limits on the dark matter velocity-weighted annihilation cross section for the nine annihilation channels. These limits were obtained using a joint-fit MLE analysis and cover a range of dark matter particle masses, extending from 200 GeV to 30 PeV. With the extended data set and the improved techniques, we were able to provide competitive ULs above about 1 TeV, even using a set of smaller but physically motivated J-factors. In addition, we examined the ultra-heavy dark matter scenario and constrained the radius of a composite ultra-heavy dark matter particle.

VII Acknowledgments

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. Tak is supported by the National Research Foundation of Korea (NRF) grant, No.2020R1A2C3011091, and No.2021M3F7A1084525, and RS-2024-00343729, funded by the Korea government (MSIT).

References

Aghanim et al. [2020] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
Ackermann et al. [2015a] M. Ackermann et al., Phys. Rev. Lett. 115, 231301 (2015a), arXiv:1503.02641 [astro-ph.HE] .
Ackermann et al. [2015b] M. Ackermann, M. Ajello, A. Albert, B. Anderson, W. B. Atwood, L. Baldini, G. Barbiellini, D. Bastieri, R. Bellazzini, E. Bissaldi, R. D. Blandford, E. D. Bloom, R. Bonino, E. Bottacini, T. J. Brandt, J. Bregeon, P. Bruel, R. Buehler, S. Buson, G. A. Caliandro, R. A. Cameron, R. Caputo, M. Caragiulo, P. A. Caraveo, C. Cecchi, E. Charles, A. Chekhtman, J. Chiang, G. Chiaro, S. Ciprini, R. Claus, J. Cohen-Tanugi, J. Conrad, A. Cuoco, S. Cutini, F. D’Ammando, A. de Angelis, F. de Palma, R. Desiante, S. W. Digel, L. Di Venere, P. S. Drell, A. Drlica-Wagner, C. Favuzzi, S. J. Fegan, A. Franckowiak, Y. Fukazawa, S. Funk, P. Fusco, F. Gargano, D. Gasparrini, N. Giglietto, F. Giordano, M. Giroletti, G. Godfrey, G. A. Gomez-Vargas, I. A. Grenier, J. E. Grove, S. Guiriec, M. Gustafsson, J. W. Hewitt, A. B. Hill, D. Horan, G. Jóhannesson, R. P. Johnson, M. Kuss, S. Larsson, L. Latronico, J. Li, L. Li, F. Longo, F. Loparco, M. N. Lovellette, P. Lubrano, D. Malyshev, M. Mayer, M. N. Mazziotta, J. E. McEnery, P. F. Michelson, T. Mizuno, A. A. Moiseev, M. E. Monzani, A. Morselli, S. Murgia, E. Nuss, T. Ohsugi, M. Orienti, E. Orlando, J. F. Ormes, D. Paneque, M. Pesce-Rollins, F. Piron, G. Pivato, S. Rainò, R. Rando, M. Razzano, A. Reimer, T. Reposeur, S. Ritz, M. Sánchez-Conde, A. Schulz, C. Sgrò, E. J. Siskind, F. Spada, G. Spandre, P. Spinelli, H. Tajima, H. Takahashi, J. B. Thayer, L. Tibaldo, D. F. Torres, G. Tosti, E. Troja, G. Vianello, M. Werner, B. L. Winer, K. S. Wood, M. Wood, G. Zaharijas, and S. Zimmer, Phys. Rev. D 91, 122002 (2015b).
Albert et al. [2017] A. Albert, B. Anderson, K. Bechtol, A. Drlica-Wagner, M. Meyer, M. Sánchez-Conde, L. Strigari, M. Wood, T. M. C. Abbott, F. B. Abdalla, A. Benoit-Lévy, G. M. Bernstein, R. A. Bernstein, E. Bertin, D. Brooks, D. L. Burke, A. Carnero Rosell, M. Carrasco Kind, J. Carretero, M. Crocce, C. E. Cunha, C. B. D’Andrea, L. N. da Costa, S. Desai, H. T. Diehl, J. P. Dietrich, P. Doel, T. F. Eifler, A. E. Evrard, A. Fausti Neto, D. A. Finley, B. Flaugher, P. Fosalba, J. Frieman, D. W. Gerdes, D. A. Goldstein, D. Gruen, R. A. Gruendl, K. Honscheid, D. J. James, S. Kent, K. Kuehn, N. Kuropatkin, O. Lahav, T. S. Li, M. A. G. Maia, M. March, J. L. Marshall, P. Martini, C. J. Miller, R. Miquel, E. Neilsen, B. Nord, R. Ogando, A. A. Plazas, K. Reil, A. K. Romer, E. S. Rykoff, E. Sanchez, B. Santiago, M. Schubnell, I. Sevilla-Noarbe, R. C. Smith, M. Soares-Santos, F. Sobreira, E. Suchyta, M. E. C. Swanson, G. Tarle, V. Vikram, A. R. Walker, R. H. Wechsler, Fermi-LAT Collaboration, and DES Collaboration, Astrophys. J. 834, 110 (2017), arXiv:1611.03184 [astro-ph.HE] .
Abdalla et al. [2018] H. Abdalla et al., JCAP 2018 (11), 037, arXiv:1810.00995 [astro-ph.HE] .
Abdallah et al. [2018] H. Abdallah, A. Abramowski, F. Aharonian, F. Ait Benkhali, E. O. Angüner, M. Arakawa, M. Arrieta, P. Aubert, M. Backes, A. Balzer, M. Barnard, Y. Becherini, J. Becker Tjus, D. Berge, S. Bernhard, K. Bernlöhr, R. Blackwell, M. Böttcher, C. Boisson, J. Bolmont, S. Bonnefoy, P. Bordas, J. Bregeon, F. Brun, P. Brun, M. Bryan, M. Büchele, T. Bulik, M. Capasso, S. Caroff, A. Carosi, J. Carr, S. Casanova, M. Cerruti, N. Chakraborty, R. C. G. Chaves, A. Chen, J. Chevalier, S. Colafrancesco, B. Condon, J. Conrad, I. D. Davids, J. Decock, C. Deil, J. Devin, P. deWilt, L. Dirson, A. Djannati-Ataï, W. Domainko, A. Donath, L. O. Drury, K. Dutson, J. Dyks, T. Edwards, K. Egberts, P. Eger, G. Emery, J.-P. Ernenwein, S. Eschbach, C. Farnier, S. Fegan, M. V. Fernandes, A. Fiasson, G. Fontaine, A. Förster, S. Funk, M. Füßling, S. Gabici, Y. A. Gallant, T. Garrigoux, F. Gaté, G. Giavitto, B. Giebels, D. Glawion, J. F. Glicenstein, D. Gottschall, M.-H. Grondin, J. Hahn, M. Haupt, J. Hawkes, G. Heinzelmann, G. Henri, G. Hermann, J. A. Hinton, W. Hofmann, C. Hoischen, T. L. Holch, M. Holler, D. Horns, A. Ivascenko, H. Iwasaki, A. Jacholkowska, M. Jamrozy, M. Janiak, D. Jankowsky, F. Jankowsky, M. Jingo, L. Jouvin, I. Jung-Richardt, M. A. Kastendieck, K. Katarzyński, M. Katsuragawa, U. Katz, D. Kerszberg, D. Khangulyan, B. Khélifi, J. King, S. Klepser, D. Klochkov, W. Kluźniak, N. Komin, K. Kosack, S. Krakau, M. Kraus, P. P. Krüger, H. Laffon, G. Lamanna, J. Lau, J.-P. Lees, J. Lefaucheur, A. Lemière, M. Lemoine-Goumard, J.-P. Lenain, E. Leser, R. Liu, T. Lohse, M. Lorentz, R. López-Coto, I. Lypova, D. Malyshev, V. Marandon, A. Marcowith, C. Mariaud, R. Marx, G. Maurin, N. Maxted, M. Mayer, P. J. Meintjes, M. Meyer, A. M. W. Mitchell, R. Moderski, M. Mohamed, L. Mohrmann, K. Morå, E. Moulin, T. Murach, S. Nakashima, M. de Naurois, H. Ndiyavala, F. Niederwanger, J. Niemiec, L. Oakes, P. O’Brien, H. Odaka, S. Ohm, M. Ostrowski, I. Oya, M. Padovani, M. Panter, R. D. Parsons, N. W. Pekeur, G. Pelletier, C. Perennes, P.-O. Petrucci, B. Peyaud, Q. Piel, S. Pita, V. Poireau, H. Poon, D. Prokhorov, H. Prokoph, G. Pühlhofer, M. Punch, A. Quirrenbach, S. Raab, R. Rauth, A. Reimer, O. Reimer, M. Renaud, R. de los Reyes, F. Rieger, L. Rinchiuso, C. Romoli, G. Rowell, B. Rudak, C. B. Rulten, V. Sahakian, S. Saito, D. A. Sanchez, A. Santangelo, M. Sasaki, M. Schandri, R. Schlickeiser, F. Schüssler, A. Schulz, U. Schwanke, S. Schwemmer, M. Seglar-Arroyo, M. Settimo, A. S. Seyffert, N. Shafi, I. Shilon, K. Shiningayamwe, R. Simoni, H. Sol, F. Spanier, M. Spir-Jacob, L. Stawarz, R. Steenkamp, C. Stegmann, C. Steppa, I. Sushch, T. Takahashi, J.-P. Tavernet, T. Tavernier, A. M. Taylor, R. Terrier, L. Tibaldo, D. Tiziani, M. Tluczykont, C. Trichard, M. Tsirou, N. Tsuji, R. Tuffs, Y. Uchiyama, J. van der Walt, C. van Eldik, C. van Rensburg, B. van Soelen, G. Vasileiadis, J. Veh, C. Venter, A. Viana, P. Vincent, J. Vink, F. Voisin, H. J. Völk, T. Vuillaume, Z. Wadiasingh, S. J. Wagner, P. Wagner, R. M. Wagner, R. White, A. Wierzcholska, P. Willmann, A. Wörnlein, D. Wouters, R. Yang, D. Zaborov, M. Zacharias, R. Zanin, A. A. Zdziarski, A. Zech, F. Zefi, A. Ziegler, J. Zorn, and N. Żywucka (H.E.S.S. Collaboration), Phys. Rev. Lett. 120, 201101 (2018).
Abdallah et al. [2020] H. Abdallah et al., Phys. Rev. D 102, 062001 (2020), arXiv:2008.00688 [astro-ph.HE] .
Abdalla et al. [2022] H. Abdalla, F. Aharonian, F. A. Benkhali, E. O. Angüner, C. Armand, H. Ashkar, M. Backes, V. Baghmanyan, V. B. Martins, R. Batzofin, Y. Becherini, D. Berge, K. Bernlöhr, B. Bi, M. Böttcher, J. Bolmont, M. d. B. de Lavergne, R. Brose, F. Brun, F. Cangemi, S. Caroff, M. Cerruti, T. Chand, A. Chen, G. Cotter, J. D. Mbarubucyeye, J. Devin, A. Djannati-Ataï, A. Dmytriiev, V. Doroshenko, K. Egberts, A. Fiasson, G. F. de Clairfontaine, G. Fontaine, S. Funk, S. Gabici, G. Giavitto, D. Glawion, J. F. Glicenstein, M.-H. Grondin, J. A. Hinton, W. Hofmann, T. L. Holch, M. Holler, D. Horns, Z. Huang, M. Jamrozy, F. Jankowsky, E. Kasai, K. Katarzyński, U. Katz, B. Khélifi, W. Kluźniak, N. Komin, K. Kosack, D. Kostunin, G. Lamanna, M. Lemoine-Goumard, J.-P. Lenain, F. Leuschner, T. Lohse, A. Luashvili, I. Lypova, J. Mackey, D. Malyshev, D. Malyshev, V. Marandon, P. Marchegiani, G. Martí-Devesa, R. Marx, G. Maurin, M. Meyer, A. Mitchell, R. Moderski, A. Montanari, E. Moulin, J. Muller, M. de Naurois, J. Niemiec, A. P. Noel, S. Ohm, L. Olivera-Nieto, E. d. O. Wilhelmi, M. Ostrowski, S. Panny, M. Panter, R. D. Parsons, G. Peron, V. Poireau, H. Prokoph, G. Pühlhofer, M. Punch, A. Quirrenbach, P. Reichherzer, A. Reimer, O. Reimer, M. Renaud, F. Rieger, G. Rowell, B. Rudak, H. R. Ricarte, E. Ruiz-Velasco, V. Sahakian, H. Salzmann, A. Santangelo, M. Sasaki, F. Schüssler, H. M. Schutte, U. Schwanke, M. Senniappan, J. N. S. Shapopi, H. Sol, A. Specovius, S. Spencer, L. Stawarz, C. Stegmann, S. Steinmassl, C. Steppa, T. Takahashi, T. Tanaka, R. Terrier, C. Thorpe-Morgan, M. Tluczykont, M. Tsirou, N. Tsuji, Y. Uchiyama, C. van Eldik, J. Veh, J. Vink, S. J. Wagner, R. White, A. Wierzcholska, Y. W. Wong, M. Zacharias, D. Zargaryan, A. A. Zdziarski, A. Zech, S. J. Zhu, S. Zouari, and N. Żywucka (H.E.S.S. Collaboration), Phys. Rev. Lett. 129, 111101 (2022).
Acciari et al. [2022] V. A. Acciari et al., Physics of the Dark Universe 35, 100912 (2022), arXiv:2111.15009 [astro-ph.HE] .
Abe et al. [2023] H. Abe, S. Abe, V. A. Acciari, T. Aniello, S. Ansoldi, L. A. Antonelli, A. Arbet Engels, C. Arcaro, M. Artero, K. Asano, D. Baack, A. Babić, A. Baquero, U. Barres de Almeida, J. A. Barrio, I. Batković, J. Baxter, J. Becerra González, W. Bednarek, E. Bernardini, M. Bernardos, A. Berti, J. Besenrieder, W. Bhattacharyya, C. Bigongiari, A. Biland, O. Blanch, G. Bonnoli, i. c. v. Bošnjak, I. Burelli, G. Busetto, R. Carosi, M. Carretero-Castrillo, G. Ceribella, Y. Chai, A. Chilingarian, S. Cikota, E. Colombo, J. L. Contreras, J. Cortina, S. Covino, G. D’Amico, V. D’Elia, P. Da Vela, F. Dazzi, A. De Angelis, B. De Lotto, A. Del Popolo, M. Delfino, J. Delgado, C. Delgado Mendez, D. Depaoli, F. Di Pierro, L. Di Venere, E. Do Souto Espiñeira, D. Dominis Prester, A. Donini, D. Dorner, M. Doro, D. Elsaesser, G. Emery, V. Fallah Ramazani, L. Fariña, A. Fattorini, L. Font, C. Fruck, S. Fukami, Y. Fukazawa, R. J. García López, M. Garczarczyk, S. Gasparyan, M. Gaug, J. G. Giesbrecht Paiva, N. Giglietto, F. Giordano, P. Gliwny, N. Godinović, J. G. Green, D. Green, D. Hadasch, A. Hahn, T. Hassan, L. Heckmann, J. Herrera, D. Hrupec, M. Hütten, R. Imazawa, T. Inada, R. Iotov, K. Ishio, I. Jiménez Martínez, J. Jormanainen, D. Kerszberg, Y. Kobayashi, H. Kubo, J. Kushida, A. Lamastra, D. Lelas, F. Leone, E. Lindfors, L. Linhoff, S. Lombardi, F. Longo, R. López-Coto, M. López-Moya, A. López-Oramas, S. Loporchio, A. Lorini, E. Lyard, B. Machado de Oliveira Fraga, P. Majumdar, M. Makariev, G. Maneva, N. Mang, M. Manganaro, S. Mangano, K. Mannheim, M. Mariotti, M. Martínez, A. Mas Aguilar, D. Mazin, S. Menchiari, S. Mender, S. Mićanović, D. Miceli, T. Miener, J. M. Miranda, R. Mirzoyan, E. Molina, H. A. Mondal, A. Moralejo, D. Morcuende, V. Moreno, T. Nakamori, C. Nanci, L. Nava, V. Neustroev, M. Nievas Rosillo, C. Nigro, K. Nilsson, K. Nishijima, T. Njoh Ekoume, K. Noda, S. Nozaki, Y. Ohtani, T. Oka, J. Otero-Santos, S. Paiano, M. Palatiello, D. Paneque, R. Paoletti, J. M. Paredes, L. Pavletić, M. Persic, M. Pihet, F. Podobnik, P. G. Prada Moroni, E. Prandini, G. Principe, C. Priyadarshi, I. Puljak, W. Rhode, M. Ribó, J. Rico, C. Righi, A. Rugliancich, N. Sahakyan, T. Saito, S. Sakurai, K. Satalecka, F. G. Saturni, B. Schleicher, K. Schmidt, F. Schmuckermaier, J. L. Schubert, T. Schweizer, J. Sitarek, V. Sliusar, D. Sobczynska, A. Spolon, A. Stamerra, J. Strišković, D. Strom, M. Strzys, Y. Suda, T. Surić, M. Takahashi, R. Takeishi, F. Tavecchio, P. Temnikov, K. Terauchi, T. Terzić, M. Teshima, L. Tosti, S. Truzzi, A. Tutone, S. Ubach, J. van Scherpenberg, M. Vazquez Acosta, S. Ventura, V. Verguilov, I. Viale, C. F. Vigorito, V. Vitale, I. Vovk, R. Walter, M. Will, C. Wunderlich, T. Yamamoto, D. Zarić, N. Hiroshima, and K. Kohri (MAGIC Collaboration), Phys. Rev. Lett. 130, 061002 (2023).
Archambault et al. [2017] S. Archambault et al., Phys. Rev. D 95, 082001 (2017), arXiv:1703.04937 [astro-ph.HE] .
Albert et al. [2018] A. Albert et al., Astrophys. J. 853, 154 (2018), arXiv:1706.01277 [astro-ph.HE] .
Abeysekara et al. [2018] A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, R. Arceo, J. C. Arteaga-Velázquez, D. Avila Rojas, H. A. Ayala Solares, A. Becerril, E. Belmont-Moreno, S. Y. BenZvi, A. Bernal, C. Brisbois, K. S. Caballero-Mora, T. Capistrán, A. Carramiñana, S. Casanova, M. Castillo, U. Cotti, J. Cotzomi, C. De León, E. De la Fuente, R. Diaz Hernandez, B. L. Dingus, M. A. DuVernois, J. C. Díaz-Vélez, K. Engel, O. Enríquez-Rivera, D. W. Fiorino, H. Fleischhack, N. Fraija, J. A. García-González, F. Garfias, A. González Muñoz, M. M. González, J. A. Goodman, Z. Hampel-Arias, J. P. Harding, S. Hernandez, A. Hernandez-Almada, F. Hueyotl-Zahuantitla, P. Hüntemeyer, A. Iriarte, A. Jardin-Blicq, V. Joshi, S. Kaufmann, R. J. Lauer, W. H. Lee, D. Lennarz, H. León Vargas, J. T. Linnemann, A. L. Longinotti, G. Luis-Raya, R. Luna-García, R. López-Coto, K. Malone, S. S. Marinelli, O. Martinez, I. Martinez-Castellanos, J. Martínez-Castro, J. A. Matthews, P. Miranda-Romagnoli, E. Moreno, M. Mostafá, L. Nellen, M. Newbold, M. U. Nisa, R. Noriega-Papaqui, R. Pelayo, J. Pretz, E. G. Pérez-Pérez, Z. Ren, C. D. Rho, N. L. Rodd, D. Rosa-González, M. Rosenberg, E. Ruiz-Velasco, B. R. Safdi, H. Salazar, F. Salesa Greus, A. Sandoval, M. Schneider, G. Sinnis, A. J. Smith, R. W. Springer, P. Surajbali, I. Taboada, O. Tibolla, K. Tollefson, I. Torres, T. N. Ukwatta, G. Vianello, L. Villaseñor, T. Weisgarber, S. Westerhoff, I. G. Wisher, J. Wood, T. Yapici, G. B. Yodh, P. W. Younk, A. Zepeda, H. Zhou, and J. D. Álvarez, J. Cosmol. Astropart. Phys. 2018, 049 (2018), arXiv:1710.10288 [astro-ph.HE] .
Albert et al. [2020] A. Albert, R. Alfaro, C. Alvarez, J. C. Arteaga-Velázquez, et al., Phys. Rev. D 101, 103001 (2020).
Tak et al. [2022] D. Tak, M. Baumgart, N. L. Rodd, and E. Pueschel, The Astrophysical Journal Letters 938, L4 (2022).
Griest and Kamionkowski [1990] K. Griest and M. Kamionkowski, Phys. Rev. Lett. 64, 615 (1990).
Smirnov and Beacom [2019] J. Smirnov and J. F. Beacom, Phys. Rev. D 100, 043029 (2019), arXiv:1904.11503 [hep-ph] .
Harigaya et al. [2016] K. Harigaya, M. Ibe, K. Kaneta, W. Nakano, and M. Suzuki, Journal of High Energy Physics 2016, 151 (2016), arXiv:1606.00159 [hep-ph] .
Geller et al. [2018] M. Geller, S. Iwamoto, G. Lee, Y. Shadmi, and O. Telem, JHEP 06, 135, arXiv:1802.07720 [hep-ph] .
Acharyya et al. [2023] A. Acharyya, A. Archer, P. Bangale, J. T. Bartkoske, P. Batista, M. Baumgart, W. Benbow, J. H. Buckley, A. Falcone, Q. Feng, J. P. Finley, G. M. Foote, L. Fortson, A. Furniss, G. Gallagher, W. F. Hanlon, O. Hervet, J. Hoang, J. Holder, T. B. Humensky, W. Jin, P. Kaaret, M. Kertzman, M. Kherlakian, D. Kieda, T. K. Kleiner, N. Korzoun, F. Krennrich, M. J. Lang, M. Lundy, G. Maier, C. E. McGrath, P. Moriarty, S. O’Brien, R. A. Ong, K. Pfrang, M. Pohl, E. Pueschel, J. Quinn, K. Ragan, P. T. Reynolds, E. Roache, N. L. Rodd, J. L. Ryan, I. Sadeh, L. Saha, M. Santander, G. H. Sembroski, R. Shang, M. Splettstoesser, D. Tak, J. V. Tucci, V. V. Vassiliev, and D. A. Williams, Astrophys. J. 945, 101 (2023), arXiv:2302.08784 [astro-ph.HE] .
Albert et al. [2022] A. Albert, S. Alves, M. André, M. Anghinolfi, G. Anton, M. Ardid, S. Ardid, J. J. Aubert, J. Aublin, B. Baret, S. Basa, B. Belhorma, M. Bendahman, F. Benfenati, V. Bertin, S. Biagi, M. Bissinger, J. Boumaaza, M. Bouta, M. C. Bouwhuis, H. Brânzaş, R. Bruijn, J. Brunner, J. Busto, B. Caiffi, D. Calvo, A. Capone, L. Caramete, J. Carr, V. Carretero, S. Celli, M. Chabab, T. N. Chau, R. Cherkaoui El Moursli, T. Chiarusi, M. Circella, A. Coleiro, R. Coniglione, P. Coyle, A. Creusot, A. F. Díaz, G. de Wasseige, C. Distefano, I. di Palma, A. Domi, C. Donzaud, D. Dornic, D. Drouhin, T. Eberl, T. van Eeden, D. van Eijk, N. El Khayati, A. Enzenhöfer, P. Fermani, G. Ferrara, F. Filippini, L. Fusco, Y. Gatelet, P. Gay, H. Glotin, R. Gozzini, R. Gracia Ruiz, K. Graf, C. Guidi, S. Hallmann, H. van Haren, A. J. Heijboer, Y. Hello, J. J. Hernández-Rey, J. Hößl, J. Hofestädt, F. Huang, G. Illuminati, C. W. James, B. Jisse-Jung, M. de Jong, P. de Jong, M. Kadler, O. Kalekin, U. Katz, N. R. Khan-Chowdhury, A. Kouchner, I. Kreykenbohm, V. Kulikovskiy, C. Lagunas Gualda, R. Lahmann, R. Le Breton, S. Lestum, D. Lefèvre, E. Leonora, G. Levi, M. Lincetto, D. Lopez-Coto, S. Loucatos, L. Maderer, J. Manczak, M. Marcelin, A. Margiotta, A. Marinelli, J. A. Martínez-Mora, B. Martino, K. Melis, P. Migliozzi, A. Moussa, R. Muller, L. Nauta, S. Navas, E. Nezri, B. Ó Fearraigh, A. Păun, G. E. Păvălaş, C. Pellegrino, M. Perrin-Terrin, V. Pestel, P. Piattelli, C. Pieterse, C. Poirè, V. Popa, T. Pradier, N. Randazzo, D. Real, S. Reck, G. Riccobene, A. Romanov, F. Sala, A. Sánchez-Losa, F. Salesa Greus, D. F. E. Samtleben, M. Sanguineti, P. Sapienza, J. Schnabel, J. Schumann, F. Schüssler, J. Seneca, M. Spurio, T. Stolarczyk, M. Taiuti, Y. Tayalati, S. J. Tingay, B. Vallage, V. van Elewyck, F. Versari, S. Viola, D. Vivolo, J. Wilms, S. Zavatarelli, A. Zegarelli, J. D. Zornoza, J. Zúñiga, and ANTARES Collaboration, J. Cosmol. Astropart. Phys. 2022, 028 (2022), arXiv:2203.06029 [astro-ph.HE] .
Cirelli et al. [2011] M. Cirelli, G. Corcella, A. Hektor, G. Hütsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, J. Cosmol. Astropart. Phys. 2011, 051 (2011), arXiv:1012.4515 [hep-ph] .
Bauer et al. [2021] C. W. Bauer, N. L. Rodd, and B. R. Webber, Journal of High Energy Physics 2021, 121 (2021), arXiv:2007.15001 [hep-ph] .
Arina et al. [2024] C. Arina, M. Di Mauro, N. Fornengo, J. Heisig, A. Jueid, and R. Ruiz de Austri, J. Cosmol. Astropart. Phys. 2024, 035 (2024), arXiv:2312.01153 [astro-ph.HE] .
Circiello et al. [2024] A. Circiello, A. McDaniel, A. Drlica-Wagner, C. Karwin, M. Ajello, M. Di Mauro, and M. Sánchez-Conde, arXiv e-prints , arXiv:2404.01181 (2024), arXiv:2404.01181 [astro-ph.HE] .
Ando et al. [2020] S. Ando, A. Geringer-Sameth, N. Hiroshima, S. Hoof, R. Trotta, and M. G. Walker, Phys. Rev. D 102, 061302 (2020), arXiv:2002.11956 [astro-ph.CO] .
Note [1] https://github.com/shinichiroando/dwarf_params.
Graus et al. [2019] A. S. Graus, J. S. Bullock, T. Kelley, M. Boylan-Kolchin, S. Garrison-Kimmel, and Y. Qi, Mon. Not. R. Astron. Soc. 488, 4585 (2019), arXiv:1808.03654 [astro-ph.GA] .
Pace and Strigari [2019] A. B. Pace and L. E. Strigari, Mon. Not. R. Astron. Soc. 482, 3480 (2019), arXiv:1802.06811 [astro-ph.GA] .
Strigari et al. [2007] L. E. Strigari, S. M. Koushiappas, J. S. Bullock, and M. Kaplinghat, Phys. Rev. D 75, 083526 (2007), arXiv:astro-ph/0611925 [astro-ph] .
Bullock et al. [2001] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin, J. R. Primack, and A. Dekel, Mon. Not. R. Astron. Soc. 321, 559 (2001), arXiv:astro-ph/9908159 [astro-ph] .
Note [2] We note that Draco II is one of the most dark-matter rich dSphs, but is not included in our deep-exposure sample because, during the data collection phase, Ando+20 had not yet been published.
Song et al. [2024] D. Song, N. Hiroshima, and K. Murase, J. Cosmol. Astropart. Phys. 2024, 087 (2024), arXiv:2401.15606 [astro-ph.HE] .
Holder et al. [2008] J. Holder et al., AIP Conference Proceedings 1085, 657 (2008), https://aip.scitation.org/doi/pdf/10.1063/1.3076760 .
Park et al. [2015] N. Park et al., in 34th International Cosmic Ray Conference (ICRC2015), International Cosmic Ray Conference, Vol. 34 (2015) p. 771, arXiv:1508.07070 [astro-ph.IM] .
Kieda et al. [2013] D. Kieda et al., in 33rd International Cosmic Ray Conference (ICRC2013), International Cosmic Ray Conference, Vol. 33 (2013) p. 1124, arXiv:1308.4849 [astro-ph.IM] .
Adams et al. [2022] C. B. Adams, W. Benbow, A. Brill, J. H. Buckley, J. L. Christiansen, A. Falcone, Q. Feng, J. P. Finley, G. M. Foote, L. Fortson, A. Furniss, C. Giuri, D. Hanna, T. Hassan, O. Hervet, J. Holder, B. Hona, T. B. Humensky, W. Jin, P. Kaaret, T. K. Kleiner, S. Kumar, M. J. Lang, M. Lundy, G. Maier, P. Moriarty, R. Mukherjee, M. Nievas Rosillo, S. O’Brien, N. Park, S. Patel, K. Pfrang, M. Pohl, R. R. Prado, E. Pueschel, J. Quinn, K. Ragan, P. T. Reynolds, D. Ribeiro, E. Roache, J. L. Ryan, M. Santander, A. Weinstein, D. A. Williams, and T. J. Williamson, Astron. Astrophys. 658, A83 (2022), arXiv:2111.04676 [astro-ph.IM] .
Fomin et al. [1994] V. P. Fomin, A. A. Stepanian, R. C. Lamb, D. A. Lewis, M. Punch, and T. C. Weekes, Astroparticle Physics 2, 137 (1994).
Cogan [2008] P. Cogan, in International Cosmic Ray Conference, International Cosmic Ray Conference, Vol. 3 (2008) pp. 1385–1388, arXiv:0709.4233 [astro-ph] .
Maier and Holder [2017] G. Maier and J. Holder, in 35th International Cosmic Ray Conference (ICRC2017), International Cosmic Ray Conference, Vol. 301 (2017) p. 747, arXiv:1708.04048 [astro-ph.IM] .
Krause et al. [2017] M. Krause, E. Pueschel, and G. Maier, Astroparticle Physics 89, 1 (2017).
Christiansen and VERITAS Collaboration [2017] J. Christiansen and VERITAS Collaboration, in 35th International Cosmic Ray Conference (ICRC2017), International Cosmic Ray Conference, Vol. 301 (2017) p. 789, arXiv:1708.05684 [astro-ph.IM] .
Note [3] Specifically, a run’s low-energy threshold is established as the higher of two values: the energy corresponding to 15% of the maximum point in the effective area curve or the energy at which the deviation between the true energy and the reconstructed energy is less than 20%.
Berge et al. [2007] D. Berge, S. Funk, and J. Hinton, Astron. Astrophys. 466, 1219 (2007), arXiv:astro-ph/0610959 [astro-ph] .
Aleksić et al. [2012] J. Aleksić, J. Rico, and M. Martinez, J. Cosmol. Astropart. Phys. 2012, 032 (2012), arXiv:1209.5589 [astro-ph.HE] .
Aleksić et al. [2014] J. Aleksić et al., J. Cosmol. Astropart. Phys. 2014, 008 (2014), arXiv:1312.1535 [hep-ph] .
Li and Ma [1983] T. Li and Y. Ma, Astrophysical Journal 272, 317 (1983).
Bergström and Kaplan [1994] L. Bergström and J. Kaplan, Astroparticle Physics 2, 261 (1994), arXiv:hep-ph/9403239 [hep-ph] .
Lin et al. [2022] Y.-H. Lin, H.-W. Hammer, and U.-G. Meißner, Phys. Rev. Lett. 128, 052002 (2022), arXiv:2109.12961 [hep-ph] .
Martinez et al. [2011] G. D. Martinez, Q. E. Minor, J. Bullock, M. Kaplinghat, J. D. Simon, and M. Geha, Astrophys. J. 738, 55 (2011), arXiv:1008.4585 [astro-ph.GA] .
McGaugh et al. [2021] S. S. McGaugh, F. Lelli, J. M. Schombert, P. Li, T. Visgaitis, K. S. Parker, and M. S. Pawlowski, Astron. J. 162, 202 (2021), arXiv:2109.03251 [astro-ph.GA] .
Geringer-Sameth et al. [2015] A. Geringer-Sameth, S. M. Koushiappas, and M. Walker, Astrophys. J. 801, 74 (2015), arXiv:1408.0002 [astro-ph.CO] .

Appendix A Deep exposure dwarfs

Based on their J-factors, we have observed three targets with deeper exposures ( $>$ 100 hours): Segue 1, Ursa Major II, and Ursa Minor. For these three dSphs, we compute the UL curves on the dark matter annihilation cross section and the null hypothesis bands for two annihilation channels, $\tau^{+}\tau^{-}$ and $b\bar{b}$ . As shown in Fig. 7, Segue 1 exhibits a deviation between the observed UL curve and the null hypothesis band. This deviation arises from Poisson fluctuations in the high-energy regime (above 10 TeV). The dark matter analysis without events beyond 10 TeV eliminates this deviation. As Segue 1 heavily influences the joint-fit UL curve, a similar deviation shows up in the joint-fit result (Fig. 3). In the case of other two dSphs, the observed UL curve is consistent with the null hypothesis band.

Appendix B Comparison Ando+20 profile with GS+15 profile

The J-factors estimated in GS+15 and Ando+20 show noticeable differences. Table 3 summarises the J-factor at 0.5^∘ for dSphs from the two studies. They differ by a factor of 0.5 to 4.5, depending on the dSph. Generally, GS+15 tends to have much higher J-factors for the dark matter-dense dSphs, compared to Ando+20. Conversely, in the case of dSphs with lower dark matter density, GS+15 reports smaller J-factors than those reported by Ando+20. This disparity is particularly decisive in our study because the J-factors of dSphs with deeper exposure (Segue 1, Ursa Major II, and Ursa Minor) from Ando+20 are approximately 2.8 and 4.3 times lower than values from GS+15.

As we lack access to the J-factor parameters from GS+15 for our complete set of 17 dSphs, an accurate computation of the extent to which this discrepancy impacts the upper limit curve for the dark matter cross section cannot be achieved. Nonetheless, assuming the two J-factors vary across all angles by a constant factor for each dSph, a rough estimation of the dark matter upper limit curve can be obtained; i.e., for each dSph, the number of expected events from the dark matter signal (as described in Equation 5) is multiplied by a constant factor based on Table 3. Note that for Draco II and Triangulum II, we assume the ratio is 1 because GS+15 did not present the J-factors for those dSphs. As a result, we find that the upper limit is lowered by a factor of 3, consistent with the J-factor difference between GS+15 and Ando+20 for the dSphs for which VERITAS has collected deep exposures.

Dwarf	log ${}_{10}J(0.5^{o})$ [GeV²/cm⁵]		ratio
	Ando+20	GS+15	GS+15 : Ando+20
Böotes	$17.77_{-0.24}^{+0.23}$	$18.24^{+0.40}_{-0.37}$	2.95
Coma Berenices	$18.37_{-0.33}^{+0.30}$	$19.02^{+0.37}_{-0.41}$	4.47
CVn I	$17.38_{-0.11}^{+0.11}$	$17.44^{+0.37}_{-0.28}$	1.15
CVn II	$17.19_{-0.47}^{+0.37}$	$17.65^{+0.45}_{-0.43}$	2.88
Hercules I	$16.93_{-0.39}^{+0.34}$	$16.86^{+0.74}_{-0.68}$	0.85
Leo I	$17.70_{-0.08}^{+0.07}$	$17.84^{+0.20}_{-0.16}$	1.38
Leo II	$17.54_{-0.10}^{+0.10}$	$17.97^{+0.20}_{-0.18}$	2.69
Leo IV	$16.56_{-0.66}^{+0.57}$	$16.32^{+1.06}_{-1.69}$	0.58
Leo V	$16.58_{-0.69}^{+0.60}$	$16.37^{+0.94}_{-0.87}$	0.62
Segue 1	$18.91_{-0.48}^{+0.39}$	$19.36^{+0.32}_{-0.35}$	2.82
Segue 2	$17.23_{-0.99}^{+0.58}$	$16.21^{+1.06}_{-0.98}$	0.10
Sextans I	$18.05_{-0.29}^{+0.25}$	$17.92^{+0.35}_{-0.29}$	0.74
Ursa Major I	$18.19_{-0.25}^{+0.22}$	$17.87^{+0.56}_{-0.33}$	0.48
Ursa Major II	$18.79_{-0.48}^{+0.36}$	$19.42^{+0.44}_{-0.42}$	4.27
Ursa Minor	$18.47_{-0.22}^{+0.20}$	$18.95^{+0.26}_{-0.18}$	3.02

Table 3: The J-factor comparison between [26] and [52]. Note that the J-factors for Draco II and Triangulum II are not included here due to their absence in [52].