# Sensitivity to sub-GeV dark matter from cosmic-ray scattering with very-
high-energy gamma-ray observatories

Igor Reis  Emmanuel Moulin  Aion Viana
Victor P. Goncalves

###### Abstract

Huge efforts have been deployed to detect dark matter (DM) in the GeV-TeV mass
range involving various detection techniques, and led to strong constraints in
the available parameter space. We compute here the sensitivity to sub-GeV DM
that can be probed from the inevitable cosmic-ray scattering onto DM particles
populating the Milky Way halo. Inelastic scattering of energetic cosmic rays
off DM would produce high-energy gamma rays in the final state, providing a
new avenue to probe the poorly-constrained so far sub-GeV dark matter mass
range. In this work we derive sensitivity forecasts for the inelastic cosmic-
ray proton - DM cross section for current and future very-high-energy gamma-
ray observatories such as H.E.S.S., LHAASO, CTA and SWGO in the 100 eV to 100
MeV mass range. These inelastic cross section constraints are converted to the
elastic proton - DM cross section to highlight further complementarity with
cosmological, collider and direct detection searches. The sensitivity computed
at 95% confidence level on the elastic cross section reaches $\sim$2$\times$
10-32 cm2 for a 100 keV DM mass for H.E.S.S.-like and $\sim$7$\times$ 10-34
cm2 for a $\sim$1 keV DM mass for LHAASO. The sensitivity prospects for CTA
and a strawman SWGO model reach $\sim$6$\times$ 10-34 cm2 and $\sim$4$\times$
10-35 cm2, for DM masses of 10 keV and 1 keV, respectively. The sensitivity
reach of the gamma-ray observatories considered here enables to probe an
uncharted region of the DM mass - cross section parameter space.

##  1 Introduction

Dark Matter (DM) pervades the universe at all scales revealing its presence
through gravitational interaction. A substantial body of evidence, such as
galaxy rotation curves, galaxy cluster collisions, the Cosmic Microwave
Background (CMB) and large-scale structure formation, implies that most of the
matter in the universe does not behave as the known baryonic matter [1, 2, 3,
4]. Various approaches have been deployed to probe the non-gravitational
interaction of DM such as its production at colliders [5], its scattering off
target nucleons or electrons in underground detectors [6], commonly referred
to as direct DM detection, or its decays or annihilations probed by indirect
detection [7]. Strong constraints have been derived in the GeV-TeV mass regime
from direct and indirect detection experiments. Indirect searches in gamma
rays set strong constraints on the velocity-weighted annihilation cross
section probing the natural scale expected for thermal-relic Weakly
interacting massive particles (WIMPs) [8, 9]. For DM particle $\chi$ with
masses $m_{\chi}$ above a few ten GeV, current constraints probe DM-nucleon
cross section expected for electroweak DM models [10]. However, the
constraints are dramatically relaxed, by orders of magnitude, for DM particle
with sub-GeV masses. For DM particle velocity of $\sim$10${}^{-3}c$ as
expected in the Milky Way (MW) halo, sub-GeV DM would produce nucleon recoils
at keV energies where standard direct detection experiments loose sensitivity.
Measurements of electronic recoils would be able to probe DM masses down to
about 1 MeV, see, for instance, Ref. [11]. To explore these sub-GeV masses,
current experimental endeavors have shifted focus towards alternative targets
and instrumentation techniques [12]. Lately, an alternative method for
detecting light DM has been explored in the context of direct detection
experiments, which involves observing DM particles that have obtained
relativistic energies via the elastic collision with cosmic rays (CRs),
constituting what has been called cosmic-ray boosted dark matter, see, for
instance, Ref. [13]. This implies standard direct detection experiments to be
able to detect light DM particles, due to the induced much larger recoil
energies. Assuming fermionic DM, kinematic measurements in dwarf spheroidal
galaxies set model-independent lower bounds on DM masses down to $\sim$0.1 keV
depending on the production scenario [14, 15]. In what follows, no specific DM
production scenario is assumed. Should one consider thermally-produced sub-GeV
DM, CMB and BBN constraints would be unavoidable (see, for instance, Refs.
[16, 17, 18]).

Inelastic collisions between CR protons and DM particles in the MW will
inevitably produce gamma rays. The basic idea, represented in Fig. 1, is that
for CRs of energy $E_{\rm p}$ colliding with DM at large center-of-mass
energies $\sqrt{s}$, the protons can be excited or break-up, producing gamma
rays of energy $E_{\gamma}$ in the final state. Over the last decades, several
models have been proposed to describe the proton - DM interactions,
considering different assumptions for the couplings and properties of the
mediator (see, e.g., Refs. [19, 20, 21, 22]). Alternatively one can derive
model-independent upper bounds on the inelastic cross section $\sigma^{\rm
inel}_{\rm p\chi}(\sqrt{s})$, which should be satisfied by all models that
intend to describe the inelastic proton - DM interaction. As the energy
dependence of the inelastic $p\chi$ cross section is still an open question,
we will investigate how such a dependence impacts our findings assuming that
$\sigma^{\rm inel}_{\rm p\chi}(\sqrt{s})$ depends on the energy in a similar
way to that observed in proton-proton, meson-proton and real photon-proton
collisions at high energies [23]. Energetic cosmic-ray protons in the MW would
consequently produce diffuse very-high-energy (VHE, E$\gtrsim$100 GeV) gamma
rays that could be eventually detected by current and future VHE gamma-ray
observatories, such as H.E.S.S. [24], LHAASO [25], MAGIC [26], VERITAS
[27],CTA [28, 29] and SWGO [30], and further used to place limits on the
inelastic cross section. In this context, observations of the Galactic Center
(GC) and the Inner Galaxy are of prime interest to search for such signatures,
as both the CR flux and DM density are measured to be higher there than
anywhere in the MW.

In recent years, the elastic proton-DM interaction strength has been
constrained severely by various direct detection experiments, as well as by
cosmological bounds. In this paper, inspired by Ref. [31], we will investigate
if the bounds derived for $\sigma^{\rm inel}_{\rm p\chi}$ can impose
additional limits on the elastic $p\chi$ cross section ($\sigma^{\rm el}_{\rm
p\chi}$). We will explore such a possibility considering distinct approaches
to relate these quantities and the constraints will be recast on the elastic
cross section on nucleon to highlight further complementary with cosmological,
collider and direct detection searches.

The paper is organized as follows. Section 2 describes the expected gamma-ray
signals from inelastic scattering of high-energy cosmic rays off dark matter
particles populating the MW halo together with gamma-ray backgrounds. In
section 3 we present current and future VHE gamma-ray observatories and their
capabilities used to compute sub-GeV DM sensitivity. Section 4 outlines the
statistical analysis framework used to compute the sensitivity to proton-DM
cross section. Results and discussion are presented in Section 5 in light of
current constraints on the elastic cross section between dark matter and
proton in the sub-GeV DM mass range. Finally, in Section 6 we summarize our
main results and conclusions.

Figure 1: Production of a gamma ray with energy $E_{\gamma}$ in inelastic CR
proton - dark matter ($p\chi$) interactions at a center - of - mass energy
$\sqrt{s}$. In the DM rest frame, one has that
$\sqrt{s}=(2E_{p}m_{\chi}+m_{p}^{2}+m_{\chi}^{2})^{\frac{1}{2}}$, where
$E_{p}$ is the incident CR proton energy and $m_{\rm p}$ and $m_{\chi}$ are
the proton and DM particle mass, respectively. The interaction is described by
the inelastic $p\chi$ cross section, $\sigma^{\rm inel}_{\rm
p\chi}(\sqrt{s})$, which is energy dependent in the general case.

##  2 Gamma-ray production from cosmic-ray scattering off dark matter

###  2.1 Expected VHE gamma-ray signal

Gamma rays of energy $E_{\gamma}$ can be produced in inelastic collisions
between CR protons and DM particles of mass $m_{\rm\chi}$ populating the MW DM
halo. The associated gamma-ray flux is a function of the cross section for the
inelastic interaction between CR protons and DM, $\sigma^{\rm inel}_{\rm
p\chi}$, and can be calculated given the functional form of the cosmic ray
flux, $d\Phi_{\rm p}/dE_{\rm p}$, and the DM distribution $\rho_{\rm DM}(r)$
in the MW. In particular, in the $\delta$-function approximation [32] and for
a power-law proton flux, the expected energy-differential flux in gamma rays
in a solid and angle $\Delta\Omega$ writes as

