Constraining the Inner Galactic DM Density Profile with H.E.S.S.

Abstract

In this short review, corresponding to a talk given at the conference “Cosmology 2023 in Miramare”, we combine an analysis of five regions observed by H.E.S.S. in the Galactic Center, intending to constrain the Dark Matter (DM) density profile in a WIMP annihilation scenario. For the analysis, we include the state-of-the-art Galactic diffuse emission Gamma-optimized model computed with DRAGON and a wide range of DM density profiles from cored to cuspy profiles, including different kinds of DM spikes. Our results are able to constrain generalized NFW profiles with an inner slope $\gamma \gtrsim 1.3$. When considering DM spikes, the adiabatic spike is completely ruled out. However, smoother spikes given by the interactions with the bulge stars are compatible if $\gamma \lesssim 0.8$, with an internal slope of ${\gamma }_{\mathrm{sp-stars}}=1.5$.

1. Introduction

In terms of the matter composition of the Universe, about $84\%$ of its matter content is in the form of a non-baryonic matter, known as Dark Matter (DM). Independent observations such as the Cosmic Microwave Background [1], galaxy rotation curves [2], and gravitational lensing studies [3], among others, have led to the same conclusion. The DM theory that can explain more successfully those pieces of evidence is known as cold DM, characterized usually by a mass of the order of ∼GeV-TeV and fast decoupling from the plasma in the primitive Universe. Within this scenario, Weakly Interacting Massive Particles (WIMPs) are some of the most studied candidates in the community in a variety of multidisciplinary approaches: in collider searches, direct detection experiments, and indirect searches [4,5]. This work focuses on the indirect detection approach, based on the annihilation of WIMP particles into SM particles, which are expected to create secondary fluxes of cosmic-rays (CRs) and gamma-rays. These fluxes can then be detected by ground-based and space telescopes, such as the Major Atmospheric Gamma-Ray Imaging Cherenkov (MAGIC), the High Energy Stereoscopic System (H.E.S.S.), the High-Altitude Water Cherenkov Observatory (HAWC) or Fermi, among others.

Such observations aim at DM-dominated targets, where the DM print should be higher than in others. Because of this, close targets with a low astrophysical gamma-ray emission are preferred. Amongst the possible candidates, dwarf galaxies [6] are one of the most studied, with a DM mass of ${10}^7$–${10}^{10}$ ${\mathrm{M}}_{\odot }$. Other common candidates are, for example, investigating the possible DM print on the Isotropic Gamma-Ray Background [7] and Galaxy Clusters. The latter, with higher DM masses, are typically further and have an important astrophysical gamma-ray background flux that must be modeled. Finally, the target on which we will focus this short review, based on [8], is the Galactic Center (GC). The GC has the advantage that is the closest target with a high DM content. The counterpart is that this target contains a high astrophysical background emission, in the form of a Galactic diffuse emission and sources.

2. The Galactic Center

Our work is based on the possible TeV DM explanation for the observed flux in the inner GC region [9,10], but with the difference of setting constraints on the DM density profile assuming the thermal relic cross-section $\langle \sigma v\rangle =2-3\times {10}^{-26}$ cm³s⁻¹ instead of setting upper limits on the annihilation cross-section. The reason for this is that there are many uncertainties, in both numerical simulations and observations, regarding the characterization of the DM density profile in the inner parsecs of the GC. Therefore, instead of assuming a DM density profile and obtaining upper limits on the annihilation cross-section, we will set the annihilation cross-section to the thermal relic and compare the results on a wide variety of DM density profiles: cuspy, cored, and spiky DM density profiles. Assuming a thermal relic cross-section and a Power-Law Galactic diffuse emission, an enhancement of $\sim 1000$ over the benchmark Navarro–Frenk–White (NFW) density profile is needed to explain the gamma-ray flux detected by H.E.S.S. in the inner parsecs of the GC, with a fitted DM mass $m_{\mathrm{DM}}\simeq 50$ TeV annihilating into the Z channel [9,10]. This region, extending up to 15 pc around Sgr A*, is coincident with the VIR defined below. We will base our work on these conclusions, but with two main novelties: including more regions extending up to ∼450 pc and modeling the Galactic diffuse emission with the gamma-optimized model computed with DRAGON [11], instead of a simple Power-Law.

