Introduction

Dark matter (DM) is an unrevealed component of the matter in the universe whose existence is widely supported by a broad set of observations (Bertone & Hooper 2018). For decades, many theoretical candidates have been considered for particle DM, of which two representative examples are ultralight axions (M χ = 1 eV) and weakly interacting massive particles (WIMPs;

). Both candidates have been hunted for with state-of-the-art experiments and observatories, and although these searches will continue to achieve important milestones-for example the long sought-after Higgsino may soon be within reach (Rinchiuso et al. 2021;Dessert et al. 2022)-so far the program has been unsuccessful (for the latest reviews, see, e.g., Gaskins 2016;Boveia & Doglioni 2018;Tao 2020).

The longstanding lack of a DM signal detection has driven theorists to look for DM candidates beyond the conventional parameter space. One such candidate is ultraheavy DM (UHDM; 10 TeV  M χ  m pl ≈ 10 19 GeV). Depending on the cosmological scenario and beyond the Standard Model (SM) theory that predicts UHDM, its abundance and properties can vary (for a broad outline, see Carney et al. 2022); e.g., WIMPzilla (Kolb et al. 1999) and Gluequark DM (Contino et al. 2019). In addition to unexplored UHDM candidates, there are models that extend the WIMP mass range beyond ∼10 TeV (e.g., von Harling & Petraki 2014; Baldes & Petraki 2017;Cirelli et al. 2019;Bhatia & Mukhopadhyay 2021). Yet there exists a general upper limit (UL) on the WIMP mass, known as the unitarity limit, which requires M χ  194 TeV (Griest & Kamionkowski 1990;Smirnov & Beacom 2019). This bound arises as the standard WIMP paradigm is associated with a thermal-relic cosmology. In this scenario, in the early universe, the DM and SM particles are in thermal equilibrium. As the universe expands and cools, the DM departs from equilibrium and its abundance is rapidly depleted by annihilations, until the expansion eventually shuts this process off and the relic abundance freezes out. The key parameter in this scenario is the DM annihilation cross section, which for point-like particles going to SM states must scale as M 2 c -by dimensional analysis. As the mass increases, the cross section generally decreases. If it becomes too small, then the DM will be insufficiently depleted by the time it freezes out, and too much DM will remain to be consistent with the observed cosmological density. Ultimately, as unitarity dictates that the cross section cannot be made arbitrarily large, this constraint translates into the stated upper bound on the DM mass.

While there is an attractive simplicity to the thermal-relic cosmology so described, as soon as we allow even minimal departures from it, the unitarity bound can be violated, allowing for the possibility that DM with even higher masses could be annihilating in the present-day universe. For example, instead of annihilating directly to SM states, the DM could produce a metastable dark state which itself decays to the SM. As shown by Berlin et al. (2016), if this dark state lives long enough to dominate the energy density of the universe, its decay to the SM will then dilute the DM density, avoiding the overproduction otherwise associated with heavy thermal DM, and allowing masses up to 100 PeV to be obtained. PeV-scale thermal DM can also be achieved if the DM is a composite state, rather than a point-like particle. Exactly such a scenario was considered by Harigaya et al. (2016), where DM with a large radius arose from a model of a strongly coupled confining theory in the dark sector. The lightest baryon in the theory plays the role of DM, which annihilates through a portal coupling to eventually produce SM states. Such a scenario can evade the unitarity bound as the annihilation cross section is no longer guaranteed to scale as M 2 c -; it can instead now be determined by the geometric size of the composite DM. Indeed, we will see that such composite DM scenarios are broadly the models that can be probed using the observational strategies considered in this work.

The self-annihilations which play a role in setting the DM abundance in the early universe can also be active today, producing an observable flux of stable SM particles such as e ± , ν e,μ,τ , and γ-rays, as well as unstable quarks, leptons, and bosons whose interaction processes can produce secondary γrays. The full energy spectrum at production can be estimated with Monte Carlo (MC) simulations of the underlying particle physics. For this purpose, PYTHIA is the most widely used program, providing an accurate prompt DM spectrum up to ( ) 10 

TeV (Sjöstrand et al. 2008), and is a central ingredient in the widely used PPPC4DMID (Ciafaloni et al. 2011;Cirelli et al. 2011). However, PYTHIA is not appropriate for studying UHDM in general, as it omits many of the interactions in the full, unbroken SM that become important as the UHDM mass becomes much larger than the electroweak scale. An alternative approach was introduced by Bauer et al. (2021), who computed the prompt DM spectrum from 1 TeV up to the Planck scale, the so-called HDMSpectrum. 4 To do so, the authors of that work mapped the calculation of the DM spectrum to the computation of fragmentation functions, which can then be computed with the DGLAP evolution in a manner that includes all relevant SM interactions, providing a better characterization of the prompt UHDM spectrum (see Bauer et al. 2021, for a discussion of earlier approaches to computing DM spectra).

