Multi-TeV dark matter density in the inner Milky Way halo: spectral and dynamical constraints

We develop a comprehensive study of the gamma-ray flux observed by the High Energy Stereoscopic System (H.E.S.S.) in 5 regions of the Galactic Center (GC). Motivated by previous works on a possible Dark Matter (DM) explanation for the TeV cut-off observed by H.E.S.S. in the innermost ∼ 15 pc of the Galaxy, we aim to constrain the DM distribution up to a radius of ∼ 450 pc from the GC. In this region, the benchmark approach (e.g. cosmological simulations and Galactic dynamics studies) fails to produce a strong prediction of the DM profile. Within our proof-of-concept analysis, we use DRAGON to model the diffuse background emission and determine upper limits on the density distribution of thermal multi-TeV Weakly Interactive Massive Particles (WIMPs), compatible with the observed gamma-ray flux. The results are in agreement with the hypothesis of an enhancement of the DM density in the GC with respect to the benchmark Navarro-Frenk-White (NFW) profile (γ = 1) and allow us to exclude profiles with an inner slope cuspier than γ ≳ 1.3. We also investigate the possibility that such an enhancement could be related to the existence of a DM spike associated with the supermassive black hole Sgr A* at the GC. We find out that the existence of an adiabatic DM spike smoothed by the scattering off of WIMPs by the bulge stars may be consistent with the observed gamma-ray flux if the spike forms on an underlying generalized NFW profile with γ ≲ 0.8, corresponding to a spike slope of γsp-star = 1.5 and spike radius of R sp-stars ∼ 25 30 pc. Instead, in the extreme case of the instantaneous growth of the black hole, the underlying profile could have up to γ ∼ 1.2, a corresponding γsp-inst = 1.4 and R sp-inst ∼ 15–25 pc. Finally, the results of our analysis of the total DM mass enclosed within the S2 orbit (updated with new GRAVITY data) are less stringent than the spectral analysis. Our work aims to guide future studies of the GC region, with both current and next generation of telescopes. In particular, the next Cherenkov Telescope Array will be able to scan the GC region with improved flux sensitivity and angular resolution.

Energy Stereoscopic System (H.E.S.S.) in 5 regions of the Galactic Center (GC).Motivated by previous works on a possible Dark Matter (DM) explanation for the TeV cut-off observed by H.E.S.S. in the innermost ∼ 15 pc of the Galaxy, we aim to constrain the DM distribution up to a radius of ∼ 450 pc from the GC.In this region, the benchmark approach (e.g.cosmological simulations and Galactic dynamics studies) fails to produce a strong prediction of the DM profile.Within our proof-of-concept analysis, we use DRAGON to model the diffuse background emission and determine upper limits on the density distribution of thermal multi-TeV Weakly Interactive Massive Particles (WIMPs), compatible with the observed gamma-ray flux.The results are in agreement with the hypothesis of an enhancement of the DM density in the GC with respect to the benchmark Navarro-Frenk-White (NFW) profile (γ = 1) and allow us to exclude profiles with an inner slope cuspier than γ ≳ 1.3.We also investigate the possibility that such an enhancement could be related to the existence of a DM spike associated with the supermassive black hole Sgr A* at the GC.We find out that the existence of an adiabatic DM spike smoothed by the scattering off of WIMPs by the bulge stars may be consistent with the observed gamma-ray flux if the spike forms on an underlying generalized NFW profile with γ ≲ 0.8, corresponding to a spike slope of γ sp−star = 1.5 and spike radius of R sp-stars ∼ 25-30 pc.Instead, in the extreme case of the instantaneous growth of the black hole, the underlying profile could have up to γ ∼ 1.2, a corresponding γ sp−inst = 1.4 and R sp-inst ∼ 15-25 pc.Finally, the results of our analysis of the total DM mass enclosed within the S2 orbit (updated with new GRAVITY data) are less stringent than the spectral analysis.Our work aims to guide future studies of the GC region, with both current and next generation of telescopes.In particular, the next Cherenkov Telescope Array will be able to scan the GC region with improved flux sensitivity and angular resolution.