$\displaystyle\frac{d\Phi_{\gamma}(E_{\gamma})}{dE_{\gamma}}=\frac{4(f_{\pi})^{\Gamma-1}}{\Gamma^{2}}\times\frac{\sigma_{\rm
p\chi}^{\rm
inel}(E_{\gamma}/f_{\pi},m_{\chi})}{m_{\chi}}\times\Phi_{0}(E_{\gamma}/{\rm
10TeV})^{-\Gamma}\times D(\Delta\Omega)\,,$ (2.1)

where $\Gamma$ and $\Phi_{0}$ are the spectral index and normalisation of the
power-law CR proton flux, respectively, and the factor $f_{\pi}$ characterizes
the fraction of the proton energy transferred to the secondary neutral pions
that decay into gamma rays. As in Ref. [32], we will assume in the following
$f_{\pi}=0.17$. The gamma-ray production for the interaction between a CR
proton and a DM particle of mass $m_{\chi}$ requires a minimum CR proton
energy to produce a pion given by $E^{\rm min}_{\rm p}=(m_{\rm
p}+m_{\rm\pi^{0}})\left[1+m_{\rm\pi^{0}}(2m_{\rm
p}+m_{\rm\pi^{0}})/(2m_{\rm\chi}(m_{\rm p}+m_{\rm\pi^{0}}))\right]$, that,
therefore, implies a minimum value for the gamma-ray energy of
$E_{\gamma}^{\rm min}\simeq E_{\rm p}^{\rm min}/f_{\rm\pi}$ at VHE [33, 34].
Finally, the $D$-factor can be defined as the integral of DM density in the MW
halo along the line-of-sight (l.o.s.) $s$ and over the solid angle
$\Delta\Omega$ as

$\displaystyle
D(\Delta\Omega)=\int_{\Delta\Omega}d\Omega\int_{l.o.s.}\rho_{\rm
DM}[r(s)]ds\,.$ (2.2)

The radial coordinate $r$ from the center of the halo is expressed as
$r=\big{(}s^{2}+r_{\odot}^{2}-2\,r_{\odot}\,s\,\cos\theta\big{)}^{1/2}$, where
$r_{\odot}$ is the distance of the observer to the GC, taken to be $r_{\odot}$
= 8.3 kpc, and $\theta$ corresponding to the angle between the direction of
observation and the GC. In the next subsection we will discuss in detail the
distinct DM density profiles considered in our analysis. However, it is
important to emphasize that they will be normalized to the same local DM
density such that $\rho$($r_{\odot}$ = 8.3 kpc) = 0.38 GeV cm-3. This value is
derived from LAMOST DR5 and Gaia DR2 [35]. As a consequence, any change in the
local density can be taken into account by rescaling the expected signal by
($\rho$($r_{\odot}$)/0.38 GeV cm-3).

In Ref. [36], the authors showed the existence of a CR proton accelerator in
the inner 10 pc of the GC, with a harder flux ($\Gamma\simeq$ 2.4) compared to
the flux measured in the local Solar neighborhood. This results in a TeV CR
proton density higher than that of the local one [36, 37]. The proton flux
normalization in regions around the GC is chosen such that the proton density
matches the measured averaged density by H.E.S.S. in the inner 200 pc of the
MW, i.e, $\omega_{\rm p}^{\rm GC}$($\geq$10 TeV) = 6$\times 10^{-3}$ eV cm-3,
such as

$\displaystyle w_{\rm p}^{\rm GC}(E_{\rm p}\geq 10\,\rm TeV)$
$\displaystyle=\frac{4\pi}{c}\int_{\rm 10\,TeV}^{\infty}\frac{d\Phi_{\rm
p}}{dE_{\rm p}}(E_{\rm p})\times E_{\rm p}\times dE_{\rm p}\,.$ (2.3)

This is a factor of about six higher than the local proton density $w_{\rm
p}^{\rm local}$ derived from proton fluxes measured by AMS [38] or DAMPE [39].

Diffusive $\gamma$-ray observations from current and future very-high-energy
gamma-ray observatories, such as H.E.S.S., LHAASO, MAGIC, VERITAS, CTA and
SWGO, can be used to place strong upper limits on the inelastic cross section.
One of our goals is to derive model-independent upper bounds on $\sigma^{\rm
inel}_{\rm p\chi}$. Although we will not assume a specific model to describe
the CR proton-DM interaction, we will investigate how our findings are
dependent on the assumption for the energy dependence of $\sigma^{\rm
inel}_{\rm p\chi}(\sqrt{s})$. Our baseline results will be derived for an
energy-independent cross section, i.e., $\sigma^{\rm inel}_{\rm
p\chi}(\sqrt{s})=\sigma^{\rm inel}_{\rm p\chi}(\sqrt{s_{0}})$, where
$\sqrt{s_{0}}$ is assume to be 1 GeV in what follows. For the energy-dependent
case, inspired by Regge theory [40], which predicts a power-law behavior for
the energy dependence of the proton-proton, meson-proton and real photon-
proton cross sections, we will assume $\sigma_{\rm p\chi}^{\rm
inel}=\sigma_{\rm p\chi}^{\rm
inel}(\sqrt{s_{0}})\times\left(s/s_{0}\right)^{\delta}$. In our analysis, we
will consider two distinct values for the power: $\delta=+0.096$ and
$\delta=-0.45$, which are the values that describe the current $pp$, $\pi p$,
$Kp$ and $\gamma p$ data for high ($\sqrt{s}\gg\sqrt{s_{0}}$) and low
($\sqrt{s}\ll\sqrt{s_{0}}$) energies [23, 41], respectively.

Finally, inspired by Ref. [31], we will investigate if the derived upper
bounds for $\sigma^{\rm inel}_{\rm p\chi}$ improve the current exclusion
limits on the elastic cross section obtained in colliders and direct detection
experiments. In order to relate the inelastic to the elastic cross section, we
assume either the relation predicted by the black disc limit for asymptotic
center of mass energies [40] such that $\sigma_{\rm p\chi}^{\rm
inel}=\sigma_{\rm p\chi}^{\rm el}$, or a relation similar to what is verified
in $pp$ collisions at the LHC [42], $\sigma_{\rm p\chi}^{\rm
inel}=8/3\,\sigma_{\rm p\chi}^{\rm el}$ .

###  2.2 Dark matter distribution in the Milky Way

As described in the previous subsection, the expected gamma-ray signal from
the CR-DM scattering depends on the assumed DM density distribution. The GC
region is an excellent target for DM searches from CR-DM scattering: (i) it is
a region where we expect high DM density as structure formation indicates that
it increases towards the centre of galaxies; and (ii) the CR density is found
to be significantly larger in the GeV-to-TeV energies [43, 36] than in the
local one. The expected signal has a weak dependence on the precise radial DM
distribution in the GC region since it scales linearly with the DM density
$\rho$ as opposed to the case of annihilating DM where it depends on the
square of $\rho$. Although the highest DM density distribution is expected in
the GC region, its precise determination in this complex environment is
challenging. Here we use canonical DM profiles for the description of the DM
density distribution in the Milky Way. The baseline DM density distribution
chosen is parametrized by the Navarro–Frenk–White (NFW) [44] profile given by

$\rho_{\rm NFW}(r)=\frac{\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}}$ (2.4)

with $\rho_{s}$ and $r_{s}$ the scale density and radius, respectively.
Alternative parametrisations are given by the Einasto [45] and Burkert [46]
density profiles, written as

$\rho_{\rm
Ein}(r)=\rho_{s}\exp\left\\{-\frac{2}{\alpha}\left[\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right]\right\\}\quad{\rm
and}\quad\rho_{\rm Bur}(r)=\frac{\rho_{s}}{(1+r/r_{s})(1+(r/r_{s})^{2})}\,,$
(2.5)

respectively, with $\alpha$ = 0.17 being a sharpness index. The values of the
parameters for each density profile are given in Tab. 1.