When γ-rays are produced from DM annihilation throughout the universe, they can propagate to the Earth and be detected. After considering the propagation effects,foot_1 the γ-ray flux at the Earth from DM annihilation can be described as

where 〈σv〉 is the velocity-averaged annihilation cross section. The prompt energy spectrum, ( ) dN E dE g , depends on the DM annihilation channel and is determined from the heavy dark matter spectrum, and ( ˆ) ln r is the DM density along the line of sight (LOS). Even though the DM annihilation process can occur anywhere that DM is present, the DM signature from DM-rich regions will be brighter. For instance, dwarf spheroidal galaxies (dSphs) in the Local Group are one of the best targets for DM study because of their high mass-tolight ratio (implying a high DM density; e.g., M/L ∼ 3400 M ☉ /L ☉ for Segue 1; Simon et al. 2011), close proximity, and absence of bright nearby background sources.

The γ-rays that could be arriving at Earth from DM annihilations would be detectable with γ-ray space telescopes and ground-based observatories, enabling indirect searches for DM. The self-annihilation of UHDM can produce γ-rays from around a TeV to above a PeV, containing the energy band in which ground-based γ-ray observatories have better sensitivity than space-based instruments. There are two classes of groundbased very-high-energy (VHE; >100 GeV) γ-ray observatories: Imaging Atmospheric Cherenkov Telescope arrays (IACTs) and extended air shower arrays (EAS). IACTs use reflecting dishes and fast cameras (generally using photomultiplier tubes; PMTs) to reconstruct the Cherenkov light stimulated by air showers triggered by TeV γ-rays as they interact with Earth's atmosphere. Current-generation EAS arrays are made of water tanks, where optical detectors (generally PMTs) in each tank directly detect Cherenkov radition from charged air shower particles. Both types of instrument can reconstruct TeV γ-rays (Funk 2015). Both have been used for indirect DM searches, with a particular focus on searches for electroweak-scale WIMPs (e.g., Aleksić et al. 2014;Archambault et al. 2017;Albert et al. 2018;Abdalla et al. 2018;Abeysekara et al. 2018;Acciari et al. 2022). In addition to those γ-ray observatories, neutrino observatories have also searched for an indirect DM signal (e.g., Aartsen et al. 2017;Albert et al. 2022).

In this paper, we explore the feasibility of detecting a UHDM annihilation signature from dSphs with current and future ground-based VHE γ-ray observatories. To this end, we use only publicly available resources. Also, we compute the expected ULs for an UHDM particle with a mass from 30 TeV to 30 PeV, assuming that the UHDM signal is not detected. We take Segue 1, one of the local classical dSphs, as our benchmark target, because it has been widely used for indirect DM searches, making it possible to place our results in the context of existing limits at lower masses (e.g., Aleksić et al. 2014;Archambault et al. 2017). Furthermore, it has good visibility (in terms of zenith angle of observation) for all of the instruments discussed in this work. We consider three instruments: the Very Energetic Radiation Imaging Telescope Array System (VERITAS; IACT), the Cherenkov Telescope Array (CTA; IACT), and the High-Altitude Water Cherenkov Observatory (HAWC; EAS array). For VERITAS and HAWC, we do not access the official instrument response functions (IRFs)foot_2 and/or observed background spectra, but rather make reasonable assumptions based on publicly available information, and introduce a VERITAS-like and a HAWC-like instrument.

The remaining discussion is organized as follows. In Section 2, we present the theoretical motivations for UHDM searches, with a particular focus on the experimentally accessible parameter space. The data acquisition and processing for each instrument is detailed in Section 3, with the methods used to calculate the projected sensitivity and ULs for each instrument outlined in Section 4. We present our results in Section 5, and the studies on the systematic and statistical uncertainties are discussed in Section 6. Our conclusions are reserved for Section 7.

Theoretical Motivation