In order to do that, we will analyze published data from different observations of the GC by H.E.S.S., defining our five regions of interest (Figure 1): Very Inner Region (VIR) [12], Ridge [13], Diffuse Region [12], Halo [14,15], and the Inner Galactic Survey (IGS) [16]. Note in the figure that only one region observes the very inner GC, the VIR, while the rest of the regions mask the inner $\sim 15$ pc. With this preamble, we will try to have an in-depth study of the different possibilities that can explain this enhancement of ∼1000, considering cuspy DM density profiles and DM spikes. Our most important region is, therefore, the VIR since it is the region where the DM enhancement is needed and, also, coincides with the peak of the DM density profile. As for the rest of the regions, given that they do not observe this peak, only constraints on its size and the outer shape of the density profile can be set.

3. Spectral Modeling

We will assume that the gamma-ray flux detected can be explained by the combination of two fluxes, the Galactic diffuse emission and the DM annihilation flux.

The gamma-ray flux produced by annihilating Majorana-type DM particles has the following form [17]:

$$\frac{d{\Phi }_{\mathrm{DM}}}{dE}=\sum _i^{\mathrm{channels}}\frac{{\langle \sigma v\rangle }_i}{2}\frac{d{N}_i}{dE}\frac{\Delta \Omega {\langle J\rangle }_{\Delta \Omega }}{4\pi m_{\mathrm{DM}}^2},$$

(1)

where ${\langle \sigma v\rangle }_i$ is the thermally averaged annihilation cross-section, $\Delta \Omega$ is the solid angle of the region observed, $m_{\mathrm{DM}}$ is the DM mass of the candidate particle, and the index i refers to the Standard Model channel, which is created in the DM annihilation, i.e., the Standard Model particle created after the annihilation. $\frac{dN}{dE}$ is the gamma-ray spectra produced by the subsequent interactions of the different primary products created in the DM annihilation, given by [18,19]. Finally, the J-factor ${\langle J\rangle }_{\Delta \Omega }$ is where the DM spatial distribution is encoded, since it is the integral of the DM density profile squared ${\rho }_{\mathrm{DM}}^2$ along the line of sight $l\left(\widehat{\theta }\right)$ [17]:

$${\langle J\rangle }_{\Delta \Omega }=\frac{1}{\Delta \Omega }{\int }_{\Delta \Omega }\mathrm{d}\Omega {\int }_{l{\left(\widehat{\theta }\right)}_{min}}^{l{\left(\widehat{\theta }\right)}_{max}}{\rho }_{\mathrm{DM}}^2\left[r\left(l\right)\right]dl\left(\widehat{\theta }\right).$$

(2)

On the other hand, the Galactic diffuse emission is created by a sea of CRs confined in the Galaxy by the turbulent magnetic field. These charged CRs interact with the interstellar medium, emitting a bright diffuse radiation, which can be in the form of gamma-rays (and other wavelengths). Because of this, it needs to be modeled since the contribution to the flux can be dominant in most regions of the GC. In order to take into account this diffuse emission, we will use the Max gamma-optimized model computed with the Diffusion Reacceleration and Advection of Galactic cosmic-rays: an Open New code (DRAGON) [11,20,21], a model that is shown to be consistent with the H.E.S.S. observations in the Galactic Ridge region [22]. This diffuse model is tuned to follow observational data such as local charged CRs, gamma-rays, and other multi-wavelength fluxes in the Galactic plane. For more information about the gamma-optimized model, see [11]. In addition to the DRAGON model, we will leave as a free parameter the renormalization of the total flux B as an $\mathcal{O}\left(1\right)$ parameter that we allow to vary in each region:

$$\frac{d{\Phi }_{\mathrm{Bg}}}{dE}=B^2\frac{d{\Phi }_{\mathrm{DRAGON}}}{dE}.$$

(3)

In order to compute our analysis, following [9,10], we will assume a modelization of the observed flux as a combination of the Galactic diffuse emission and the DM annihilation flux, coming from a thermal WIMP particle annihilating into the Z channel with the benchmark thermal relic cross-section $\langle \sigma v\rangle \simeq 2.2\times {10}^{-26} {\mathrm{cm}}^3{\mathrm{s}}^{-1}$, which yields the observed DM relic abundance in the Universe for a TeV candidate [23]. The choice of this channel is because the best fit is obtained for the VIR when compared to other channels [8,9,10]. This approach is very conservative, but since our final aim is to set constraints on the DM density profile in the inner Galaxy ($r<100$ pc), the fit values, obtained with a ${\chi }^2$ analysis, will serve us as upper limits on the J-factors. This follows the idea that, if a DM density profile gives higher J-factors than the fit ones, the theoretical gamma-ray flux would be greater than that observed, and therefore, it is excluded.