| Einasto [9] | NFW [47] | cNFW [47] | Burkert [48]
---|---|---|---|---
$\rho_{\rm s}$ (GeVcm-3) | 0.077 | 0.50 | 0.24 | 0.93
$r_{\rm s}$ (kpc) | 20.0 | 15.5 | 23.8 | 12.6
Table 1: DM profile parameters used in this work. Scale density $\rho_{\rm s}$
and radius $r_{\rm s}$ values are given for the Einasto, NFW, cNFW and Burkert
profile parametrizations. All the profiles are normalized such that $\rho$(r=
r⊙) = 0.38 GeVcm-3 assuming r⊙ = 8.3 kpc [35].

In addition to these canonical profiles, we make use a profile parametrization
derived from a determination of the MW mass profile using Gaia DR2
measurements of the rotation curve together with a modeling of the baryonic
components in the GC [49]. The profile is contracted compared to the standard
NFW profile by the presence of baryons. The method followed in Ref. [47] does
not allow for a precise determination within the inner kpc of the MW and we
will conservatively assume a cored density with 1 kpc [47]. Such a profile is
hereafter referred as to the cNFW profile. However, the DM density could still
increase towards the inner GC as shown in several realizations of the DM
distribution of MW-like galaxies from realistic N-body cosmological
simulations including baryonic effects [50, 51]. We extracted a representative
DM distribution from Ref. [51], where various hydrodynamical simulations using
FIRE-2 zoom-in simulations on MW-like galaxies are performed including
radiative heating and cooling for gas, stellar feedback from OB stars, type Ia
and type II supernovae, radiation pressure, and star formation effects. A DM
profile for representative DM-only FIRE-2 simulation is also extracted.

Figure 2: Differential D-factor profiles, expressed in TeV cm-2sr-1 versus the
angular distance from the GC, $\theta$. Various DM profiles are considered:
Einasto (dark blue solid line) , NFW (purple dashed line), cNFW (burgundy
dotted line), Burkert (red dotted-dashed line) as well as those obtained from
FIRE-2 cosmological simulations with (orange dotted line) and without (yellow
dash-double-dotted line) baryonic feedback. See text for more details.

Figure 2 shows the differential D-factors versus angular distance from the GC
for the selected benchmark DM distributions in the inner MW halo. When
integrated within the inner 3∘ of the GC, a factor of $\sim$4.4 is expected
between the D-factor, and hence the signal expectation, computed for the two
extreme cases considered here, i.e., the Einasto and Burkert profiles.

###  2.3 VHE gamma-ray backgrounds in the regions of interest

An irreducible background arises from the limited discrimination of the
instruments considered here between true gamma-ray and cosmic-ray hadrons and
electrons, hereafter referred as the residual background. Its determination
can be performed via careful Monte Carlo simulations of the response of the
instrument or via measurements in blank-field observations in the same
observational conditions as for the signal measurement (see Sec. 3 for more
details). The energy-differential fluxes per solid angle for the residual
background together with the expected signal are shown in Fig. 3 for
H.E.S.S.-like observations as well as CTA, and strawman SWGO expected
observations.

Figure 3: Energy-differential gamma-ray fluxes expected for the DM-induced
signal (solid lines) and residual background (dashed lines), respectively, in
the GC region for H.E.S.S.-like (dark blue lines), CTA (purple lines) and SWGO
(red lines) ROIs. For SWGO, strawman perfomances are extracted from Ref. [52].
The expected signals assume a DM density distribution following the NFW
profile populated by DM particles of mass $m_{\rm\chi}$ = 10 keV, and an
inelastic $p\chi$ cross section taken to be constant such as $\sigma^{\rm
inel}_{\rm p\chi}$ = 10-32 cm2. The energy-differential flux of the GDE for
the SA0 (R12) model and averaged over the inner 3∘ is also plotted (dotted-
dashed yellow line). See text for further details.

The central region of the MW is a very rich and complex region of the sky in
VHE gamma rays. In addition to the bright sources such as the central source
HESS J1745-290 [53, 36], the supernova remnant HESS J1745-303 [54] as well as
various emitters lying close by in the Galactic plane [55], a faint diffuse
emission has been detected, including VHE diffuse emission attributed to PeV
hadrons accelerated in the vicinity of the supermassive black hole Sgr A*
[36]. In order to avoid complex modelling of these emissions, we will mask
Galactic latitudes $|b|<$0.3∘ in the analysis of the inner GC region for
H.E.S.S.-like, CTA and SWGO observations.

Energetic cosmic rays interact with the ambient target materials of the
interstellar medium as well as interstellar photon fields producing a diffuse
flux of gamma rays known as the Galactic diffuse emission (GDE). Such gamma
rays are produced from a combination of Bremsstrahlung from cosmic-ray hadrons
as well as the decay of neutral pions produced from these hadrons colliding
with the interstellar gas, and from cosmic-ray electrons scattering off
ambient photon fields via inverse Compton scattering (ICS). The GDE
contribution dominates the photons seen by Fermi-LAT at MeV-to-GeV energies
[56]. With the current and expected deep observations of the GC region, GDE is
expected to be a relevant background component at TeV energies. In the present
work, we will include a realistic model of the expected TeV-scale GDE
including the $\pi^{0}$ and ICS components from the most recent release of the
CR propagation framework GALPROP (version 57) [57]. We will make use of the
SA100 (R12) model as a model for the GDE predictions with GALPROP extracted
from Ref. [58]. A detailed study of the GDE impact on the sensitivities to
sub-GeV DM is out of the scope of the present paper, however, we will show how
much the sensitivity can be improved when the GDE is subtracted from the
measured residual background for the H.E.S.S.-like and CTA observations. The
expected energy-differential fluxes for the GDE including both the $\pi_{0}$
and ICS components for the SA0 (R12) GDE model is shown in Fig. 3.

##  3 Current and near-future VHE gamma-ray observatories

###  3.1 H.E.S.S.

The High Energy Stereoscopic System (H.E.S.S.) consists of an array of five
Imaging Atmospheric Cherenkov Telescopes (IACTs) located in the Khomas
Highlands of Namibia near the tropic of Capricorn at an altitude of 1800 m
[24]. It is the only operating array of IACTs in the Southern hemisphere.
H.E.S.S. is sensitive to gamma rays from a few ten GeV to several tens TeV and
reaches an energy resolution of $\sigma$/E $\simeq$ 10% and an angular
resolution of 0.06∘ at 68% containment radius.

Its location allows the observation of the inner few degrees of the MW under
favorable conditions. The GC can be observed from March to September each year
with observational zenith angles lower than 30∘. Observations at low zenith
angle are usually preferred since they provide lower energy threshold and
control of observational systematic uncertainties given that atmosphere layers
are thinner along the line of sight. However, observations at high zenith
angles increase the effective area at higher energies to the price of a raise
in the energy threshold and degraded event reconstruction performances and
potentially higher systematic uncertainties (see, for instance, Refs. [59, 60]
for further discussions).

In order to assess the constraints on DM from CR scattering from current IACTs
in the central region of the Milky Way, we will make use of the the most up-
to-date H.E.S.S.-like observations of the inner Galactic halo with the full
five-telescope array [9]. These observations are taken thanks to the Inner
Galaxy Survey program which provides the highest exposure in VHE gamma rays so
far, for a total of 546 hours of high quality data distributed over a few
degrees around the GC. In the following, we consider a region of interest
(ROI) of 3∘ in radius around the GC with an homogeneous time-exposure of 500
hours. The ROI is further divided into 27 concentric rings of 0.1∘ width. The
energy-dependent acceptance and residual background flux for five-telescope
H.E.S.S. observations of the inner halo of the Milky Way is extracted from
Ref. [9]. We make use of 67 logarithmically-spaced energy bins in the 170 GeV
- 67 TeV energy range. This enables us to obtain realistic instrument response
for the region of interest we consider.

###  3.2 LHAASO