Theoretical arguments for DM have often downplayed the UH mass regime. The prejudice against heavier masses arises from the so-called unitarity limit of Griest & Kamionkowski (1990), which is based on the following "bottom-up" argument. The naive expectation is that DM annihilation rates for pointlike particles will scale as v C M 2 s á ñ ~c , where M χ is the particle mass and C is a dimensionless parameter. For a thermal relic, this cross section is what depletes the DM abundance away from its equilibrium value once the temperature of the universe drops below M χ , and so we expect Ω χ ∝ 1/〈σv〉. Accordingly, for too-large M χ , DM cannot destroy itself with enough vigor, and the universe overcloses. One can boost the size of C, but only up to an amount allowed by unitarity. DM as a simple self-annihilating thermal relic is only possible for masses up to ∼194 TeV (Smirnov & Beacom 2019). We show this UL in Figure 1; 194 TeV is an updated value of the conservative bound from Griest & Kamionkowski (1990) (those authors used Ω χ h 2 = 1, as opposed to the current measurement of Ω χ h 2 = 0.12 given by Aghanim et al. 2020).

To derive M χ  194 TeV, one assumes that the annihilation rate saturates the unitarity limit (〈σv〉 ∝ 1/v; see Equation (2) with J = 0) for the entire relevant history of the DM. A rate that scales inversely with velocity is typically found only at low velocities and in the presence of a long-range force, as in the celebrated case of Sommerfeld enhancement. As discussed below, it is difficult to model-build a scenario where the cross section is maximally large, but where the DM continues to behave as a simple elementary particle. Typically, bound-state and compositeness effects will enter in this limit. For such reasons, in Griest & Kamionkowski (1990), the authors felt the above cross-section scaling was overly conservative. Instead, they assumed that the cross section was dominantly S-wave (〈σv〉 ∝ v 0 ) but with a maximum value still set by unitarity (as given in Equation ( 2)). Using this, and assuming Ω χ h 2 = 1, they derived the well-known UL of 340 TeV. Repeating their calculation for Ω χ h 2 = 0.12, the bound is reduced to M χ  116 TeV. Nevertheless, we will adopt the more conservative value of 194 TeV in our results. It involves the fewest assumptions about the early universe, but amounts to assuming that DM finds a way to annihilate at the limiting cross-section value throughout the era that set its relic abundance.

The presence of additional structure in either the DM particles themselves or the final states they capture into can weaken even this conservative limit, though. For example, if capture into bound states is possible, then selection rules can open up annihilation channels into higher partial waves. The total relic abundance of DM is necessarily set by the sum over all channels, but each partial wave respects the limit from unitarity unitarity:

As discussed by Bottaro et al. (2022), even for the straightforward scenario of thermal relics that are just multiplets of the electroweak group SU(2) L , this allows DM consistent with unitarity up to ∼325 TeV. It would seem uncontroversial to analyse the full regime that allows this simple scenario.

To relax the bound farther, as mentioned above, the unitarity limit of roughly 100 TeV assumes a point-like particle. This was explicitly recognized in the classic 1990 reference on the matter. If, however, DM is a composite particle, then the relevant dimensionful scale that sets the annihilation rate can be its geometric size, R, which may be much larger than its Compton wavelength, ∼1/M χ . It is thus possible to realize a thermal-relic scenario for masses ?100 TeV (e.g., the example of Harigaya et al. (2016) discussed above). 7 For pointing telescopes like VERITAS, H.E.S.S., or CTA to have a discovery advantage, one needs a scenario, like compositeness, with non-negligible DM annihilation, since the resulting flux will scale like ρ 2 . Bound-state particles with a heavy constituent, whether obtained as thermal relics or by a more complicated cosmology, provide a means to get annihilation rates of

, where C unitary is the largest factor consistent with quantum mechanics in a single partial wave. One may therefore consider this as a generalization of the "sum over partial waves" loophole we first mentioned in the boundstate capture scenario. As we see in Figure 1, there is a large region of parameter space beyond the point-like unitarity limit. Furthermore, we project that the limits from CTA exceed those from HAWC out to several PeV, and are primed for testing these models.

The generic possibility of a geometric cross section for composite particles can be seen with atomic (anti)hydrogen, as pointed out by Geller et al. (2018), whose arguments we briefly recap. In a hydrogen-antihydrogen collision, an interaction with a geometric cross section is the "rearrangement" reaction, which produces a protonium ( pp) + positronium (e + e -) final Figure 1. A comparison of our estimated limits for annihilation to t t against various theoretical benchmarks. The black solid curve refers to the standard thermal-relic cross section (2.4 × 10 -26 cm 3 s -1 ; Steigman et al. 2012), and the region shaded in gray is the conventional parameter space associated with a point-like thermal relic. For Segue 1, the J = 0 partial-wave unitarity limit on a point-like annihilation cross section is shown in orange-irrespective of the early universe cosmology, point-like particles can only annihilate at a rate below this. Composite states are not so restrictive, however, and can annihilate up to the various composite unitarity bounds. For a detailed discussion, see Section 2. state. Partial wave by partial wave, unitarity is naturally respected. However, summing over all allowed angular momenta gives