Introduction
Many observations, such as galaxy rotation curves [1], gravitational lensing [2] and the cosmic microwave background [3] have led to the conclusion that non-baryonic Dark Matter (DM) constitutes about the 84% of the total mass content of the universe, yet its nature is still unknown.Among other candidates [4], cold DM is able to explain most of the astrophysical and cosmological evidence, with one natural candidate being the Weakly Interactive Massive Particles (WIMPs).Many WIMP candidates are expected to have been produced thermally in the early Universe, similarly to the particles of the Standard Model (SM), and usually constitute cold DM.Obtaining the correct abundance of DM today via thermal production requires a self-annihilation cross-section between ⟨σv⟩ ≃ 5.2×10 −26 cm 3 s −1 at ≈ 0.3 GeV and ⟨σv⟩ ≃ 2.2 × 10 −26 cm 3 s −1 above ≈ 10 GeV [5].Many efforts have been addressed to detect WIMPs, being focused mainly on searches at colliders [6], direct detection experiments [7] and indirect searches [8,9].SM particles are accelerated in colliders, where they could produce, among other particles, DM particles which can be detected as a "missing energy" in the original process.For the direct detection experiments, the transferred energy from the DM particle to SM particles via elastic collisions is expected to be detectable.In this work, we focus on indirect searches, which rely on the detection of secondary fluxes of astroparticles produced in the annihilation/decay of DM in astrophysical targets.Those fluxes can be observed in experiments and observatories such as High Energy Spectroscopic System (H.E.S.S.), Major Atmospheric Gamma-Ray Imaging Cherenkov (MAGIC), High-Altitude Water Cherenkov Observatory (HAWC), Fermi, etc., allowing to set constraints on a broad range of WIMP masses and the associated parameter space (annihilation cross-section or decay time).
Astrophysical targets of interest for indirect detection of DM are traditionally dwarf galaxies [10][11][12][13][14], the Galactic Center (GC) [15][16][17][18][19][20], or Galaxy Clusters [21,22].The most important features of an astrophysical target are its DM component, its distance, and all the possible uncertainties relating to the modeling of the object, like the background flux.Dwarf galaxies are desirable targets since they are DM-dominated systems, with a DM mass of 10 7 -10 10 M ⊙ and a negligible astrophysical background, that are close to the earth (∼ 0.5 Mpc).Focusing on specific targets, dwarf spheroidal galaxies have the counterpart of its dynamics being pressure-dominated, making it difficult to estimate its DM density profile.Dwarf irregular galaxies are, on the other hand, more distant objects but dominated by rotation, alleviating the problem (e.g., see [10,11]).Galaxy Clusters are DM-dominated as well, although they are distant objects with high astrophysical backgrounds, so propagation effects need to be modeled.Finally, the GC is the closest source to the Earth with the highest expected DM annihilation flux.However, the GC has a rich astrophysical background, both from sources and diffuse galactic emission, whose modeling can be a challenging task.
Besides the background modeling, one of the highest sources of uncertainty in the indirect detection of DM is modeling the DM density distribution profile in the astrophysical targets.The DM distribution has been a subject of debate for years.On one hand, DM-only numerical simulations favor steeper profiles, generally well described by the benchmark Navarro-Frenk-White profile (NFW, [23]).Hydrodynamical simulations, which also include the effect of baryons [24], seem to favor a generalized NFW profile with a different slope contracting the inner halo with higher densities [25][26][27] or with the existence of the bar in the Galaxy that could also affect the shape of the halo [27].On the other hand, observations of the rotation curves of dwarf galaxies prefer cored profiles [11,12].Regarding our Galaxy, the outer shape of the profile is mostly constrained by the observed rotation curves [1], whose results put constraints on the parameters that define it.Nonetheless, difficulties in obtaining dynamical measurements of the baryonic matter (due to the fact that we are inside the Milky Way) translate into difficulties when determining the distribution of the DM at a distance ≲ 2.5-3 kpc from the GC [28,29].Despite this issue, upper limits on the DM mass distribution in the very central regions of the GC (< 10 −2 pc) can be obtained by observing the orbits of the S stars surrounding the Super Massive Black Hole (SMBH) Sgr A* [30][31][32][33].When considering parsec scales, the uncertainty on the DM distribution is even worse.In fact, a possible enhancement of the DM density (hereafter, DM spike) could be associated with the growth of the central SMBH Sgr A*.This spike may have different characteristics depending on the evolution history of the SMBH, e.g. the adiabatic growth of the SMBH [34,35], the interaction with the stars surrounding the central spike [36,37], the extreme case of the instantaneous growth of the SMBH [38] or the effects of a rotating Kerr BH [39].The biggest question regarding DM spikes is their stability and survival in galaxies such as the Milky Way.
In this work, we aim to set new constraints on the DM distribution in our Galaxy.We develop a comprehensive study of compatibility between the analysis of the gamma-ray spectra observed by H.E.S.S. in 5 different concentric regions in the GC [19,[40][41][42][43][44][45], the dynamical constraints both in the outer [1,46] and inner region [33,47] of the Galaxy, and the possibility of having a detected multi-TeV DM signature in the inner GC [17,18,48].Despite the benchmark approach in DM indirect searches, which aims to set constraints on the DM annihilation cross-section by selecting a couple of possible (typically, one core and one cusp) DM density distribution profiles (see e.g.[49,50]), our main focus is to set constraints on the DM density profile within in a radius of 450 pc from the GC, by assuming a thermal DM candidate which could explain the gamma-ray spectra detected by H.E.S.S. in the inner 15 pc of the Galaxy [17,18].Our hypothesis is based on the fact that when considering a multi-TeV DM candidate with a thermal cross-section ⟨σv⟩ = 3 × 10 −26 cm 3 s −1 , a boost factor is needed in the DM density profile (assumed to be a benchmark NFW) to explain the gamma-ray flux observed by H.E.S.S. in the source J1745-290 at the Very Inner Region (VIR) of the Galaxy [17,18,51].Furthermore, the angular dimension of such local enhancement of the DM density could favor the disentangling between the existence of a cusp profile or a DM spike.This can be summarized as, from a morphological point of view, that the HESS J1745-290 data would be compatible with the scenario of a continuous DM annihilation signal.
The manuscript is organized as follows.In Section 2 we introduce the gamma-ray spectra observed by H.E.S.S. in 5 concentric regions at the GC and its possible interpretation as DM signal.In Section 3 we discuss the new background modeling, the spectral analysis and the results obtained for the spectral constraints on the DM distribution.We compare these results with DM density distribution profiles obtained by both dynamical studies and simulations in Section 4. We verify the compatibility of our results with the dynamical constraints obtained at sub-parsec scale by the analysis of the S2 star orbit in Section 5. Finally, the conclusions can be found in Section 6.

The Galactic Center region
As anticipated in the introduction, the GC region is still a hot topic for the indirect detection of DM.Due to the complexity of the region, distinguishing the DM signal and the astrophysical background is a very hard issue [52].So far, there are two main research directions dealing with a prospective DM signal at the GC: the first one is the GeV excess detected by Fermi-LAT [15,[53][54][55], the second one is the TeV cut-off by H.E.S.S. [9,17,18,[56][57][58].While  the first hypothesis has been extensively studied in the literature, we focus our study on the multi-TeV DM candidate for the gamma-ray cut-off detected by H.E.S.S.In this section, we introduce the collection of H.E.S.S. gamma-ray data in the GC region, the indirect detection of WIMP particles, the multi-TeV WIMP candidate and the open issue about the DM density distribution profile.

H.E.S.S. spectral data
Thanks to its position in the southern hemisphere, the H.E.S.S. telescope is currently the only Imaging Air Cherenkov Telescope capable of observing the GC at TeV energy with good angular and energy resolution.For this reason, we collect the H.E.S.S. published data of 5 concentric regions of the GC, with the aim to perform a spatial and spectral study of the gamma-ray flux in the region.The collection of 5 regions analyzed in this work is shown in Figure 1: the Very Inner Region (VIR) [40,48], the Ridge [42,43], the Diffuse emission region [40], the Halo [44,45] and the Inner Galaxy Survey (IGS) [19].General spectral and spatial information of each region is given by the H.E.S.S. collaboration, by fitting the gamma-ray flux with the following function: ) . (2.1) These values are summarized in the following lines: H.E.S.S. J1745-290 VIR.The very inner part of the GC region θ < 0.1 • (Figure 1, green region), with ∆Ω = 9.57×10 −5 sr.A power-law with Φ VIR = 2.55±0.41×10−12 TeV −1 cm −2 s −1 , Γ point = 2.14 ± 0.12 and a preference for an exponential cut-off at E cut = 10.7 ± 4.1 TeV is fitted [40].
H.E.S.S. Ridge.The region |b| < 0.  [19].For the Halo and IGS regions, in the ON-OFF analyses no excess is found.Instead, a Test Statistics analysis is performed to compute the corresponding upper limits.
In other words, a cut-off in the gamma-ray flux (compatible with a TeV DM annihilation signal, [40]) has been detected by the H.E.S.S. collaboration only in the VIR.From a morphological point of view, we hypothesize that it would be compatible with the scenario of a continuous DM annihilation signal.If a DM cusp exists at the GC, we may expect that in the inner region the annihilation component of the gamma-ray flux becomes more important over the background than in the external regions, where the DM detection could be locally suppressed by the astrophysical background.In the following sections, we study in-depth the gamma-ray flux in each region in order to model the DM density distribution in the Galaxy of the possible multi-TeV WIMP candidate, which well fits the gamma-ray cut-off in the VIR.

Gamma-ray flux from WIMP annihilation
The gamma-ray flux produced by WIMP annihilation has the following form: Where ⟨σv⟩ i is the thermally averaged annihilation cross-section for the WIMP particle, m DM is the mass of the DM particle, and ∆Ω is the solid angle of the observed region.In a model-independent approach, the WIMP particle is considered to annihilate in only one i-th channel of the Standard Model (SM).Also, by assuming a specific WIMP candidate with different branching ratios, the annihilation could happen in a combination of several SM channels.In particular, brane-world DM is an interesting possibility to naturally produce thermal multi-TeV WIMP candidates [59][60][61][62].Both those possibilities have been studied in [17,18,57,61].Motivated by those previous studies, we focus on the ZZ channel, which gives one of the best fits to the H.E.S.S. data in the inner 15 pc of the GC.The factor dN i dE is the differential flux of secondary particles, in this case gamma rays, produced by the possible hadronizations, decays, annihilations and other interactions that can create them starting from the primary products of DM annihilation, indeed the SM annihilation channel.They are computed using PPPC4DMID [63,64].It can be shown that including interactions on the secondary particles created by DM annihilation such as electrons/positrons (Inverse Compton, Bremsstrahlung, and synchrotron radiation) the total DM spectra can change, but since we are in an energy range of ∼ 10 2 -10 5 GeV this flux is a few orders of magnitude below the primary gamma-ray flux [65,66], so it is a good approximation to neglect them.Note that, for kinematic reasons, the energy of the gamma rays generated by DM annihilation cannot exceed its own mass.
Finally, the J-factor ⟨J⟩ ∆Ω (also called astrophysical factor in the literature) is where the information about DM distribution profile is contained [67]: 3) The J-factor is the integral of the DM density profile squared ρ DM (r), along the line of sight l, and averaged over the solid angle ∆Ω.In these coordinates, l( θ) min and l( θ) max are the edges of the regions observed in the direction given by θ (if the region corresponds to the integration starting at the position of the Sun, l( θ) min = 0).r is the radial coordinate from the center of the GC (taken in Sgr A*), and can be related to the line of sight l with the expression , where R ⊙ = 8.277 kpc is the distance from the Sun to the Sgr A* ( [30], hereafter GRAVITY2021).Assuming a virial radius R vir for the DM halo, we have that l( θ

The multi-TeV WIMP candidate
According to [17,18], we assume here that the total gamma-ray flux detected by H.E.S.S. in the GC region is a combination of the gamma-ray flux produced by self-annihilating DM particles dΦ DM dE (Equation 2.2) and a background component This hypothesis is well motivated by the difficulty with modeling the astrophysical background emission in the region, due to the abundance of known and unknown astrophysical sources because of the particularly bright diffuse emission caused by the cosmic-ray interactions with the interstellar medium.We allow the normalization of the background to vary for matching the data and allow significant fluctuations of this normalization in the very small regions (such as the VIR) with respect to the Ridge region, motivated by a combination of uncertainties on modeling the column density of the molecular gas and on the properties of the CR sea at small scales.
By assuming a benchmark NFW profile, a boost factor of ∼ 10 3 in the J-factor was needed in order to explain the observed signal in the inner 15 pc of the GC, with a DM mass m DM ≃ 50 TeV [17,18].On the one hand, such a boost factor may be compatible with the characteristics expected by an enhancement of the DM distribution in the region due to the The parameters of the profiles can be found in Table 1.
presence of the SMBH Sgr A* [51].Depending on the DM distribution profile in the Galaxy (cusp or core), both the slope and the radius of the DM spike change, by giving a footprint on the inner DM distribution and, indeed, on the morphology of the expected gamma-ray emission due to the annihilation of WIMPs [51].On the other hand, the spectral cut-off strictly depends on the DM mass and the background component.To derive these results, the authors used a background given by a simple power law ( [17,18].Our work is based on these results but with the novelty of using the state-of-the-art background model computed with the public numerical codes DRAGON and HERMES [68,69] and extending the analysis up to 5 concentric regions in the GC.

The DM density profile
As introduced in Section 2.2, high uncertainty exists in determining the DM density distribution profiles in the Galaxy.Of course, this fact has a high impact on studying prospective DM annihilation signals.In this work, we consider a wide range of DM density profile models, from cuspy to cored profiles and from both cosmological simulations and dynamical observations.Here, we adopt the Navarro-Frenk-White (NFW) generalized formalism [70]: where ρ s is a normalizing factor and r s is the scale radius of the profile.This expression allows us to recover the benchmark NFW profile resulting from DM-only N-body simulations, Profile  2.5 for the DM-only simulation, the hydrodynamical simulations GARR, GARR-I, GARR-I300, GARR-II300, ERIS, MOLL and EAGLE (see [51] and the references within) and the observational models McMillan17 [46] and Benito20 [1].Also, it is shown the local DM density ρ⊙ of each simulation, with R⊙ = 8.277 kpc (GRAVITY2021 [30]).
Under this formalism, we can also describe two models obtained from observations: 1) the set of NFW-like models for values of γ between 0 and 1.5 [46] (hereafter, McMillan17); 2) the 2σ constraints on the DM profile in the Milky Way obtained by the most recent analysis of the galactic rotation curve [1] (hereafter, Benito201 ).For both the observational models, we have α = 1 and β = 3.The parameters of each DM profile are given in Table 1 and the DM profiles themselves are shown in Figure 2. From this figure, it is easily visible that -due to the high uncertainty in the data set -the range of parameters obtained for the DM profiles via the study of the galactic kinematic cover all the different profiles obtained via different cosmological simulations.For a better comparison of the density distribution profile parameters, see Section A in the Appendix.
Let us remark that neither simulations nor the study of the Milky Way rotation curve allows us to set constraints on the DM profile at parsec scales.In this work, we aim to set interesting upper limits on the DM density distribution in such an unconstrained region via the comprehensive analysis of the gamma-ray spectra detected by H.E.S.S. in the previously introduced 5 regions of the inner Galaxy (within 450 pc), by following the hypothesis of the existence of a DM component.Interestingly, this approach allows us to also set constraints not only on benchmark outer DM profiles (i.e. a core or cusp profile) but also on more complex models, which include possible modifications in the inner part of the outer benchmark DM profile, such as: 1) the existence of a DM-only adiabatic spike [34,35], 2) the effect of the interactions between stars and an adiabatic spike [36,37], 3) the extreme case of an instant spike [38], and 4) the perturbation due to a rotating Kerr SMBH [39].Furthermore, the matter distribution in the Galaxy below ≲ 10 −2 pc can be constrained by dynamical constraints of the S stars [33].This fact will represent an additional independent cross-check to the prospective DM distribution modeled by the analysis of the gamma-ray spectra, as explained in Section 5.