4. Spatial Modeling

As introduced in Equation (2), the DM density profile plays a big role in DM indirect searches, as the flux is proportional to the DM density profile squared. For the purpose of studying constraints on the DM density profile, we focus on comparing the fit J-factors obtained in each region with different models coming form N-body simulations (DM-only, GARR, GARR-I, GARR-I300, GARR-II300, ERIS, MOLL, and EAGLE) and dynamical studies of the Milky Way (McMillan17 [24] and Benito20 [2]). To model it, we used the generalized NFW profile (Equation (4) [25]), which contains different key parameters, allowing us to obtain a wide range of profiles, from cuspy profiles to cored profiles (see

Table 1. Parameters of the different DM density profiles considered in this work, defined in Equation (4). We base our models on simulations and dynamical observations of the Milky Way. As for the simulations, we consider DM-only and hydrodynamical simulations (GARR, GARR-I, GARR-I300, GARR-II300, ERIS, MOLL, and EAGLE) (see [26] and the references within). For the phenomenological models, we follow McMillan17 [24] and Benito20 [2]. Also, the local DM density ${\rho }_{\odot }$ of each simulation is shown, evaluated at $R_{\odot }=8.277$ kpc [27].

Profile	$\mathit{\gamma }$	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{\rho }}_{\mathit{s}}\left({\mathit{M}}_{\odot }{\mathbf{kpc}}^{-3}\right)$	${\mathit{r}}_{\mathit{s}}\left(\mathbf{kpc}\right)$	${\mathit{\rho }}_{\odot }\left(\mathbf{GeV}{\mathbf{cm}}^{-3}\right)$
DM-only	1	1	3	$5.38\times {10}^6$	$21.5$	$0.28$
GARR-I	$0.59$	1	$2.70$	$4.97\times {10}^8$	$2.3$	$0.35$
GARR-I300	$1.05$	1	$2.79$	$1.01\times {10}^8$	$4.6$	$0.35$
GARR-II300	$0.02$	$0.42$	$3.39$	$2.40\times {10}^{10}$	$2.5$	$0.35$
ERIS	1	1	3	$2.25\times {10}^7$	$10.9$	$0.36$
MOLL	$8\times {10}^{-9}$	$2.89$	$2.54$	$4.57\times {10}^7$	$4.4$	$0.31$
EAGLE	$1.38$	1	3	$2.18\times {10}^6$	$31.2$	$0.35$
McMillan17	0–$1.5$	1	3	$1.2\times {10}^8$–$5.3\times {10}^5$	$6.8$–$59.9$	$0.33$–$0.43$
Benito20	$0.1$–$1.3$	1	3	$1.8\times {10}^8$–$2.5\times {10}^6$	$7.0$–$40.0$	$0.41$–$0.71$

for the characterization of all the DM density profiles used in this work).

$${\rho }_{\mathrm{halo}}\left(r\right)=\frac{{\rho }_s}{{\left(\frac{r}{r_s}\right)}^{\gamma }{(1+{\left(\frac{r}{r_s}\right)}^{\alpha })}^{\frac{\beta -\gamma }{\alpha }}},$$

(4)

with $r_s$ being the scale radius of the profile and ${\rho }_s$ the normalization factor, given by the value of the local DM density ${\rho }_{\odot }$; $\alpha$ and $\beta$ define the external slope of the profile and $\gamma$ the internal slope. The characterization of external/internal is mainly determined by the scale radius $r_s$, meaning that, when $r\ll r_s$, the profile is proportional to ${\rho }_{\mathrm{halo}}\propto 1/r^{\gamma }$. As can be seen in Table 1, the values of the scale radius are of the order of kpc, greater than the size of the regions we are considering, so the most important parameters for our analysis are the internal slope $\gamma$ and the renormalization ${\rho }_s$. However, the latter parameter can be related to the local DM density ${\rho }_{\odot }$ directly (see Table 1), so we will focus our analysis on the $\gamma$ parameter.