The Large High Altitude Air Shower Observatory (LHAASO) is a large air shower
observatory currently operating in southwest China [25], set at an altitude of
4,410 meters. It is composed of the kilometer square array (KM2A), the water
Cherenkov detector array (WCDA), and the wide field of view Cherenkov
telescope array (WFCTA). LHAASO detects gamma rays from 100 GeV up to 1 PeV.
The performances in terms of acceptance, energy and angular resolutions depend
on the considered subarray. In what follows, we will make use of the
observations taken by KM2A. It consists of a full-duty extensive air shower
array with gamma/hadron separation capability by using both electromagnetic
particle detectors and underground muon detectors. At 30 TeV, the angular and
energy resolutions of KM2A are about 0.4∘ (68% containment radius) and 20%,
respectively, for zenith angles below 20∘ [61].

Given its location in the Northern hemisphere, LHAASO is not suited for
observations of the GC region. However, the searched signal is expected to be
higher at smaller galactocentric distances. We will therefore consider the
region of interest around $15^{\circ}<b<45^{\circ}$ and
30${}^{\circ}<l<60^{\circ}$ with a solid angle of 0.274 sr, as described in
Ref. [62]. Due to the coordinates of the region of interest, for LHAASO, the
local flux of CRs (as measured in the Solar neighborhood) will be used in the
signal calculations instead the flux in the GC region. Four alternative ROIs
were considered with the same solid angle and declination for the estimate of
the residual background. In what follows, we will assess the expected
sensitivity of LHASSO for 340 days of observations with 1/2-KM2A and 230 days
with 3/4-KM2A using the energy-dependent acceptance and background rates with
6 logarithmically-spaced energy bins in the 100 TeV -1.5 PeV energy range
provided in Ref [62].

###  3.3 CTA

The Cherenkov Telescope Array (CTA) [28, 29] is the next-generation of IACTs
with one array in the southern hemisphere at the Paranal site (Chile) and one
in the northern hemisphere in La Palma (Spain). It will observe the gamma-ray
sky from 20 GeV up to 300 TeV with a factor up-to-ten increase in sensitivity,
and angular resolution improved by a factor 2-to-3 compared to currently
operating IACTs, and an energy resolution down to as low as 5% at TeV
energies. As the GC region is a prime target in the context of CR-DM
scattering, for this work we will use as baseline the expected performance of
the so-called Omega configuration array in the Southern hemisphere. This full-
scope configuration is expected to be composed of 4 large-sized, 25 medium-
sized, and 70 small-sized telescopes. However, the current initial phase of
the CTAO, at Paranal Observatory, will consist of 14 medium-sized Telescopes
and 37 small-sized Telescopes, known as the Alpha configuration. The
performances of these two configurations will be compared.

We will consider a circular ROI of 5∘ in radius around the GC, excluding
Galactic latitudes between $\pm$0.3∘, in order to avoid complex modeling of
the astrophysical sources in the GC plane region, with an homogeneous time-
exposure of 500 hours. The ROI is further divided into 47 concentric annuli of
0.1∘ width each. The Instrument Response Functions (energy-dependent effective
area, angular and energy resolutions) computed from Monte Carlo simulations
from the Omega configuration array are taken from the publicly available
prod3b-v2 IRF library [63] using 92 logarithmically-spaced energy bins between
50 GeV and 140 TeV. In particular the IRF file South_z20_average_50h is used
as optimized for the detection of a point-like source at zenith angle of 20∘.
As for the Alpha configuration array, the IRFs ate taken from the publicly
available prod5-v0.1 library [64] using the same 92 logarithmically-spaced
energy bins. The file Performance-prod5-v0.1-South-20deg is used as optmized
for the detection of a point-like source at zenith angle 20∘.

###  3.4 SWGO

Southern Wide-field Gamma-ray Observatory (SWGO) is the project to build a
water Cherenkov particle detector with a wide field of view and high duty-
cycle currently in the research and development phase, and planned to be
located in South America at an altitude of 4.4 km or higher and latitude
between 10 and 30 degrees south. This gamma-ray observatory is planned to be
made of detector units which could be deployed as an array or assembled in a
building covering an energy range between a few hundred GeV up to the PeV
scale [65].

In this work, we use the instrument response functions of the Strawman design,
as described in the Science Case white paper for the SWGO [66], which are
publicly available [52]. We consider ten years of observation time with SGWO
looking at the inner 10∘ of the inner MW, excluding Galactic latitudes
$|b|<0.3^{\circ}$ to assess the SWGO sensitivity. The ROI is further divided
into of 49 concentric annuli of 0.2∘ width each, following the procedure of
Ref. [67]. We make use of 46 logarithmically-spaced energy bins in the 200 GeV
- 670 TeV energy range.

##  4 Statistical analysis framework

###  4.1 Expected number of signal and background photons

For a given density profile, the expected number of signal events, $N^{\rm
S}_{\rm ij}$ in the $i$th ROI with a solid angle $\Delta\Omega_{\rm i}$, and
$j$th energy bin having a width $\Delta E_{\rm j}$, is given by

$N^{S}_{ij}=T_{{\rm obs},i}\int_{E_{j}-\Delta E_{j}/2}^{E_{j}+\Delta
E_{j}/2}dE\int_{-\infty}^{\infty}dE^{\prime}\,\frac{d\Phi^{S}_{ij}}{dE^{\prime}}(\Delta\Omega_{i},E^{\prime})\,A_{\rm
eff}^{\gamma}(E^{\prime})\,G(E-E^{\prime}).$ (4.1)

Here E and E’ are the true and reconstructed energies, respectively,
d$\Phi^{S}_{ij}$/dE’ is the predicted gamma-ray flux from CR-DM scattering, as
given in Eq. (2.1), Tobs,i corresponds to the observation time in the $i$th
ROI, A${}^{\gamma}_{\rm eff}$ is the energy-dependent effective area to gamma
rays, and $G$ stands for the finite energy resolution of the instrument which
is modeled as a Gaussian with width specified by $\sigma$/E. The number of
background events can be obtained straightforwardly from Eq. (4.1) from the
substitution of $d\Phi_{ij}^{\rm S}/dE$ by the energy-differential background
flux, $d\Phi_{ij}^{\rm B}/dE$, for both residual and GDE backgrounds.

The IRFs used for each observatory are mentioned in Sec. 3. They enable us to
provide a realistic description of the survey we consider for the computation
of the sensitivity. Homogeneous exposure over all the ROIs are considered
unless otherwise stated. Given the limited field of view of IACTs, a pointing
strategy should be defined with a grid of pointing positions in order to reach
such an exposure. However, the optimization of the pointing strategy is beyond
the scope of this work.

###  4.2 Analysis technique and sensitivity computation

The statistical analysis and the computation of the expected sensitivity of
the observatories considered here (see Sec. 3) are performed following the
definition of a log-likelihood ratio test statistic, a well-establish
procedure in VHE gamma rays (see, for instance, Refs. [68, 69, 9]). The TS
exploits the expected spectral and spatial features of the searched signal.
The 2-dimensional binned likelihood function for a given comparison between a
theoretical model and observations writes as

$\mathcal{L}_{ij}(N_{ij}^{S},N_{ij}^{B},\bar{N}_{ij}^{S},\bar{N}_{ij}^{B}\,|\,N_{ij}^{\rm
ON},N_{ij}^{\rm OFF})=\textrm{Pois}[(N^{S}_{ij}+N^{B}_{ij}),N^{\rm
ON}_{ij}]\times\textrm{Pois}[(\bar{N}^{S}_{ij}+\bar{N}^{B}_{ij}),N^{\rm
OFF}_{ij}]\,,$ (4.2)

with $\textrm{Pois}[\lambda,n]=e^{-\lambda}\lambda^{n}/n!$. $N^{\rm ON}_{ij}$
and $N^{\rm OFF}_{ij}$ are the number of measured events in the ON and OFF
regions, respectively, in the spectral bin $i$ and the spatial bin $j$.
$N^{S}_{ij}+N^{B}_{ij}$ and $\bar{N}^{S}_{ij}+\bar{N}^{B}_{ij}$ denote the
expected numbers of events in the signal and background regions, respectively.
The background is assumed to be measured from extragalactic blank-field
observations or computed from Monte Carlo simulations (see Sec. 3) such that
$\bar{N}^{S}_{ij}$ can be considered negligible compared to $N^{S}_{ij}$.