Spectral analysis
Our analysis of the spectral data is based on a χ 2 statistical approach: where E 2 i y HESS i and E 2 i y error i are, respectively, the differential gamma-ray flux observed by H.E.S.S. and its uncertainty within the E i energy bin; E 2 i dΦ total dE (E i ) is the differential flux computed following Equation 2.4 and N is the total amount of data points observed in the region.Therefore, the ddof of this analysis would be the number of data points N minus the number of the free parameters fitted.

Background model
Beyond the uncertainty in modeling the DM density distribution profile, the second main issue in studying the GC region is the background flux modeling.Due to the presence of a very bright diffuse emission and several point sources, plus possibly a population of unresolved sources, separating the prospective DM annihilation signal from the astrophysical background emission is a hard task.
The diffuse background is due to the presence of a diffuse charged cosmic rays (CR) "sea" in the Galaxy, effectively confined by its turbulent magnetic field [75][76][77].The presence of this emission over a wide range of energies constitutes an unavoidable challenge whenever one tries to pinpoint any emission of exotic origin [78].Charged CRs interact with several components of the interstellar medium (interstellar gas, low-energy photons, regular and turbulent magnetic field) and emit a bright, diffuse radiation from the radio domain all the way up to the multi-TeV gamma-ray band as a consequence of a variety of non-thermal processes (mainly synchrotron, bremsstrahlung, neutral pion decay and Inverse Compton scattering).We model such background emission by means of the Diffusion Reacceleration and Advection of Galactic cosmic rays: an Open New code (DRAGON) [68,69,79,80], developed within the HERMES publicly available computational framework for the line of sight integration over galactic radiative processes, which creates sky maps in the HEALPix-compatible format [81].The diffuse model we adopt here is described in detail in a recent paper [82] and is tuned on different sets of local charged cosmic-ray, multi-wavelength and gamma-ray data along the Galactic plane.It features a harder diffusion coefficient (hence, a harder propagated proton spectrum) compared to the locally measured one, as hinted by Fermi-LAT data and first introduced in [83].This feature was shown to be consistent with the bulk of the emission observed by the H.E.S.S. collaboration in the Galactic Ridge region [52].Indeed, given this choice for the background, in Equation 2.4 we use where dΦ DRAGON dE is obtained for each of the 5 regions by applying to the DRAGON sky map the same mask as in the H.E.S.S. data, and B is an O(1) normalizing factor that is left as a free parameter.The reason why the overall normalization is allowed to vary is the existence of an intrinsic, unavoidable uncertainty in the model itself, which ultimately stems from our poor knowledge of the conversion factor between the CO emissivity and the molecular gas column density (see e.g. the recent review [84]).The diffuse emission scales linearly with this poorly constrained parameter which we leave free to vary together with the exotic component originating from DM annihilation.