where k i is the initial momentum, R is the size of the particle, and J max is set by angular momentum conservation and the classical value (k i R). 8 Importantly, a parametric enhancement in the cross section has been achieved by saturating each partial-wave bound up to J max . Whatever partial-wave protonium is captured into, it will ultimately decay down the spectroscopic ladder until reaching the lowest allowed energy state, at which point it annihilates. For a generic scenario with the dark sector charged under the SM, the entire process of capture, decay, and annihilation is prompt on observational timescales. An UH dark-hydrogen thus provides a proof of concept for a "detection-throughannihilation" scenario. The argument for geometric scaling generalizes, though, to include states bound by strong dynamics (Jacoby & Nussinov 2007;Kang et al. 2008). Thus, DM may be more like an UH B-meson (as studied by Geller et al. 2018), or a gluequark (adjoint fermion with color neutralized by a cloud of dark gluons; Contino et al. 2019), heavy-light baryon (Harigaya et al. 2016), etc. For a complete scenario, one would necessarily need an explanation for why these heavy-constituent composites came to be the DM with the right abundance. Nonetheless, the physics behind their ability to annihilate with an effective rate far above the pointparticle unitarity limit is straightforward. Therefore, models with dynamics not too different from the SM can realize annihilating particle DM all the way to the Planck scale, and should be tested.

With the above in mind, in Figure 1, we outline basic theoretical aspects of the parameter space we will consider (see Albert et al. 2022). First, we see that the majority of the mass range probed is above the conventional unitarity limit. Next, the curve we label as "Partial-Wave Unitarity" represents the largest present-day annihilation cross section consistent with the same point-particle unitarity constraints that, when applied in the early universe, constrains M χ  194 TeV. In particular, we require

, where we take v rel ∼ 2 × 10 -5 as an approximate value for the average velocity between DM particles in nearby dwarf galaxies (Martinez et al. 2011;McGaugh et al. 2021). 9 Composite states can readily evade this bound, although as shown by Griest & Kamionkowski (1990), even these systems eventually hit a "composite unitarity" bound, which requires

, and which for large masses reduces to the result in Equation (3). We show this result for different values of R in Figure 1, and note that for M χ = R -1 , these results reduce to the point-like unitarity limit.

Data Reduction

VERITAS-like Instrument

VERITAS is an array of four imaging atmospheric Cherenkov telescopes located in Arizona, USA (Weekes et al. 2002). One of the VERITAS scientific programs is to search for indirect DM signals from astrophysical objects such as dSphs and the Milky Way Galactic Center (Zitzer et al. 2017). Since it has a similar sensitivity to other IACT observatories like MAGIC and H.E.S.S. (Aharonian et al. 2006;Park et al. 2015;Aleksić et al. 2016), we adopt VERITAS as representative of current-generation IACTs.

For our analysis, we take the published IRFs and observed ON and OFF regionfoot_5 counts from Archambault et al. (2017). The size of the ON region was 0.03 deg 2 , and the OFF region was defined by the crescent background method (Zitzer et al. 2013). The relative exposure time between the ON and OFF regions (α) was 0.131. From 92.0 hr of Segue 1 observations, the number of observed events from the ON (N on ) and OFF regions (N off ) was 15895 and 120,826, respectively. We introduce a reference instrument, denoted "VERITAS-like," whose observables are limited to total N on , total N off , and α (see the Appendix for a comparison between the VERITAS and VERITAS-like constraints on the DM annihilation cross section). In addition, we scale down the N on and N off values to a nominal observation time of 50 hr.

CTA

CTA is a next-generation ground-based IACT array, which is expected to have about 10 times better point-source sensitivity when compared with current IACT observatories, in addition to a broader sensitive energy range, stretching from 20 GeV to 300 TeV, and two to five times better energy and angular resolutions (Bernlohr et al. 2013). The observatory will be made up of two arrays, providing full-sky coverage: one in the Northern Hemisphere (CTA-North; La Palma in Spain) and the other in the Southern Hemisphere (CTA-South; Atacama Desert in Chile). CTA will be equipped with tens of telescopes. In this study, we consider the CTA-North array, from which our target, Segue 1, can be observed. CTA will broaden our understanding of the extreme universe, including the nature of DM (Consortium et al. 2019), and will be able to probe longpredicted, but so far untested candidates like Higgsino DM (Rinchiuso et al. 2021).