The full likelihood function is given by the product of the likelihood
functions (Eq. (4.2)) over all spatial and energy bins. The remaining nuisance
parameter is the background rate and it is estimated by maximization of the
profile likelihood. Therefore, the likelihood will depend only on the DM
particle properties, which are parameterized here by $\sigma^{\rm inel}_{\rm
p\chi}$ and mχ. For a given DM mass mχ, the likelihood is then used to define
the log-likelihood ratio test statistic (TS) for setting upper limits such as

$\text{TS}(m_{\chi})=-2\ln\frac{\mathcal{L}(\sigma^{\rm inel}_{\rm
p\chi},m_{\chi})}{\mathcal{L}(\widehat{\sigma^{\rm inel}_{\rm
p\chi}},m_{\chi})}\,,$ (4.3)

where $\widehat{\sigma^{\rm inel}_{\rm p\chi}}$ corresponds to the value of
the inelastic cross section which maximizes the likelihood for a given
$m_{\chi}$. In the large-statistics limit, this TS follows a $\chi^{2}$
distribution with a one degree of freedom. Assuming this limit holds, we can
set one-sided 95% upper limits on $\sigma^{\rm inel}_{\rm p\chi}$ by solving
Eq. (4.3) for the cross section above the best fit where TS = 2.71. This is
the procedure we use in the following to compute our limits.

Figure 4: Log-likelihood ratio TS profile curves for H.E.S.S.-like (dark blue
solid line), CTA (purple dashed line), strawman SWGO (orange dotted line) and
LHAASO (yellow dotted-dashed line) observations as a function of the inelastic
cross section $\sigma^{\rm inel}_{\rm p\chi}$ assuming a NFW DM profile and a
DM mass of 1 MeV. The black dashed horizontal line represents TS = 2.71 used
to derive 95% C.L. value on $\sigma^{\rm inel}_{\rm p\chi}$.

The mean expected upper limits and confidence intervals can be computed by
applying the above-mentioned procedure to a large number of Monte Carlo
realizations of mock datasets of $N^{\rm ON}_{ij}$ and $N^{\rm OFF}_{ij}$ and
analyzing the distribution of the obtained limits. This procedure is largely
used in IACT analyses, and in particular in H.E.S.S. [68, 69, 9]. However, in
order to avoid the computational intensive generation of datasets, one can
take the mean expected background as a single ‘representation’ of the dataset,
known as the Asimov dataset [70], in order to compute the mean expected
sensitivity. The Asimov procedure can be also used to compute the $N$-sigma
containment interval [70] from ${\rm TS}=(\Phi^{-1}[0.95]\pm N)^{2}$, where
$\Phi^{-1}$ is the inverse of the cumulative distribution function of a normal
distribution with mean $\mu=0$ and standard deviation $\sigma=1$. We will use
the Asimov procedure throughout this work for the derivation of the mean
expected limits together with their 1$\sigma$ and 2$\sigma$ containment bands.
Figure 4 shows the TS profile obtained for the different observatories as a
function of the free parameter $\sigma^{\rm inel}_{\rm p\chi}$ for a given DM
mass of 1 MeV and a NFW DM density distribution profile for the MW.

Figure 5: Sensitivity curves for H.E.S.S.-like (blue line), CTA (purple
lines), strawman SWGO (orange line) and LHAASO (yellow line) observations on
inelastic scattering of CR protons off DM. The solid and dashed purple lines
corresponds to the sensitivity expected for the Omega (baseline) and Alpha
configurations, respectively. The sensitivity is expressed as 95% C.L. mean
expected upper limit. Here we assumed a NFW DM profile and an energy-
independent inelastic cross section normalized at $\sqrt{s}$ = 1 GeV. The
sensitivity reach for H.E.S.S.-like observations when the GDE component is
subtracted is shown as purple dashed line.

##  5 Results and discussions

The sensitivity is expressed in terms of mean expected upper limit computed at
95% C. L. (see Sec. 4.2). We choose the NFW profile as the baseline DM density
profile and the fraction of the proton energy going into neutral pions to be
$f_{\pi}$ = 0.17, excepted when explicitly stated otherwise.

Figure 5 shows the sensitivity for H.E.S.S.-like, CTA (Omega and Alpha
configurations), SWGO and LHAASO, on the inelastic cross section versus DM
mass, assuming a NFW profile and an energy-independent inelastic cross section
(see Sec. 2.1). The strongest sensitivities reached for each observatory are
the following: for H.E.S.S.-like, $\sigma_{\rm p\chi}^{\rm inel}$ =
5.3$\times$10-32 cm2 for a mass $m_{\rm\chi}$ = 108 keV; for CTA, $\sigma_{\rm
p\chi}^{\rm inel}$ = 1.6$\times$10-33 cm2 for a mass $m_{\rm\chi}$ = 13 keV;
for SWGO, $\sigma_{\rm p\chi}^{\rm inel}$ = 1.0$\times$10-34 cm2 for a mass
$m_{\rm\chi}$ = 1.4 keV; and for LHAASO, $\sigma_{\rm p\chi}^{\rm inel}$ =
2.0$\times$10-33 cm2 for a mass $m_{\rm\chi}$ = 0.8 keV. The expected
sensitivity from CTA Alpha configuration (purple dashed line) shows a reduced
sensitivity by a factor up to 3 in the 10 - 100 keV mass range compared to the
Omega configuration sensitivity. Above 100 keV masses, the Alpha CTA
configuration is about 50% more sensitive than the baseline configuration due
to an increased ratio between signal and square root of the expected
background. In order to understand why the minimum values of $m_{\chi}$
reached for these observatories are different, it is important to remember
that the gamma-ray production only occurs if the CR proton has enough energy
to produce (at least) one neutral pion in the interaction with a DM of mass
$m_{\chi}$. Such threshold defines the value of $E_{\rm p}^{\rm min}$, which
is inversely proportional to $m_{\chi}$, and also $E_{\gamma}^{\rm min}$ [see
discussion after Eq. (2.1)]. As the accessible gamma-ray energy ranges of the
experiments are distinct, it implies that they probe different minimum values
of $m_{\chi}$. In particular, since LHAASO is capable of detecting the
highest-energy gamma rays, it can probe the smallest DM masses. However, its
sensitivity above 100 keV masses is lower than that of the other experiments
due to the smaller CR flux in its ROIs compared to the flux in the GC region.
This decrease in sensitivity occurs despite the fact that the observed region
and effective area are larger compared to other experiments. When subtracting
the GDE component to the overall background, the H.E.S.S.-like sensitivity
improves by 10%. A slight dependency is obtained with the DM mass. The
residual background becomes significantly larger than the GDE component at the
highest energies, which therefore implies that the GDE impact becomes
negligible at the lowest masses.

Figure 6 shows the effect of the energy dependency for the inelastic $p\chi$
cross section on the sensitivity curve for H.E.S.S.-like observations. The
solid line represents the energy-independent inelastic cross section between
DM and CR protons, respectively. For the energy-dependent case, the dashed and
dotted-dashed lines correspond to the $\propto E^{0.096}$ and $\propto
E^{-0.45}$ cases, respectively. The energy dependence of the inelastic cross
section induces a specific behaviour of the mass-dependency of the
sensitivity, reaching a factor up to $\sim$ 50 for a DM mass of 100 MeV. In
the energy-independent case, the sensitivity scales approximately with
$m_{\chi}$ in the high-mass limit.

Figure 6: Comparison between the sensitivity reached for H.E.S.S.-like
observations for the energy-dependent and independent inelastic cross sections
$\sigma^{\rm inel}_{\rm p\chi}$, respectively. The sensitivity is shown as a
solid purple line for the constant cross section, and a band between the
energy-dependent power law cases for the inelastic cross section with spectral
indexes of 0.096 (increasing power law, orange dashed line) and -0.45
(decreasing power law, dotted-dashed orange line), respectively. At the lowest
masses, it implies a difference of 60% while a factor of up to $\sim 50$ is
obtained at the highest masses.