Methodology
We assume a thermal WIMP particle with the benchmark annihilation cross-section ⟨σv⟩ ≃ 2.2 × 10 −26 cm 3 s −1 necessary to explain the observed DM relic abundance in the universe [5].Indeed we fit the following equation to the data of each region: To perform the fits, the expected annihilation gamma-ray signal dN/dE given by [63,64] is convoluted to a Gaussian distribution (with σ/E = 0.10) to take into account the energy resolution of the H.E.S.S. telescope [49].

VIR
In Figure 3 (left panel) we show the fit of the DM plus DRAGON background component (Equation 2.4) to the H.E.S.S. data in the VIR.Following previous analyses [17,18], the inclusion of a background component in the analysis of the VIR region is required in order to well fit the low energy H.E.S.S. data, and in agreement with high energy Fermi-LAT data, whose angular resolution above ∼ 10 GeV is θ ∼ 0.1 • , comparable to H.E.S.S..Although we do not include Fermi-LAT data in our analysis our background model DRAGON and renormalization factor B 2 are in agreement with the gamma-ray spectra detected by Fermi-LAT at energy E > 10 GeV (see in the Conclusions Figure 12, upper left panel).The fitted parameters, i.e. the DM mass m DM , background renormalization B and astrophysical J-factor ⟨J⟩ ∆Ω are shown in Table 2: the good quality of the fits show an agreement of the data with the hypothesis of a DM component to the gamma-ray signal, produced by a m DM = 36 +4 −6 TeV DM particle annihilating in the ZZ channel at GC, with a J-factor ⟨J⟩ ∆Ω = 2.7 +0.6 −0.3 × 10 28 GeV 2 cm −5 .This result is compatible with previous works [17,18].Minor differences may be  2. Right Panel: 2σ confidence level region (in grey) for the J-factor ⟨J⟩∆Ω and B 2 parameter space.The red dot is the best fit of the χ 2 analysis, with χ 2 min /ddof = 0.89.Note that the 2σ region is compatible with ⟨J⟩∆Ω = 0, indeed no DM signal.The best fit for the Ridge region (black cross) is compatible with the Diffuse best fit within the 2σ confidence level.
⟨J⟩ ∆Ω (GeV 2 cm −5 ) 2.7 +1.0 −0.6 × 10 28 2.5 +1.0 −0.9 × 10 27 1.1   3 and 4).The DM mass mDM corresponds to the best fit in the VIR and it has been kept fixed in the rest of the regions.The B parameter is the normalization of the DRAGON background model, and ⟨J⟩∆Ω is the J-factor fitted.We also show the boost factor to the reference value for an NFW profile, i.e.J VIR DM-only = 2.6 × 10 25 GeV 2 cm −5 , J Ridge DM-only = 2.5 × 10 24 GeV 2 cm −5 and J Diff DM-only = 4.2 × 10 24 GeV 2 cm −5 .Finally, the χ 2 /ddof of each fit is shown with the solid angle ∆Ω of each region.The uncertainties of the values correspond to a 2σ confidence level.
explained by: 1) newly updated data by the H.E.S.S. collaboration in the region [40,48]; 2) the inclusion of the DRAGON background model instead of the simple power-law [17,18]; 3) the use of a different Monte Carlo event generator software to model the gamma-ray annihilation flux [85].In Table 2 we also give the ⟨J⟩ ∆Ω /J DM-only es reference value: this value represents the enhancement factor (when compared to a benchmark NFW profile) required to the DM density distribution to be the origin of such a gamma-ray flux.This boost factor represents the first hint to assume the existence of a more cuspy profile or a DM spike, as we discuss in the next Section 4. Following this idea, it is straightforward to expect that the DM signal could manifest only locally over the background, remaining undetected in outer regions.Indeed, we keep the best fit of the DM mass obtained in this region to develop the fits in the Ridge, Diffuse, Halo and IGS, where the DM component is expected to be subdominant.