The CTA IRFs and background distributions as a function of energy, as well as official analysis tools,foot_6 are publicly available (Deil et al. 2017; Cherenkov Telescope Array Observatory & Cherenkov Telescope Array 2021). We assume the alpha configuration (prod5 v0.1). In the alpha configuration, the CTA-North array consists of four Large-Sized Telescopes (LSTs) and nine Medium-Sized Telescopes (MSTs). 12 To compare with the VERITAS-like instrument, we use the same observation conditions; the size of the ON region is set to 0.03 deg 2 with an α of 0.131. = c , comparable to or larger than the incoming particle's binding energy, E b . If E i = E b , then only the Swave will contribute, and the cross section becomes σ ∼ R/k i . Since this involves just a single partial wave, we therefore cannot use a sum with many terms to exceed the point-particle unitarity limit. 9 We note that the location of the partial-wave unitarity bound strongly depends on the system observed. A search for DM annihilation within the Milky Way, for instance, would depend on a higher relative velocity, v rel ∼ 10 -3 , given the larger mass of our galaxy as compared to its satellites. This would lower the "Partial-Wave Unitarity" curve shown in Figure 1 by roughly two orders of magnitude.

HAWC-like Instrument

HAWC, located at Sierra Negra, Mexico, is a γ-ray and cosmic-ray observatory. The instrument constitutes 300 water tanks. Each tank contains about 1.9 × 10 5 L of water with four PMTs. After applying γ/hadron separation cuts, observed γray events are divided into analysis bins ( hit  ) based on the fraction of the number of PMT hits. HAWC observes twothirds of the sky on a daily basis and has found many previously undetected VHE sources (Albert et al. 2020). In addition, they have studied 15 dSphs within its field of view to search for DM annihilation and decay signatures (Albert et al. 2018;Abeysekara et al. 2018).

The IRFs and observed background spectrum for Segue 1 are not publicly available, so we introduce a "HAWC-like" reference instrument based on reasonable assumptions. A data set including 507 days of observations of the Crab Nebula is publicly available (Abeysekara et al. 2017), 13 and the declination angle (decl.) of the Crab Nebula is not significantly different from that of Segue 1 (Δdecl. ≈ 6°). Since the decl. is expected to be one of the key factors determining the shape of the IRFs and background rate, we assume that the background rate and IRFs should be similar for observations of Segue 1 and the Crab Nebula (see the Appendix for a comparison between the HAWC and HAWC-like constraints on the DM annihilation cross section). With the help of the Multi-Mission Maximum Likelihood framework (3ML; Vianello et al. 2015), we acquire the IRFs and background rate for each hit  (a total of nine bins) as used by Abeysekara et al. (2017). We set the radius of an ON region to 0°. 2, and the background is calculated from a circular region with a 3°radius around the Crab Nebula, providing an α of 0.04/9 (∼0.004).

Analysis Methods

Ingredients for Estimating the UHDM Signal

To compute the γ-ray annihilation flux at the Earth, as given in Equation (1), we need two ingredients: the photon spectrum for each DM annihilation channel and the DM density profile of the selected target, Segue 1. As stated, we use the HDMSpectrum (Bauer et al. 2021) to calculate the expected DM signal because it provides an accurate spectrum for the full mass range we consider. The annihilation of UHDM produces γ-rays of energies equal to or less than M χ . We compute the fraction of the produced energy flux (F

) that is observable and the number of expected γ-ray events (N dE dN dE ò µ ); i.e., the energy flux and γ-ray counts distributions within the energy band of the current and future VHE γray observatories (E 100 TeV). In this work, we consider nine annihilation channels: three charged leptons (e + e -, μ + μ -, and τ + τ -), two heavy quarks ( t t and bb), three gauge bosons (W + W -, ZZ, and γγ), and one neutrino ( ē e n n ). For the DM density profile, we take a generalized version of the Navarro-Frenk-White (NFW) profile, which is a function of five parameters (Hernquist 1990;Zhao 1996;Geringer-Sameth et al. 2015)

where the choice of (α, β, γ) = (1, 3, 1) recovers the original NFW profile (Navarro et al. 1997) and r s is the scale radius of the DM halo. The so-called J-factor is defined as the integral of the squared DM density along the LOS within a region of interest (ROI)