Figure 7:  H.E.S.S.-like sensitivity expressed as 95% C.L. mean upper limit
versus DM mass for the elastic CR - DM scattering process. Left panel:
H.E.S.S.-like sensitivity expressed in terms of 95% C.L. mean expected upper
limit for different DM density profiles considered here: Einasto (dark blue
solid line), NFW (purple dashed line), cNFW (burgundy dotted line), Burkert
(red dotted-dashed line), Fire-2 (orange loosely dotted line) and Fire-2 - DMO
(yellow up-ticked line), respectively. The difference between the strongest
(Einasto profile) and weakest (Burkert profile) curves is given by a
multiplicative factor of $\sim 4$. Right panel: The sensitivity expressed as
the mean expected upper limit computed at 95% C.L. (dark blue solid line) is
shown for the H.E.S.S.-like observations together with the included 68% (green
area) and 95% (yellow area) confidence level bands, respectively, assuming the
NFW DM profile,

In what follows, we investigate if the constraints derived for $\sigma^{\rm
inel}_{\rm p\chi}$ can improve the current limits in the elastic proton - DM
cross section, which are usually derived considering direct detection and
cosmological bounds. As explained in Sec. 2, we assume that the relation
between these cross sections can be expressed by $\sigma^{\rm inel}_{\rm
p\chi}=\kappa\cdot\sigma^{\rm el}_{\rm p\chi}$, with $\kappa=8/3$ or
$\kappa=1$, motivated by the current results for $pp$ collisions at the LHC
and by the black disc limit, respectively. The left panel of Fig. 7 shows our
results for $\sigma^{\rm el}_{\rm p\chi}$ on the H.E.S.S.-like sensitivity,
derived for $\kappa=8/3$, considering different DM density profiles. A factor
of about 4 is obtained between the strongest sensitivity obtained for the
Einasto profile, and the weakest one for the Burkert profile. The right panel
of Fig. 7 shows the sensitivity expressed in terms of mean expected upper
limits computed at 95% C. L. together with the 68% and 95% confidence level
statistical containment bands. We have verified that these results are
modified by a factor 2.6 if the value $\kappa=1$ is considered in the
calculations. Moreover, we also have investigated the impact of a smaller
value of $f_{\pi}$, which can occur with the increasing of the neutral pion
multiplicity at high energies. Assuming a decrease of $f_{\pi}$ by a factor 3
would degrade the sensitivity by a factor of $\sim$4\. Decreasing $f_{\rm\pi}$
would modify the accessible mass range to higher masses. We have verified that
the dependence of the expected gamma-ray flux on $\kappa$ and $f_{\pi}$
translates into a straightforward scaling of the sensitivity [see Eq. (2.1)].
A change by a factor of about 9 is obtained between the scenario with $\kappa$
= 8/3 and $f_{\pi}$ = 0.17, and the scenario with $\kappa$ = 1 and $f_{\pi}$ =
0.17/3 at larger masses. Surely, the uncertainty on these predictions can be
reduced by assuming a more detailed description of the CR proton-DM
interaction and for the photon production in the final state, which we plan to
perform in a forthcoming study.

In order to put these results into a wider context and compare them with
direct detection, colliders and cosmological constraints, we will use the
energy-independent case for the inelastic cross section as well as the
relationship between the inelastic and elastic cross sections, discussed in
Sec. 2.1, given by a factor 8/3. Figure 8 shows the sensitivity for sub-GeV DM
computed for the observatories considered here together with the excluded
regions from other probes such as Cosmological, using CMB [71, 72, 73],
Lyman-$\alpha$ [74] and Milk Way subhalo (MWS) [74] measurements, LHC [75],
and direct detection [76, 13, 77] searches. The sensitivity reached for each
experiment goes as follows: for H.E.S.S.-like, $\sigma_{\rm p\chi}^{\rm el}$ =
1.9$\times$10-32 cm2 for a mass $m_{\rm\chi}=$ 108 keV; for CTA, $\sigma_{\rm
p\chi}^{\rm el}$ = 6.3$\times$10-34 cm2 for a mass $m_{\rm\chi}=$ 13 keV; for
a strawman SWGO model, $\sigma_{\rm p\chi}^{\rm el}$ = 3.9$\times$10-32 cm2
for a mass $m_{\rm\chi}=$ 1.4 keV; for LHAASO, $\sigma_{\rm p\chi}^{\rm el}$ =
7.6$\times$10-34 cm2 for a mass $m_{\rm\chi}=$ 0.8 keV.

Figure 8: Current constraints on the elastic cross section between DM and
protons. VHE gamma-ray constraints assuming no energy dependence of the
elastic cross-section (see Sec. 2.1) obtained in the present work are shown
for current (H.E.S.S-like and LHAASO) and future (CTA Omega, strawman SWGO)
VHE gamma-ray observatories. Direct detection experiment constraints derived
by CRESST [76] and Xenon1T [13, 22], neutrino searches with Super-K [77], and
colliders [75] searches are shown as well as from cosmological probes such, as
CMB [71], Lyman-$\alpha$ [78, 72, 74], and Milky-Way subhalos (MWS) [79, 73,
74].

This highlights the potential of this novel approach to probe a region of the
parameter space currently unavailable to the more traditional collider and
direct detection searches. Using a set of conservative values for $f_{\pi}$
and $\kappa$ would degrade the sensitivity up to one order of magnitude,
leaving this method well positioned to cover uncharted regions of the
parameter space. We forecast here the sensitivity of current on ongoing VHE
observatories in a model-independent way using approximate relationship
between the inelastic and elastic cross section. A more accurate comparison
would require the definition of an underlying model for the dark matter nature
and properties such as coupling strength and mediator mass.

##  6 Summary

In this work, we have derived predictions for the constraints on the inelastic
cross section between CR protons and DM particles populating the DM halo of
the Milky Way based on the latest H.E.S.S. observations of the inner 3∘ of the
MW, and LHAASO observations in the region of the sky defined by
15${}^{\circ}<$ b $<$ 45∘ and 30${}^{\circ}<$ l $<$ 60∘ using the PeV CR
proton source located in the GC and realistic DM distributions in the Milky
Way. The impact of our modelling assumptions is investigated considering a
variety of viable DM profiles and distinct assumptions for the energy
dependence of the inelastic $p\chi$ cross section. Our results indicate that
VHE gamma-ray observations enable to probe a DM mass range poorly covered so
far and provide model - independent upper bounds on the value of $\sigma^{\rm
inel}_{\rm p\chi}$. The strongest limits come from SWGO in the few keV range,
reaching an inelastic cross section of 1.4$\times$10-34 for a DM mass of 1.5
keV. The sensitivity for H.E.S.S.-like observations probe inelastic cross
sections down to $5.3\times 10^{-32}$ cm2 for a 108 keV DM mass. We also have
forecast the sensitivity reach of the future observations of CTA in the
central region of the Milky Way, which demonstrates the potential of these
observatories.

The constraints derived on the inelastic cross section $\sigma^{\rm
inel}_{p\chi}$, assuming no specific underlying DM model, have been recast in
terms of the elastic cross section $\sigma^{\rm el}_{\chi p}$ and upper bounds
for this quantity have been presented considering the distinct VHE gamma-ray
observatories. Moreover, our results have been compared with other probes such
as cosmological, collider and direct detection searches, which has allowed us
to highlight further synergies of the present approach. These results indicate
that current and future VHE gamma-ray observations open a new window for
constraining inelastic and elastic cross sections in the poorly charted DM
mass range from hundred eV to hundred MeV. Further detailed comparisons
between the various probes in this DM mass range will require the choice of a
specific DM model, which will be carried on in a future work, for the
derivation of the relationship between inelastic and elastic cross section
such as couplings, mediator mass, and properties of the DM particles.

## Acknowledgments