Ridge
In Figure 3 (right panel) we show the data observed by H.E.S.S. in the Ridge, fitted to the Equation 2.4.Here, we fixed the DM mass to the best-fit value obtained in the analysis of the VIR region, and we fit here only two free parameters: the background renormalization B and the astrophysical J-factor ⟨J⟩ ∆Ω .The best-fit parameters are given in Table 2.We have also verified that without fixing the m DM to the VIR value, we get a similar result with higher uncertainty.

Diffuse
As for the VIR and the Ridge regions, we perform a χ 2 analysis (Equation 3.1) of the Diffuse region data.In Figure 4 (right panel) we show the best-fit parameters (red dot) and the 2σ confidence level (the grey region in the same figure).The best-fit parameters are also shown in Table 2.The 2σ region is compatible with ⟨J⟩ ∆Ω = 0, meaning that the DM gamma-ray component is subdominant to the background signal in this region, and the result could be considered as an upper limit to the J-factor.This result is also compatible with the hypothesis that most of the DM-related gamma-ray signature is hidden in this region, being the contribution of the VIR covered by the applied masked.

Halo and IGS
In the Halo [44,45] and IGS [19] regions, no signal is detected by the H.E.S.S. collaboration.In both regions, the flux in the ON region is compatible with the flux in the OFF control region, adopted for background rejection.Indeed, in this case, we adopt a different procedure in order to set upper limits on the DM density distribution and astrophysical J-factor, and we ask for a 2σ detection of the theoretical gamma-ray flux, i.e.Φ th = Φ DM + Φ DRAGON .We determine the background given by DRAGON for the same regions and masks, and we scan the J-factor and B parameter space for a thermal WIMP mass of 36 TeV annihilating in the ZZ channel, consistently with the fit in the VIR, Ridge and Diffuse (see Figure 5 and Table 3).Indeed, we have [11,86]: UL (cm 2 s −1 ) 0.13 0.02 ⟨J⟩ UL ∆Ω (GeV 2 cm −5 ) 2.5 × 10 26 1.7 × 10 25 ∆Ω (sr) 5.97 × 10 −4 6.38 × 10 −3 Table 3. Upper limits within 2σ confidence level on the B 2 and J-factor parameters for the Halo and IGS regions, calculated in the respective solid angle ∆Ω.The value of the DM particle mass is fixed, in agreement with the VIR best-fit value, i.e. 36 TeV.The reported values are the maximum values allowed within the grey area in Figure 5.For comparison, J Halo DM-only = 1.8 × 10 24 GeV 2 cm −5 and J IGS DM-only = 3.6 × 10 21 GeV 2 cm −5 .
where A eff is the H.E.S.S. effective area, t exp is the exposure time, ∆Ω is the solid angle and Φ HESS is the integrated flux observed by H.E.S.S., compatible with the background component.