The set of five NFW parameters (α, β, γ, ρ s , and r s ) is obtained by fitting the observed kinematic data of the dSphs. Limited data produce large uncertainties in estimates of the J-factor, which propagate as a systematic uncertainty when estimating the DM cross section (see Section 6). In a thorough study, Geringer-Sameth et al. (2015) obtained a number of parameter sets that adequately describes the data. Among more than 6000 sets for Segue 1, we take one that approximates the median of the J-factor (see Table 1).

Figure 2 shows the expected number of γ-ray photons under the conditions stated below (left panel) and the ratio of observable energy flux to the total energy flux (right panel) for the nine annihilation channels. For the expected counts distribution, we assume that the effective area is 10 10 cm 2 , the exposure time is 50 hr, the J-factor is 10 18 GeV 2 cm -5 sr -1 , and the DM cross section is 10 -23 cm 3 s -1 . This result implies that the current and future observatories, whose sensitive energy ranges extend to 100 TeV, can observe a large portion of the produced γ-rays and/or energy flux from UHDM annihilation, up to M χ of a few PeV. For the γγ channel, the majority of the energy remains in the sharp spectral feature at E γ ∼ M χ , and so the energy flux ratio sharply drops once the mass is above 100 TeV and the continuum component becomes dominant. This sharp decrease is not clearly visible in the expected count level because the emission at E γ ∼ M χ produces only about 10% of the total counts in the high-mass regime.

Projected Sensitivity Curves

To explore the feasibility of detection, we compare expected γ-ray counts from UHDM self-annihilation to background counts. The number of expected signal counts (N s ) is obtained by forward-folding Equation (1) with the IRFs

where unprimed and primed quantities represent observed (strictly speaking, reconstructed) and true quantities, respectively. The function

, W ¢ W¢ , refers to an IRF consisting of three sub-functions: effective area, energy bias, and pointspread function. Assuming that the number of ON region events is N on = N s + αN off , we calculate the significance of the UHDM signal by using the so-called Li & Ma significance ( ;  Note. The maximum angular distance, max q , is given by the location of the furthest member star, which is an estimate of the size of Segue 1.

Li & Ma 1983)

Finally, for each annihilation channel, we find a set of values of M χ and 〈σv〉 for which  = 5σ.

Expected UL Curves

To estimate the UL on the UHDM annihilation cross section for a given M χ , we perform a maximum likelihood estimation (MLE). Since we cannot access the energy distribution of background events for the VERITAS-like instrument, we use a simple likelihood analysis using the total N on and N off counts,

where the nuisance parameter b represents the expected background rate. This likelihood function is expected to be less sensitive compared to a full likelihood function incorporating event-wise energy information, especially at high masses, as it does not utilize any features present in the DM spectrum; see Aleksić et al. (2012) for a full discussion of this hindrance.

For CTA and the HAWC-like instrument, we perform a binned likelihood analysis

We calculate the expected UL with the assumption that the ON region does not contain any signal from UHDM selfannihilation but only Poisson fluctuations around α × N off ; i.e., we can randomly sample N on from the Poisson distribution of αN off . For the binned likelihood analysis, we can apply the Poisson fluctuations to each background bin to get the binned ON-region data. With the synthesized ON-region data, we perform an MLE analysis and calculate the UL on the DM cross section for a given M χ . Throughout this paper, UL refers to the one-sided 95% confidence interval, which is obtained from the profile likelihood ( ln 1.35

). We repeat the process of calculating the expected limit to get the median or the containment band for the 95% UL.

Results

Here, we present two sets of analysis results: sensitivity curves and expected ULs, as functions of the UHDM particle mass. Since above a few tens of petaelectronvolts the energy flux ratio for all annihilation channels is less than 10% (Figure 2), we perform the analyses for UHDM masses from 30 TeV up to 30 PeV. Note that all of the following results are based on assumed exposure times of 50 hr for the VERITASlike instrument and CTA-North, and 507 days for the HAWClike instrument.