We also consider different kinds of DM spikes coming from the growth of the Super Massive Black Hole Sgr A*. For this case, we follow the adiabatic spike formalism described in [28,29], where the DM density profile is highly enhanced inside $R_{\mathrm{sp}}$, changing the slope of the profile up to ${\gamma }_{\mathrm{sp}}=(9-2\gamma )/(4-\gamma )$ (for the outer part of the density profile, it behaves as Equation (4) without any modification). This kind of spike is created assuming that Sgr A* has grown adiabatically, i.e., slowly compared to the typical timescales of the system, in the center of the Galaxy and without suffering big mergers during the last $\sim 10$ Gyr [30]. For a more realistic scenario, we have considered the dynamical interactions of the DM particles with the stars of the bulge [31,32]. In this case, which we call star spike, we have the same profile as in the adiabatic spike, but with the difference in the internal slope to ${\gamma }_{\mathrm{star}}=1.5$ inside the radius of influence of Sgr A* ($r_{\mathrm{b}}=2$ pc). For this work, we present all the J-factors in the regions defined in Figure 1 computed with all the profiles in Table 1, and also considering the two types of DM spikes.

5. Results

Following the different models considered, here, we present the fit values from the gamma-ray analysis (

Table 2. Fit values of the gamma-ray flux analysis for the 5 regions of interest. Note that the DM mass value $m_{\mathrm{DM}}$ has been left fixed with the fit value obtained in the VIR, ${36}_{-10}^{+7}$ TeV. For the case of the Halo and IGS, we present the 2$\sigma$ UL. As for the rest of the regions, the fit values and 2$\sigma$ errors are presented. See the text for more details.

Parameters	VIR	Ridge	Diffuse	Halo	IGS
$B^2$	${9.2}_{-0.9}^{+0.8}$	${0.3}_{-0.1}^{+0.2}$	${0.8}_{-0.6}^{+0.2}$	$0.13$	$0.02$
${\langle J\rangle }_{\Delta \Omega }$ (${\mathrm{GeV}}^2$${\mathrm{cm}}^{-5}$)	${2.7}_{-0.9}^{+1.0}\times {10}^{28}$	${2.5}_{0.9}^{1.0}\times {10}^{27}$	${1.1}_{-1.1}^{+3.4}\times {10}^{27}$	$2.5\times {10}^{26}$	$1.7\times {10}^{25}$

), obtained with a ${\chi }^2$ analysis, and upper limits on the J-factors for the three approaches presented (Figure 2): generalized NFW (upper row), adiabatic spike (second row) and star spike (last row). Following our approach, the fit value for the DM mass is $m_{\mathrm{DM}}={36}_{-6}^{+4}$ TeV (1$\sigma$ uncertainty), obtained in the VIR since it is where the DM signal should be more prominent as it is where the DM density profile peaks. However, it has been tested that leaving it as a free parameter in the Ridge region, the same results are obtained, but with greater error bars. For the Diffuse Region, Halo, and IGS, the background diffuse emission is dominant since these regions are defined such that the inner part of the GC, the VIR, is masked. This translates into the fact that, in the Diffuse Region, the 2$\sigma$ region for the fit J-factor is compatible with 0; hence, in the figures, it is labeled as an upper limit. For the Halo and IGS, similar upper limits are obtained for the J-factor, so we only show the results for the Halo in the figures.

As we can see, in Figure 2, we can exclude the computed J-factors above the fit ones. This reason is because a greater J-factor would yield a greater gamma-ray flux than that observed, ruling out, then, the underlying DM density profile. Studying the VIR, for the generalized NFW approach, the crossing point is at $\gamma \simeq 1.3$, so cuspier profiles are ruled out. For the adiabatic spike, all slopes are ruled out (except with $\gamma \ll 1$). Finally, for the star spike, we can rule out profiles with slopes $\gamma \gtrsim 0.8$. However, in the rest of the regions, since they do not contain the peak of the DM density profile, no profiles can be excluded. Thanks to this, we can set constraints for the size of the spike at ${\theta }_{sp}\lesssim {\theta }_{Diff}=(0.15–0.45°)$. This constraint is independent of the spike model used and also serves as a consistency check between our assumptions, the observed gamma-ray flux, and the size and slope of the cusp of the inner DM density profile: a DM cuspy profile in the inner pc of the GC dominating the gamma-ray flux and, for the outer part, a domination of the Galactic diffuse emission.