Results
So far, we have determined the upper limits on the J-factor for the prospective thermal multi-TeV DM candidate of m DM = 36 +4 −6 TeV in 5 concentric regions of the GC observed by H.E.S.S.These J-factors have been determined as the amplitude required for the DM component to fit the H.E.S.S. data.It is well known that the J-factor is, by definition, the integral of the squared DM density distribution profile along the line of sight (Equation 2.3).Indeed, within our hypothesis, we can set upper limits on the radial DM density distribution profile in the Galaxy by the study of the gamma-ray spectra.
In Figure 6 we show the J-factors obtained by the analysis of the gamma-ray spectra.We compare our results with the J-factors obtained for the same regions by assuming the benchmark DM-only NFW profile (orange dashed lines) and the Benito20 model (brown region).From the analysis of the VIR, we find out J-factors higher than the benchmark NFW profile, in agreement with the range of J-factors obtained by integrating the DM profile obtained in Benito20, extrapolated up to those regions.In particular, a boost factor of ∼ 1000 over the benchmark NFW profile is required to explain the fitted values.In the Diffuse and Ridge regions, the best fit of the gamma-ray spectra favors a DM distribution profile steeper than the extrapolation of the Benito20's profile across the GC.Note that, for the Diffuse, as explained in the right panel of Figure 4, the lower limit of the 2σ confidence level is compatible with ⟨J⟩ ∆Ω = 0 and, therefore, the results are shown as upper limits.Finally, in the Halo and IGS outer regions the spectral upper limits are not that restrictive.
It is worth mentioning that the dynamical data studied in Benito20 is able to give constraints on the DM distribution only to r ≳ 2.5 kpc from the GC.Indeed, our study represents a complementary analysis to solve the issue of the DM distribution in the Galaxy.Our results are the first upper limits on the DM distribution for the innermost Galactic region within r ≲ 400 pc, in agreement with the previous constraints on the DM distribution in the outer region of the Galaxy.In the following section, we study the possibility to disentangle within our study a cuspy DM profile from a DM spike.In fact, the enhancement of the DM density distribution in the VIR (when compared to the benchmark NFW profile) could be GC, where r is the distance from the GC in pc and b represents the galactic latitude (deg).Each region is represented by the vertical grey dotted lines, corresponding to their boundaries.For the VIR and Ridge regions, the best fit and the 1σ and 2σ uncertainties are represented; for the Diffuse region, the best fit and the 2σ uncertainty are shown, where the lower limit of the 2σ confidence level is compatible with ⟨J⟩∆Ω = 0, as explained in the right panel of Figure 4. Finally, for the Halo and IGS regions, the values correspond to the 2σ upper limit.For comparison, we show the J-factors computed with the DM-only profile (orange dashed lines) and with the Benito20 DM profile [1] (brown area).Here, the DM mass is mDM = 36 +4 −6 TeV and we assume a thermal annihilation cross section ⟨σv⟩ ≃ 2.2 × 10 −26 cm 3 s −1 .Due to the non-spherical symmetry of the Ridge region, we show with the dot-dashed lines the extension of the region in the galactic longitude, up to 1 • .explained by the interaction of the WIMP particles with both the gravitational potential of the SMBH Sgr A* and the baryonic component of the region.

The innermost dark matter density profile
In this section, we compare our results with several DM density distribution profiles, obtained by both cosmological simulations and studies of Galactic dynamics.First, we consider the generalized NFW profile (Equation 2.5), by extrapolating the profiles to the innermost GC region without any modifications to the inner part.Second, we investigate the possibility that the DM profile could be modified in the central region due to the presence of the SMBH Sgr A*.We consider different models: 1) a DM spike formed adiabatically [34,35]; 2) the effect of the interaction with the stars of the bulge with the adiabatic spike [36,37]; 3) a DM spike formed instantaneously [38]; and 4) we review the case for the rotating Kerr BH.As an example, in Figure 7 we show all those models applied to the benchmark NFW profile.Furthermore, it is well known that the DM density distribution in the innermost GC region cannot increase indefinitely [34].Instead, in the WIMP framework, it is typically expected to reach a maximum plateau due to the growth of the annihilation rate with the TeV, with their uncertainties represented in grey: 1σ and 2σ for the VIR and Ridge, 2σ for the Diffuse and the Halo 2σ upper limit.The colored dots correspond to DM density models from N-body simulations and the colored regions to the observational models: McMillan17 [46] in blue and Benito20 [1] in brown.The vertical gray dashed line indicates when the integrated J-factors cross the fitted values.Above them, the predicted gamma-ray flux exceeds the observed flux by H.E.S.S. and, therefore, this part of the parameter space is excluded.first analysis we can rule out any DM distribution profile which could produce a gamma-ray flux above the H.E.S.S. data.Indeed, we can exclude generalized NFW profiles with γ ≳ 1.3.

Adiabatic Dark matter Spike
The existence of the SMBH Sgr A* in the GC with a mass M BH = 4.297 ± 0.012 × 10 6 M ⊙ [30] has been recently demonstrated [87].The possibility to have an enhancement of the DM density due to the gravitational interaction with the SMBH Sgr A* has been largely studied in the literature.Among other models, a DM spike is considered to be formed adiabatically if the SMBH has grown in the center of the Galaxy, very slowly compared to the typical dynamical timescales of the region, and without any occurrence of big mergers in the Milky Way for the last ∼ 10 Gyr [88] (a conservative assumed age for the SMBH Sgr A* [33]).for different DM density profiles: DM-only simulation, the Hydrodynamical simulations GARR, GARR-I, GARR-I300, GARR-II300, ERIS, MOLL and EAGLE (see [51] and the references within) and the range of values given by the observational models McMillan17 [46] and Benito20 [1].* Larger Rinst corresponds to smaller γ.
The DM adiabatic spike density profile [34,35] can be described with the modification of the inner part as ρ sp (r) ≃ ρ R (1 − 2Rs r ) 3 ( r Rsp ) −γsp .In this description, both the slope and the radius of the spike depend on the inner slope γ of the underlying DM density profile: the spike slope is γ sp = (9 − 2γ)/(4 − γ) and the spike extends up to radius ; where α γ is computed numerically and ρ R = ρ halo (R sp ) [34].In Table 4 we show the γ sp and R sp for all the DM profiles considered in this work.Regarding DM density spikes with DM self-annihilating particles, it is important to stress that the DM density can only reach a maximum value ρ sat .This modifies the total DM density profile as follows: where R S is the Schwarzschild radius and ρ halo (r) is the underlying DM density distribution profile.In Figure 7 we show the adiabatic spike (dashed line) associated with the benchmark NFW profile (solid line), as an example.In the literature, the saturation factor ρ sat is usually taken as a central plateau ρ sat = m DM /(t BH ⟨σv⟩), where t BH is the age of the central SMBH, with a value ∼ 10 10 yr [33,51].Here, we relax this plateau with the expression ρ sat (r) = m DM /(t BH ⟨σv⟩)( r Rsat ) −1/2 , where this −1/2 slope appears as a correction when taking into account anisotropies in the velocity distribution of the DM particles in the GC [89], being R sat the radius at which ρ(R sat ) = ρ sat .Furthermore, we follow the relativistic approach [35] and assume that the DM density distribution vanishes at r < 2R S instead of r < 4R S as in the Newtonian analysis [34].
In Figure 8 (lower row), we compare the J-factor obtained by the spectral analysis with the J-factor expected in the case of having an adiabatic DM spike on different underlying DM profiles.We find out that, by analyzing the VIR region (first panel), we can exclude all profiles except very shallow ones like MOLL.Nonetheless, the latter enhancement would be indistinguishable in the gamma-ray spectra.Indeed, it would be compatible only with an astrophysical origin of the VIR spectral feature.The Ridge, Diffuse and Halo regions (rest of the panels) do not allow us to exclude any profile.In fact, all those 3 regions do not include the inner region inside θ ≲ 0.15 • (Figure 1).A negligible contribution of the prospective DM spike to the Diffuse region is in agreement with a DM spike angular extension smaller or compatible with the inner and outer radius observed by H.E.S.S., i.e. θ sp ≲ θ Diff = 0.15 • -0.45 • or R sp ≲ R Diff = 22-65 pc, also in agreement with the result obtained by the spectral analysis in Section 3.2.