Figure 3 shows the sensitivity curves for the nine UHDM annihilation channels (e + e -, μ + μ -, t t , bb, W + W -, ZZ, γγ,foot_9 and ē e n n ) for the VERITAS-like (50 hr; left panel), CTA-North (50 hr; middle panel), and HAWC-like (507 days; right panel) instruments. Considering the annihilation of an UHDM particle with an M χ of 1 PeV via the τ + τ -channel, the HAWC-like instrument is likely to reach an  of 5σ with the smallest cross section; specifically, the VERITAS-like instrument is expected to detect UHDM for a cross section of ∼ 5 × 10 -19 cm 3 s -1 , CTA-North for ∼ 4 × 10 -19 cm 3 s -1 , and the HAWC-like instrument for ∼ 1 × 10 -19 cm 3 s -1 . However, this sensitivity depends on the annihilation channel and the UHDM mass, not to mention the exposure time. For example, for an M χ of 100 TeV, CTA-North shows, in general, better sensitivity compared Figure 2. The number of expected γ-ray events (left) and relative ratio between the observable and total γ-ray energy flux (right). The expected counts are computed assuming an effective area of 10 10 cm 2 , 50 hr of exposure time, a J-factor of 10 18 GeV 2 cm -5 sr -1 , and 〈σv〉 = 10 -23 cm 3 s -1 . The observable energy flux is defined as the integrated γ-ray energy flux up to 100 TeV, and for reference in the black dashed curve we show a value of 10%. The portion of the observable UHDM signal from M χ > 100 TeV decreases progressively as M χ increases. The various line styles refer to the classes of annihilation channel: charged leptons (solid), quarks (dashed), gauge bosons (dotted), and ē e n n (dashed-dotted).

to the other instruments. For the γγ channel, a discontinuity in the sensitivity lines can be seen because, as explained earlier, the line-like contribution (E γ ∼ M χ ) falls outside the sensitive energy range. It is worth noting when comparing limits from CTA-North and the VERITAS-like instrument that while the effective area of CTA-North is 4-5 times larger than that assumed for the VERITAS-like instrument, the size of the signal regions differ between the two analyses. The analysis for the VERITAS-like instrument considers a larger ON region. Next, we estimate the ULs on the UHDM annihilation cross section as a function of UHDM particle mass for the same annihilation channels for the three instruments (Figure 4). The curves represent the median value from 100 realizations generated at each mass. With the assumed observation conditions (e.g., livetime), CTA-North shows the most constraining ULs at lower masses (M χ < 1 PeV), whereas the HAWC-like instrument provides more stringent ULs at higher masses. Note that the UL on the DM cross section is expected to decrease as we increase the exposure time, v t 1 UL s á ñ µ . As expected from the relative sensitivity between VERITAS and CTA-North, the UL curves from CTA-North are about 10 times lower than those from the VERITAS-like instrument.

In the case of the γγ annihilation channel, a discontinuity in the UL curve is again observed at 100 TeV, most strongly for CTA-North. In contrast to the VERITAS-like instrument, it is possible for the CTA-North instrument to perform the full binned likelihood analysis by comparing the signal and background energy distributions, which lowers the UL curve (see the Appendix). Note that in the case of the γγ annihilation channel, the two distributions differ clearly compared to those of the other channels. In the case of the HAWC-like instrument, the energy dispersion matrix for the highest energy bin is relatively broad, which smooths out the discontinuity.

Discussion of the Statistical and Systematic Uncertainties

Here we briefly discuss the impact of the statistical and systematic uncertainties on the presented UL curves. For these studies, we consider a single annihilation channel ( t t ) for simplicity, although the results are representative of what we expect for the additional channels.

Due to the Poisson fluctuations in the observed counts, statistical uncertainties are inevitable. For this study, we compute the 68% containment band of the expected UL curves for a large number of MC realizations (10,000), using the method described in Section 4.3. Figure 5 shows the statistical uncertainty band for 68% (shaded region) and 95% (dashed lines) containment. This figure implies that the Poisson fluctuations can result 45%-55% of the statistical uncertainty (at the 1σ level) across all masses for the three instruments: VERITAS-like (∼45%), CTA-North (∼53%), and HAWClike (∼54%).