We warmly thank Alessandro Montanari for providing GDE emission fits data of
Figures 3 and 4 of Ref. [58] for the SA10 (R12) model. This work was conducted
in the context of the CTA Consortium111https://www.cta-
observatory.org/consortium_acknowledgments/ and the SWGO Collaboration. This
research has made use of the CTA instrument response functions provided by the
CTA Consortium and Observatory, see https://www.cta-
observatory.org/science/cta-performance for more details, and the SWGO
instrument response functions provided by the SWGO Consortium (see [52]). This
work is supported by the "ADI 2021" project funded by the IDEX Paris
ANR-11-IDEX-0003-02, and by FAPESP, process numbers 2019/14893-3, 2021/02027-0
and 2021/01089-1. A.V. is supported by CNPq grant 314955/2021-6. V.P.G. was
partially supported by CNPq, FAPERGS and INCT-FNA (Process No. 464898/2014-5).
A.V. and I.R. acknowledge the National Laboratory for Scientific Computing
(LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer
(http://sdumont.lncc.br).

## References

  * [1] E.W. Kolb and M.S. Turner, _The Early Universe_ , vol. 69 (1990), 10.1201/9780429492860.
  * [2] D. Clowe, M. Bradac, A.H. Gonzalez, M. Markevitch, S.W. Randall, C. Jones et al., _A direct empirical proof of the existence of dark matter_ , _Astrophys. J. Lett._ 648 (2006) L109 [astro-ph/0608407].
  * [3] K. Garrett and G. Duda, _Dark Matter: A Primer_ , _Adv. Astron._ 2011 (2011) 968283 [1006.2483].
  * [4] S. Dodelson, _The Real Problem with MOND_ , _Int. J. Mod. Phys. D_ 20 (2011) 2749 [1112.1320].
  * [5] F. Kahlhoefer, _Review of LHC Dark Matter Searches_ , _Int. J. Mod. Phys. A_ 32 (2017) 1730006 [1702.02430].
  * [6] M. Schumann, _Direct Detection of WIMP Dark Matter: Concepts and Status_ , _J. Phys. G_ 46 (2019) 103003 [1903.03026].
  * [7] L.E. Strigari, _Dark matter in dwarf spheroidal galaxies and indirect detection: a review_ , _Rept. Prog. Phys._ 81 (2018) 056901 [1805.05883].
  * [8] Fermi-LAT collaboration, _Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data_ , _Phys. Rev. Lett._ 115 (2015) 231301 [1503.02641].
  * [9] H.E.S.S. collaboration, _Search for Dark Matter Annihilation Signals in the H.E.S.S. Inner Galaxy Survey_ , _Phys. Rev. Lett._ 129 (2022) 111101 [2207.10471].
  * [10] XENON collaboration, _First Dark Matter Search with Nuclear Recoils from the XENONnT Experiment_ , _Phys. Rev. Lett._ 131 (2023) 041003 [2303.14729].
  * [11] XENON collaboration, _Emission of single and few electrons in XENON1T and limits on light dark matter_ , _Phys. Rev. D_ 106 (2022) 022001 [2112.12116].
  * [12] M. Battaglieri et al., _US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report_ , in _U.S. Cosmic Visions: New Ideas in Dark Matter_ , 7, 2017 [1707.04591].
  * [13] T. Bringmann and M. Pospelov, _Novel direct detection constraints on light dark matter_ , _Phys. Rev. Lett._ 122 (2019) 171801 [1810.10543].
  * [14] J. Alvey, N. Sabti, V. Tiki, D. Blas, K. Bondarenko, A. Boyarsky et al., _New constraints on the mass of fermionic dark matter from dwarf spheroidal galaxies_ , _Mon. Not. Roy. Astron. Soc._ 501 (2021) 1188 [2010.03572].
  * [15] A. Boyarsky, O. Ruchayskiy and D. Iakubovskyi, _A Lower bound on the mass of Dark Matter particles_ , _JCAP_ 03 (2009) 005 [0808.3902].
  * [16] C. Boehm, M.J. Dolan and C. McCabe, _A Lower Bound on the Mass of Cold Thermal Dark Matter from Planck_ , _JCAP_ 08 (2013) 041 [1303.6270].
  * [17] N. Sabti, J. Alvey, M. Escudero, M. Fairbairn and D. Blas, _Refined Bounds on MeV-scale Thermal Dark Sectors from BBN and the CMB_ , _JCAP_ 01 (2020) 004 [1910.01649].
  * [18] P.F. Depta, M. Hufnagel, K. Schmidt-Hoberg and S. Wild, _BBN constraints on the annihilation of MeV-scale dark matter_ , _JCAP_ 04 (2019) 029 [1901.06944].
  * [19] D.E. Soper, M. Spannowsky, C.J. Wallace and T.M.P. Tait, _Scattering of Dark Particles with Light Mediators_ , _Phys. Rev. D_ 90 (2014) 115005 [1407.2623].
  * [20] G. Guo, Y.-L.S. Tsai, M.-R. Wu and Q. Yuan, _Elastic and Inelastic Scattering of Cosmic-Rays on Sub-GeV Dark Matter_ , _Phys. Rev. D_ 102 (2020) 103004 [2008.12137].
  * [21] X. Chu, J.-L. Kuo and J. Pradler, _Dark sector-photon interactions in proton-beam experiments_ , _Phys. Rev. D_ 101 (2020) 075035 [2001.06042].
  * [22] J. Alvey, T. Bringmann and H. Kolesova, _No room to hide: implications of cosmic-ray upscattering for GeV-scale dark matter_ , _JHEP_ 01 (2023) 123 [2209.03360].
  * [23] S. Donnachie, H.G. Dosch, O. Nachtmann and P. Landshoff, _Pomeron physics and QCD_ , vol. 19, Cambridge University Press (12, 2004).
  * [24] https://www.mpi-hd.mpg.de/HESS/.
  * [25] https://english.ihep.cas.cn/lhaaso.
  * [26] http://www.magic.iac.es.
  * [27] https://veritas.sao.arizona.edu.
  * [28] https://www.cta-observatory.org.
  * [29] CTA Consortium collaboration, _Science with the Cherenkov Telescope Array_ , WSP (11, 2018), 10.1142/10986, [1709.07997].
  * [30] https://www.swgo.org/SWGOWiki/doku.php.
  * [31] R.H. Cyburt, B.D. Fields, V. Pavlidou and B.D. Wandelt, _Constraining strong baryon dark matter interactions with primordial nucleosynthesis and cosmic rays_ , _Phys. Rev. D_ 65 (2002) 123503 [astro-ph/0203240].
  * [32] F.A. Aharonian and A.M. Atoyan, _On the emissivity of $\pi$0̂-̂decay gamma radiation in the vicinity of accelerators of galactic cosmic rays._, _Astron. Astrophys._ 309 (1996) 917.
  * [33] F.W. Stecker, _Cosmic gamma rays_ (1971).
  * [34] T.K. Gaisser, R. Engel and E. Resconi, _Cosmic Rays and Particle Physics: 2nd Edition_ , Cambridge University Press (6, 2016).
  * [35] R. Guo, C. Liu, S. Mao et al., _Measuring the local dark matter density with LAMOST DR5 and Gaia DR2_ , _Mon. Not. Roy. Astron. Soc._ 495 (2020) 4828 [2005.12018].
  * [36] H.E.S.S. collaboration, _Acceleration of petaelectronvolt protons in the Galactic Centre_ , _Nature_ 531 (2016) 476 [1603.07730].
  * [37] D. Gaggero, D. Grasso, A. Marinelli, M. Taoso and A. Urbano, _Diffuse cosmic rays shining in the Galactic center: A novel interpretation of H.E.S.S. and Fermi-LAT $\gamma$-ray data_, _Phys. Rev. Lett._ 119 (2017) 031101 [1702.01124].
  * [38] AMS collaboration, _Precision Measurement of the Proton Flux in Primary Cosmic Rays from 1 GV to 1.8 TV with the Alpha Magnetic Spectrometer on the International Space Station_ , in _25th European Cosmic Ray Symposium_ , 12, 2016 [1612.08562].
  * [39] DAMPE collaboration, _Measurement of the cosmic p+He energy spectrum from 46 GeV to 316 TeV with the DAMPE space mission_ , 2304.00137.
  * [40] P.D.B. Collins, _An Introduction to Regge Theory and High-Energy Physics_ , Cambridge Monographs on Mathematical Physics, Cambridge Univ. Press, Cambridge, UK (5, 2009), 10.1017/CBO9780511897603.
  * [41] M.J. Menon and P.V.R.G. Silva, _A study on analytic parametrizations for proton–proton cross-sections and asymptotia_ , _J. Phys. G_ 40 (2013) 125001 [1305.2947].
  * [42] M. Broilo, V.P. Gonçalves and P.V.R.G. Silva, _Model of diffractive excitation in $pp$ collisions at high energies_, _Phys. Rev. D_ 101 (2020) 074034 [2003.04768].
  * [43] Fermi-LAT collaboration, _Development of the Model of Galactic Interstellar Emission for Standard Point-Source Analysis of Fermi Large Area Telescope Data_ , _Astrophys. J. Suppl._ 223 (2016) 26 [1602.07246].
  * [44] J.F. Navarro, C.S. Frenk and S.D.M. White, _The Structure of cold dark matter halos_ , _Astrophys. J._ 462 (1996) 563 [astro-ph/9508025].
  * [45] J. Einasto, _On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters_ , _Trudy Astrofizicheskogo Instituta Alma-Ata_ 5 (1965) 87.
  * [46] A. Burkert, _The Structure of dark matter halos in dwarf galaxies_ , _Astrophys. J. Lett._ 447 (1995) L25 [astro-ph/9504041].
  * [47] A. Montanari, E. Moulin and N.L. Rodd, _Toward the ultimate reach of current imaging atmospheric Cherenkov telescopes and their sensitivity to TeV dark matter_ , _Phys. Rev. D_ 107 (2023) 043028 [2210.03140].
  * [48] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci et al., _PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection_ , _JCAP_ 03 (2011) 051 [1012.4515].
  * [49] M. Cautun, A. Benitez-Llambay, A.J. Deason, C.S. Frenk, A. Fattahi, F.A. Gómez et al., _The Milky Way total mass profile as inferred from Gaia DR2_ , _Mon. Not. Roy. Astron. Soc._ 494 (2020) 4291 [1911.04557].
  * [50] P.F. Hopkins et al., _FIRE-2 Simulations: Physics versus Numerics in Galaxy Formation_ , _Mon. Not. Roy. Astron. Soc._ 480 (2018) 800 [1702.06148].
  * [51] D. McKeown, J.S. Bullock, F.J. Mercado, Z. Hafen, M. Boylan-Kolchin, A. Wetzel et al., _Amplified J-factors in the Galactic Centre for velocity-dependent dark matter annihilation in FIRE simulations_ , _Mon. Not. Roy. Astron. Soc._ 513 (2022) 55 [2111.03076].
  * [52] https://github.com/harmscho/SGSOSensitivity.
  * [53] H.E.S.S. collaboration, _Spectrum and variability of the Galactic Center VHE gamma-ray source HESS J1745-290_ , _Astron. Astrophys._ 503 (2009) 817 [0906.1247].
  * [54] H.E.S.S. collaboration, _Exploring a SNR/Molecular Cloud Association Within HESS J1745-303_ , _Astron. Astrophys._ 483 (2008) 509 [0803.2844].
  * [55] HESS collaboration, _The H.E.S.S. Galactic plane survey_ , _Astron. Astrophys._ 612 (2018) A1 [1804.02432].
  * [56] Fermi-LAT collaboration, _Fermi-LAT Observations of the Diffuse Gamma-Ray Emission: Implications for Cosmic Rays and the Interstellar Medium_ , _Astrophys. J._ 750 (2012) 3 [1202.4039].
  * [57] T.A. Porter, G. Johannesson and I.V. Moskalenko, _The GALPROP Cosmic-ray Propagation and Nonthermal Emissions Framework: Release v57_ , _Astrophys. J. Supp._ 262 (2022) 30 [2112.12745].
  * [58] A. Montanari, O. Macias and E. Moulin, _TeV gamma-ray sensitivity to velocity-dependent dark matter models in the Galactic Center_ , _Phys. Rev. D_ 108 (2023) 083027 [2309.09691].
  * [59] VERITAS collaboration, _VERITAS Observations of the Galactic Center Region at Multi-TeV Gamma-Ray Energies_ , _Astrophys. J._ 913 (2021) 115 [2104.12735].
  * [60] MAGIC collaboration, _Search for Gamma-Ray Spectral Lines from Dark Matter Annihilation up to 100 TeV toward the Galactic Center with MAGIC_ , _Phys. Rev. Lett._ 130 (2023) 061002 [2212.10527].
  * [61] LHAASO collaboration, _Observation of the crab nebula with lhaaso-km2a - a performance study_ , _Chinese Physics C_ 45 (2021) 025002.
  * [62] LHAASO collaboration, _Constraints on heavy decaying dark matter from 570 days of lhaaso observations_ , _Phys. Rev. Lett._ 129 (2022) 261103.
  * [63] C.T.A. Observatory and C.T.A. Consortium, _CTAO Instrument Response Functions - version prod3b-v2_ , Aug., 2021.  10.5281/zenodo.5163273.
  * [64] C.T.A. Observatory and C.T.A. Consortium, _CTAO Instrument Response Functions - prod5 version v0.1_ , Sept., 2021.  10.5281/zenodo.5499840.
  * [65] SWGO collaboration, _The Southern Wide-field Gamma-ray Observatory_ , _PoS_ ICRC2023 (2023) 963 [2309.04577].
  * [66] A. Albert et al., _Science Case for a Wide Field-of-View Very-High-Energy Gamma-Ray Observatory in the Southern Hemisphere_ , 1902.08429.
  * [67] A. Viana, H. Schoorlemmer, A. Albert, V. de Souza, J.P. Harding and J. Hinton, _Searching for Dark Matter in the Galactic Halo with a Wide Field of View TeV Gamma-ray Observatory in the Southern Hemisphere_ , _JCAP_ 12 (2019) 061 [1906.03353].
  * [68] H.E.S.S. collaboration, _Search for dark matter annihilations towards the inner Galactic halo from 10 years of observations with H.E.S.S_ , _Phys. Rev. Lett._ 117 (2016) 111301 [1607.08142].
  * [69] H.E.S.S. collaboration, _Search for $\gamma$-Ray Line Signals from Dark Matter Annihilations in the Inner Galactic Halo from 10 Years of Observations with H.E.S.S._, _Phys. Rev. Lett._ 120 (2018) 201101 [1805.05741].
  * [70] G. Cowan, K. Cranmer, E. Gross and O. Vitells, _Asymptotic formulae for likelihood-based tests of new physics_ , _European Physical Journal C_ 71 (2011) 1554 [1007.1727].
  * [71] V. Gluscevic and K.K. Boddy, _Constraints on Scattering of keV–TeV Dark Matter with Protons in the Early Universe_ , _Phys. Rev. Lett._ 121 (2018) 081301 [1712.07133].
  * [72] W.L. Xu, C. Dvorkin and A. Chael, _Probing sub-GeV Dark Matter-Baryon Scattering with Cosmological Observables_ , _Phys. Rev. D_ 97 (2018) 103530 [1802.06788].
  * [73] T.R. Slatyer and C.-L. Wu, _Early-Universe constraints on dark matter-baryon scattering and their implications for a global 21 cm signal_ , _Phys. Rev. D_ 98 (2018) 023013 [1803.09734].
  * [74] M.A. Buen-Abad, R. Essig, D. McKeen and Y.-M. Zhong, _Cosmological constraints on dark matter interactions with ordinary matter_ , _Phys. Rept._ 961 (2022) 1 [2107.12377].
  * [75] N. Daci, I. De Bruyn, S. Lowette, M.H.G. Tytgat and B. Zaldivar, _Simplified SIMPs and the LHC_ , _JHEP_ 11 (2015) 108 [1503.05505].
  * [76] J.I. Collar, _Search for a nonrelativistic component in the spectrum of cosmic rays at Earth_ , _Phys. Rev. D_ 98 (2018) 023005 [1805.02646].
  * [77] Super-Kamiokande collaboration, _Search for Cosmic-Ray Boosted Sub-GeV Dark Matter Using Recoil Protons at Super-Kamiokande_ , _Phys. Rev. Lett._ 130 (2023) 031802 [2209.14968].
  * [78] R. Murgia, A. Merle, M. Viel, M. Totzauer and A. Schneider, _”Non-cold” dark matter at small scales: a general approach_ , _JCAP_ 11 (2017) 046 [1704.07838].
  * [79] E.O. Nadler, V. Gluscevic, K.K. Boddy and R.H. Wechsler, _Constraints on Dark Matter Microphysics from the Milky Way Satellite Population_ , _Astrophys. J. Lett._ 878 (2019) 32 [1904.10000].