Effect of stars on the spike
The previous approach only takes into account the consequences of the inner DM density profile when considering the adiabatic growth of the SMBH, without considering any other interactions with other elements of the galaxy.As a correction, we can study the same case with the addition of the interactions with the stars of the bulge [36,37].In particular, with this approach that we name stars-spike in this paper, we are including the scattering of DM particles off of stars, which causes the dynamical heating of DM particles and, consequently, flattens the existing spike.Another process is the DM capture by the stars, but this process is secondary unless the cross-section of WIMP-baryons is sufficiently large.In summary, when taking into account the effect of the stars, inner profiles are less cuspy than the adiabatic spike case [36,37], lowering the internal slope to γ star = 1.5 inside the radius of influence of Sgr A*, i.e. r b = 2 pc.In Figure 7 we show the change of the profile (dotted line) with respect to the adiabatic case (dashed line).
As a consequence, the values of the J-factors are also generally smaller (upper row of Figure 9).In this case, for the VIR region, we can see in the figure that we can explain the fitted values with a stars-spike, formed on an underlying generalized NFW profile, with a slope γ ∼ 0.8.As a consequence, we can rule out any DM distribution profile that produces gamma-ray flux above the H.E.S.S. data, with a slope γ ≳ 0.8.When studying the Ridge, Diffuse and Halo regions no profile can be ruled out since the effect of the spike is irrelevant in those regions, meaning that the angular extension for the stars-spike cannot extend in the Diffuse and Ridge region, θ star ≲ θ Diff = 0.15 • -0.45 • or R star ≲ R Diff = 22-65 pc, sharing the same spatial constraint as in the adiabatic case.
Actually, the star-spike scenario is an appealing one.In fact, the possibility of the existence of a DM density enhancement in a form of an adiabatic spike in the Milky Way remains unclear, mostly due to the uncertainty on its stability and the evolution history of the Galaxy.However, the small dimensions of the possible spike, R sp ≲ 22-65 pc or θ ≲ 0.15 • -0.45 • , make it very difficult to prove the existence of the DM spike in our galaxy or even in other galaxies so far.

Effect of the Kerr Black Hole
So far, we have considered in our analysis a non-rotating SMBH.Nonetheless, the recent observation of the Sgr A* shadow by the Event Horizon Telescope shows evidence to conclude that the SMBH at the center of our Galaxy is a rotating Kerr BH [87,90].In fact, the shadow diameter depends on the metric and its properties, e.g. the BH spin.As a result, we investigate whether the SMBH spin could affect the detection of a DM signature.As an example, assuming an Hernquist profile ((α, β, γ) = (1, 4, 1)) an extra maximum factor of 2 appears in the J-factor -depending on the value of the spin-when compared to the non-rotating case.Nonetheless, this extra factor mostly disappears when considering the annihilation factor ρ sat [39].Also, Kerr BH shows asymmetric DM spikes, the asymmetry being dependent on the direction of rotation [39,91].However, this asymmetry is only important in the inner ∼ 10 −6 pc (or 5R S ), so it would not be visible by the current gamma-ray telescopes.Another important factor is that the annihilation products are sensitive to the SMBH spin: 1) in general higher spins generate higher-energy annihilation products and 2) the annihilation products from the DM particles close to the BH get redshifted, as this population moves closer and closer to it [91].Hence, because of the small differences between the rotating and non-rotating cases, for simplicity, we will focus only on the non-rotating case.

Instant spike
Finally, we consider the extreme case that the SMBH is formed instantly as a consequence of a non-adiabatic process [38] and the consequences for the DM density profile (dot-dashed line in Figure 7).In this case, the spike has a less steep profile with an extension of R inst and slope γ inst , where typically have values of ∼ 50 pc (∼ 0.15 • ) and slopes of ∼ 4/3.The form of the instant spike density profile is given by: Where γ inst and ρ inst are computed numerically following the prescription given in [38] for each profile, and R inst is taken such that ρ inst = ρ halo (R inst ).In Table 4 we show the values of the R inst and γ inst parameters for each profile considered in this work.
We remark that the "instantaneous" spike and the "adiabatic" spike represent the two relevant opposite limits that have been studied analytically; however, a realistic formation scenario of such DM overdensities may proceed through a non-trivial combination of adiabatic and non-adiabatic phases.A comprehensive description of this process under realistic assumptions is beyond the scope of the present work.
In Figure 9 (lower row) we show the results of the integrated J-factors compared to the fitted values.As we can see, a cuspy profile of γ ∼ 1.2 can explain the fitted J-factors in the VIR region.Note in the instant case only it has been represented γ up to 1.25, this is because at higher γ the instantaneous spike is very similar to the profile without modification.As we can see, profiles with γ ≳ 1.2 produce gamma-ray flux above the measured H.E.S.S. data in these regions, so we can rule them out.Studying the Ridge, Diffuse and Halo regions from the figure, we see similar results as the rest of the formalisms used (generalized NFW, adiabatic spike or instant spike) with no profile being ruled out.As a conclusion, the same spectral constraint can be extracted as before but now with the instantaneous case: θ inst ≲ θ Diff = 0.15 • -0.45 • .