A major systematic uncertainty, beyond that inherent in IRFs, is the present uncertainty in the DM density profile assumed for Segue 1. A DM density profile estimated from insufficient and possibly inaccurate kinematic observations will inevitably have a large uncertainty. Also, it depends on assumptions and approximations made in the modeling-for instance, the assumption of a NFW profile with exact spherical symmetry-can also lead to systematic uncertainties. A Bayesian approach to find accurate proper profile parameters has been helpful for reducing these effects by considering the physical information of an object (e.g., Ando et al. 2020). In addition to the aforementioned uncertainties, unclear memberships of stars can also enlarge this uncertainty. The stellar sample selection when fitting the DM density profile affects the J-factor significantly, such that any ambiguity in the sample selection, possibly due to contamination from foreground stars or stellar streams, can overestimate the J-factor. The magnitudes of the systematic uncertainties are different from dSph to dSph, and depend on the definition of the DM density profile. For further discussion on this uncertainty, see Bonnivard et al. (2015aBonnivard et al. ( , 2015b) ) As mentioned earlier, Geringer-Sameth et al. (2015) provide more than 6000 viable parameter sets for Segue 1, and we compute 10 4 expected UL curves by randomly sampling the parameter set. In this work, we use the parameter sets to estimate the systematic uncertainty of the expected UL curve due to uncertainty in the J-profile. Note that in this study, we do not include the Poisson fluctuations of the simulated ON region counts; i.e., N on,i is equal to αN off,i . Finally, we take ULs corresponding to the 68% and 95% containment for each mass (Figure 5). This figure implies that, for Segue 1, the J-factor can increase or decrease the UL curve by a factor of 2 (1σ level) across all masses, regardless of the instrumental properties, at the level of the statistical uncertainties seen in Figure 5. Note that Bonnivard et al. (2016) claimed that the J-factor may be overestimated by about two orders of magnitude due to the stellar sample selection bias. However, the accurate prediction of the Segue 1 J-profile is beyond the scope of this paper.

Summary and Outlook

In this work, we have explored the potential of current and future γ-ray observatories to extend the search for DM beyond the unitarity bound. Our results allow one to determine whether the discovery of an UHDM candidate of a given mass and annihilation cross section is within reach. Furthermore, we provide an estimate of the constraints that can be derived on the UHDM annihilation cross section by current and future γ-ray observatories, assuming a non-detection.

Returning to Figure 1, we can place our obtained limits in the context of theoretical constraints on the allowed annihilation cross section of UHDM. All instruments considered can probe realistic cross sections for composite UHDM particles whose annihilation respects partial-wave unitary. For the given exposure times (50 hr for CTA-North and the VERITAS-like instrument, and 507 days for the HAWC-like instrument), CTA-North is projected to provide the most constraining limits, probing scales down to R = (10 GeV) -1 for UHDM with a mass around 300 TeV. At higher masses, above 1 PeV, HAWC-like limits become the most constraining, reaching scales around R = (1 GeV) -1 at 10 PeV. The VERITAS-like limits, while less constraining, are worse than those of CTA-North or the HAWC-like instrument by less than or equal to an order of magnitude for the entire mass range (with a slight advantage over the HAWC-like instrument at masses below 100 TeV). For other current IACTs like MAGIC and H.E.S.S., we expect similar results to that from the VERITAS-like instrument since they have similar sensitivities. In the case of the Large High Altitude Air Shower Observatory (LHAASO; Cao et al. 2019), as we do not have full access to data and IRFs, we could not perform the same analysis. However, we expect roughly 2-10 times better results than those of HAWC, considering the simulation study by He et al. (2019) in the mass range of 1-100 TeV and their relative effective areas. However, such results can vary depending on the annihilation channel considered, background rejection, details of the IRFs such as the point-spread function, and the observing conditions. This work draws attention to the exploration of DM beyond the conventional parameter range. The results we have derived are indicative, using reasonable assumptions about the data and IRFs for current-generation instruments, as well as realistic exposure times for current and future instruments. We hope that this work illustrates the interest and feasibility of searches for UHDM with current-generation γ-ray instruments, and the value of considering such searches for future observatories such as CTA and the Southern Wide-field Gamma-ray Observatory (SWGO; a proposed next-generation EAS observatory; Huentemeyer et al. 2019). The phase space that can be probed, in terms of DM particle mass and annihilation cross section, is a relevant one for models predicting composite UHDM. This parameter space is currently unconstrained, but Figure 5. Left: statistical uncertainties on the expected 95% limits. Each uncertainty band is obtained from 10 4 realizations for the t t annihilation channel. Right: systematic uncertainty on the same expected limits, resulting from uncertainties in the J-factor estimation. Each uncertainty band is obtained from 10 4 realizations for the t t annihilation channel. In both figures, the shaded region refers to 68% containment, and dashed lines are 95% containment.

could be probed with archival data sets from current-generation γ-ray instruments, including HAWC, VERITAS, and other IACTs.

Our work has benefited from discussions with Michael Geller, Diego Redigolo, and Juri Smirnov. We would like to thank Alex Geringer-Sameth, Savvas M. Koushiappas, and Matthew Walker for providing the parameter sets for the J-