S2 star dynamical constraints
In the previous sections, we have presented the constraints on the DM density profile obtained by the spectral study of a collection of gamma-ray observations of the GC by H.E.S.S. and the predictions obtained via both cosmological simulations and galactic dynamics of the outer region of the Galaxy.Nonetheless, we find out that different scenarios could explain the gamma-ray flux observed in the very inner region of the Galaxy, e.g. a more enhanced NFW profile with γ ∼ 1.3 or an adiabatic DM spike formed on a shallower profile.In order to disentangle this dichotomy, we performed another independent analysis with the study of the S2 star orbit [33,92].
We can set further constraints on the inner DM density distribution by considering the dynamical constraints obtained by the study of the inner S2 star precession: GRAVITY2020 [31], GRAVITY2021 [30] and Do+2019 [32].Indeed, the extended mass within the orbit can modify its predicted precession up to a point where constraints can be extracted.In particular, we update the results obtained in [33], by including new observational data (the GRAVITY and Do+2019 data). 2ndeed, we have computed the limits in the parameters (γ sp and R sp ) from the mass enclosed within the orbits of the stars, for a DM mass m DM = 36 TeV, in agreement with the best fit of the VIR.In Figure 10 we show the exclusion region for such a parameter space.In the left panel, we show the constraints integrating the McMillan17 profile with the BH adiabatic spike prescription given by [34,35] and computing the maximum R sp allowed by the independent and complementary to the gamma-ray analysis, although it is less constraining.On the one hand, beyond the adiabatic spike, a less spiky profile cannot be excluded from the study of S2 orbit.On the other hand, the adiabatic spike is already excluded from the cuspy-like profile by the spectral analysis of this work.
To conclude, our analysis seems to discard the possibility of having an adiabatic DM spike in the GC.Also, it reinforces the possibility of a DM profile with slope γ ∼ 1.3.The latter possibility should be further checked with an in-depth study of a region between 0.1 • −1 • deg, where the DM density is expected to increase gradually (Figure 12, lower right panel).Our work represents a proof-of-concept paper in order to guide the study of the GC with both current and next-generation of gamma-ray telescopes.In particular, the Cherenkov Telescope Array (CTA) [16] will be able to scan the GC region with a flux sensitivity of ∼ 1 order of magnitude better than H.E.S.S. and angular resolution of ∼ 0.05 • -0.03 • , instead of the ∼ 0.1 • of H.E.S.S. values given by numerical simulations (Table 1) and the blue and brown regions correspond to the McMillan17 [46] and Benito20 [1] models (1σ and 2σ regions), given by observations.

Figure 2 .
Figure2.DM density distribution profiles considered in this work.Profiles represented with a solid line correspond to cosmological simulations: DM-only (orange line), GARR-I (gray), GARR-I300 (light blue), GARR-II300 (red), ERIS (purple), MOLL (green) and EAGLE (blue).The blue (McMillan17) and brown (Benito20) regions cover the range of possible profiles which are compatible with the study of Galactic dynamics (see main text for details).To ensure an easy reading, we kept this color code along the manuscript.All the profiles are extrapolated up to the inner radius of the Galaxy, defined by the Schwarzschild radius (see Section 4.2 for further details).The grey horizontal dashed line represents the central plateau due to DM annihilation ρsat = mDM/(tBH⟨σv⟩) = 1.4 × 10 20 M⊙kpc −3 for a mass of mDM = 36 TeV (see text for further details).The parameters of the profiles can be found in Table1.

6 Figure 3 .
Figure 3. Best fit (magenta line) to the H.E.S.S. data of the VIR (left panel) and Ridge (right panel).The background (green-dashed line) and DM component (red-dashed line) are also shown.The best-fit parameters (Table2) show an agreement with a gamma-ray signal produced by the annihilation of a 36 +4 −6 TeV WIMP particle and a background component, taking into account the different background component as modeled by DRAGON in the different regions (see text for further details).The dark grey and light gray regions correspond to the 1σ and 2σ uncertainty, respectively.

Figure 4 .
Figure 4. Left Panel: same as Figure 3 for the Diffuse region.The light grey band is the 2σ confidence level.The best-fit parameters are shown in Table2.Right Panel: 2σ confidence level region (in grey) for the J-factor ⟨J⟩∆Ω and B 2 parameter space.The red dot is the best fit of the χ 2 analysis, with χ 2 min /ddof = 0.89.Note that the 2σ region is compatible with ⟨J⟩∆Ω = 0, indeed no DM signal.The best fit for the Ridge region (black cross) is compatible with the Diffuse best fit within the 2σ confidence level.

0B 2 (cm −2 s − 1 )Figure 5 .
Figure 5. 2σ upper limits on the J-Factor obtained by the analysis of the Halo region (left panel) and IGS region (right panel).The B 2 parameter corresponds to the renormalization of the DRAGON background model and the grey area corresponds to the 2σ upper limits.The upper limits are computed assuming the thermal relic cross-section ⟨σv⟩ = 2.2 × 10 −26 cm 3 s −1 and with a DM particle of mDM = 36 +4 −6 TeV, the best-fit value in the VIR, annihilating in the ZZ channel.

Figure 6 .
Figure 6.J-factor ⟨J⟩∆Ω obtained by the study of the gamma-ray spectra in 5 different regions of the

Figure 8 .
Figure 8. J-factor ⟨J⟩∆Ω as a function of the density slope γ, for the approaches: generalized NFW (first row) and DM adiabatic spike (second).We show, for each row the regions VIR (first panel from left to right), Ridge (second panel), Diffuse (third panel) and Halo (fourth panel).The horizontal lines are the fitted values from the gamma-ray spectra, assuming the thermal relic cross-section ⟨σv⟩ ≃ 2.2 × 10 −26 cm 3 s −1 and a fitted DM mass mDM = 36 +4 −6 TeV, with their uncertainties represented in grey: 1σ and 2σ for the VIR and Ridge, 2σ for the Diffuse and the Halo 2σ upper limit.The colored dots correspond to DM density models from N-body simulations and the colored regions to the observational models: McMillan17[46] in blue and Benito20[1] in brown.The vertical gray dashed line indicates when the integrated J-factors cross the fitted values.Above them, the predicted gamma-ray flux exceeds the observed flux by H.E.S.S. and, therefore, this part of the parameter space is excluded.

Figure 9 .
Figure 9. Same as Figure8, but now for the approaches stars-spike (first row) and instant spike (second row).

Figure 13 .
Figure 13.Parameters of the different DM density profiles used in this work.The dots correspond to the

Table 1 .
Parameters of different DM density profiles as in Equation

Table 2 .
Best-fit parameters for the VIR, Ridge and Diffuse regions (Figures

Table 4 .
Parameters used for the adiabatic spike (Equation 4.1) and instantaneous spike (Equation 4.2)