Testing dark matter with Cherenkov light -- prospects of H.E.S.S. and CTA for exploring minimal supersymmetry

TL;DR

This study assesses the capabilities of H.E.S.S. and CTA gamma-ray telescopes in detecting neutralino dark matter within the p9MSSM framework, highlighting their potential to probe regions beyond direct detection limits and improve sensitivity to TeV-scale dark matter.

Contribution

The paper provides an updated analysis of CTA and H.E.S.S. sensitivities to neutralino dark matter, incorporating recent theoretical developments and a comprehensive parameter space scan.

Findings

H.E.S.S. nearly reaches the thermal WIMP cross section.

CTA will improve sensitivity by at least an order of magnitude.

Both telescopes can detect dark matter below the neutrino floor.

Abstract

We provide an updated and improved study of the prospects of the H.E.S.S. and Cherenkov Telescope Array (CTA) experiments in testing neutralino dark matter in the Minimal Supersymmetric Standard Model with nine free parameters (p9MSSM). We include all relevant experimental constraints and theoretical developments, in particular the calculation of the Sommerfeld enhancement for both present-day annihilation and the relic abundance. We perform a state-of-the-art analysis of the CTA sensitivity with a log-likelihood test ratio statistics and apply it to a numerical scan of the p9MSSM parameter space focusing on a TeV scale dark matter. We find that, assuming Einasto profile of dark matter halo in the Milky Way, H.E.S.S. has already been able to nearly reach the so-called thermal WIMP value, while CTA will go below it by providing a further improvement of at least an order of magnitude. Both H.E.S.S. and CTA are sensitive to several cases for which direct detection cross section will be below the so-called neutrino floor, with H.E.S.S. being sensitive to most of the wino region, while CTA also covering a large fraction of the ~1 TeV higgsino region. We show that CTA sensitivity will be further improved in the monochromatic photon search mode for both single-component and underabundant dark matter.

Figures

Tables

Table 1. Table 1: Parameters of the Einasto, NFW and Cored Einasto DM profiles at the GC considered in this work.

Profiles	Einasto (E)	NFW	Cored Einasto (CE)
${\rho }_s$ (GeVcm-3)	0.079	0.307	0.079
$r_s$ (kpc)	20.0	21.0	20.0
${\alpha }_s$	0.17	$-$	0.17
$r_c$ (kpc)	$-$	$-$	3.0

Table 2. Table 2: Ranges of the p9MSSM parameters used in our scans. All masses and trilinear couplings are given in TeV.

Parameter	Range
bino mass
wino mass
gluino mass	$M_3=20$
trilinear couplings
pseudoscalar mass
$\mu$ parameter
3rd gen. left soft squark mass
3rd gen. right up soft squark mass
3rd gen. right down soft squark mass	$m_{{\tilde{b}}_R}=20$
1st/2nd gen. soft squark masses	$m_{{\tilde{Q}}_{1,2}}=m_{{\tilde{d}}_R,{\tilde{s}}_R}=20$
soft slepton masses
soft slepton masses	$m_{{\tilde{e}}_R,{\tilde{\mu }}_R}=m_{{\tilde{L}}_{1,2}}=20$
ratio of Higgs doublet VEVs
Nuisance parameter	Central value, error
Top pole mass $m_t$ (GeV)	(173.34, 0.76)ATLAS:2014wva

Table 3. Table 3: Selected benchmark points characterized by different properties. The main annihilation mechanism at freeze-out corresponds to the final state with the largest branching ratio. The theoretical value of the cross section σ ​ v 0 𝜎 subscript 𝑣 0 \sigma v_{0} is given as well as the relic density Ω χ ​ h 2 subscript Ω 𝜒 superscript ℎ 2 \Omega_{\chi}h^{2} . CTA sensitivity is reported for the simplified and full scheme computation. The sensitivity is expressed as 95% C.L. upper limits. The LLRTS value is derived according to Eq. ( 2 ) for the given σ ​ v 0 𝜎 subscript 𝑣 0 \sigma v_{0} in the two computation schemes.

Benchmark	BM1	BM2	BM3	BM4	BM5	BM6	BM7
points
$m_{\chi }$	1099	1765	1840	531	1516	2288	997
$\left[\mathrm{GeV}\right]$
Branching	$W^{+}{W}^{-}$ 0.23	$b\overline{b}$ 0.35	$W^{+}{W}^{-}$ 0.64	$b\overline{b}$ 0.85	$t\overline{t}$ 0.17	${\tau }^{+}{\tau }^{-}$ 0.26	${\tau }^{+}{\tau }^{-}$ 0.39
fractions	$b\overline{b}$ 0.23	$W^{+}{W}^{-}$ 0.29	$hA$ 0.14	${\tau }^{+}{\tau }^{-}$ 0.14	$b\overline{b}$ 0.16	$\gamma {\tau }^{+}{\tau }^{-}$ 0.22	$t\overline{t}$ 0.37
	$t\overline{t}$ 0.21	ZZ 0.24	ZH 0.14	$t\overline{t}$ 0.01	$hA$ 0.16	$b\overline{b}$ 0.21	$b\overline{b}$ 0.22
	ZZ 0.20	${\tau }^{+}{\tau }^{-}$ 0.05	$\gamma W^{+}{W}^{-}$ 0.08		ZH 0.16	$\gamma {\mu }^{+}{\mu }^{-}$ 0.14	$\gamma {\tau }^{+}{\tau }^{-}$ 0.01
	Zh 0.06	$\gamma W^{+}{W}^{-}$ 0.04			$W^{+}{H}^{-}$ 0.16	$\gamma e^{+}{e}^{-}$ 0.13
	${\tau }^{+}{\tau }^{-}$ 0.04	Zh 0.03			$W^{-}{H}^{+}$ 0.157	$t\overline{t}$ 0.03
${\Omega }_{\chi }{h}^2$	0.10	0.16	0.13	0.13	-	-	-
Main	${\chi }_1^{\pm }$	$\chi \chi \to \text{SM}$	${\chi }_1^{\pm }$, ${\chi }_2$	$\chi \chi \to \text{SM}$	${\chi }_1^{\pm }$, ${\chi }_2$	$\chi \chi \to \text{SM}$	slepton
mechanism	coann.	$A$-funnel	coann.	t-channel	coann.	t-channel	coann.
$\sigma v_0$	$1.98\cdot {10}^{-26}$	$8.18\cdot {10}^{-27}$	$1.08\cdot {10}^{-26}$	$1.12\cdot {10}^{-26}$	$7.63\cdot {10}^{-28}$	$4.54\cdot {10}^{-31}$	$2.14\cdot {10}^{-32}$
$\left[{\mathrm{cm}}^3{s}^{-1}\right]$
$\sigma v_0^{95\%\text{C.L. (simplified)}}$	$7.64\cdot {10}^{-27}$	$7.82\cdot {10}^{-27}$	$6.59\cdot {10}^{-27}$	$7.29\cdot {10}^{-27}$	$8.84\cdot {10}^{-27}$	$7.82\cdot {10}^{-27}$	$6.65\cdot {10}^{-27}$
$\left[{\mathrm{cm}}^3{s}^{-1}\right]$
$\sigma v_0^{95\%\text{C.L. (full calc.)}}$	$7.92\cdot {10}^{-27}$	$8.05\cdot {10}^{-27}$	$6.06\cdot {10}^{-27}$	$7.47\cdot {10}^{-27}$	$9.34\cdot {10}^{-27}$	$6.50\cdot {10}^{-27}$	$4.87\cdot {10}^{-27}$
$\left[{\mathrm{cm}}^3{s}^{-1}\right]$

Equations

$${\mathcal{L}}_{ijk}\left(s_{ijk},b_{ijk}\right)=\text{Pois}\left(s_{ijk}+b_{ijk},m_{ijk}\right)\text{Pois}\left({\alpha }_{jk}{b}_{ijk},n_{ijk}\right),$$

$${\rm LLRTS}=-2\ln\frac{\mathcal{L}\big{(}s_{ijk},\hat{\hat{b}}_{ijk})}{\mathcal{L}\big{(}\hat{s}_{ijk},\hat{b}_{ijk})}\,,$$

$$\frac{\text{d}\Phi^{\rm DM}_{\gamma}}{\text{d}E}\big{(}\Delta\Omega,E\big{)}=\frac{\sigma v_{0}}{8\,\pi\,m_{\rm DM}^{2}}\frac{\text{d}N_{\gamma}(E)}{\text{d}E}\,\times\,J\big{(}\Delta\Omega\big{)}\,,$$

$$J\big{(}\Delta\Omega\big{)}\equiv\int_{\Delta\Omega}\text{d}\Omega\int_{0}^{\infty}\text{d}s\,\rho_{\text{DM}}\big{(}r(s,\theta)\big{)}^{2}\,.$$

$$\begin{gathered} {\rho }_{\mathrm{E}}(r)={\rho }_{\mathrm{s}}\exp \left\{-\frac{2}{{\alpha }_{\mathrm{s}}}\left[{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{{\alpha }_{\mathrm{s}}}-1\right]\right\} \\ \mathrm{and}{\rho }_{\mathrm{NFW}}(r)={\rho }_{\mathrm{s}}{\left[\frac{r}{r_{\mathrm{s}}}{\left(1+\frac{r}{r_{\mathrm{s}}}\right)}^2\right]}^{-1}, \end{gathered}$$

$$s_{ijk}=T_{\text{obs}}{\int }_{\Delta E_i}dE\frac{d{\Phi }_{\gamma }^{\text{DM}}}{dE}\left(\Delta {\Omega }_{jk},E\right)A_{\mathrm{eff}}^{\gamma }\left(E\right)\mathcal{G}\left(m_{\mathrm{DM}}-E\right)$$

$$n_{ijk}=T_{\mathrm{obs}}{\int }_{\Delta E_i}{\int }_{\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{jk}}}dEd\Omega \frac{d{\Gamma }^{\mathrm{CR}}}{dEd\Omega }\left(\Omega ,E\right),$$

$${\Omega }_{\chi }{h}^2=0.120\pm \mathrm{0.001.}$$

$$\text{BR}\left(B\to X_s\gamma \right)$$

$$\text{BR}\left(B_s^0\to {\mu }^{+}{\mu }^{-}\right)$$

Full text

1 Introduction

Dark matter (DM) is the dominant component of matter in the universe but its nature remains unknown. Dark matter in the form of weakly interacting massive particles (WIMPs) has attracted a great deal of interest in the last decades, and a worldwide experimental effort is underway to unveil its fundamental properties (for recent reviews see, e.g.,Arcadi:2017kky ; Roszkowski:2017nbc ). WIMP candidates appear in many extensions of the Standard Model (SM), among which a notable example is supersymmetry (SUSY). A multi-faceted approach has been developed to search for WIMP DM that exploits the complementarity of direct detection strategies, in which one attempts to detect WIMPs scattering off the target nuclei, indirect detection, which seeks detecting products of WIMP self annihilations, and the production at high-energy colliders.

The so-far null experimental searches carried out at colliders and in underground laboratories for the direct detection of DM have pushed the WIMP mass scale into the TeV range. In the SUSY framework, this is in agreement with expectations for the scale of soft SUSY-breaking consistent with the discovery at the LHC of a 125$\,\mathrm{GeV}$ Higgs boson. To be able to study TeV WIMPs at colliders requires center-of-mass energy beyond that of the Large Hadron Collider (LHC), and direct detection faces the lowered number density of DM particles due to the larger DM mass.

It is in this mass regime where indirect detection with gamma rays may also play a major role. While a continuous part of the gamma-ray flux is expected to drop for energies close to the DM mass, pronounced line-like features that appear there provide a distinctive signature of TeV DM over astrophysical backgrounds. The quest for such spectral features further motivates the searches carried out with instruments with large effective area at TeV energies, such as ground-based arrays of Imaging Atmospheric Cherenkov Telescopes (IACTs). The currently operating IACTs H.E.S.S.Abdallah:2018qtu , MAGICAleksic:2013xea and VERITASArchambault:2017wyh , as well as the Fermi-LATAckermann:2015lka experiment on satellite and the ground-based water tank array HAWCAbeysekara:2017jxs , have done deep observation campaigns in the Galactic Center (GC) of the Milky Way and nearby dwarf galaxy satellites of the Milky Way. The next-generation Cherenkov Telescope Array (CTA) is expected to start data taking within a decade.

The GC region is arguably the most promising astrophysical environment to detect DM annihilation signals in very high energy (VHE, $E\gtrsim 100\,\mathrm{GeV}$) gamma rays due to its relative proximity and the expected large accumulation of DM. However, the GC region is known to be also populated with numerous standard astrophysical emitters in VHE gamma rays. The detection of sharp spectral features expected from TeV DM annihilations would then be key to provide convincing signatures against the smoother energy spectra exhibited by astrophysical backgrounds.

The purpose of this work is to improve and update on previous papersFowlie:2013oua ; Cahill-Rowley:2014twa ; Roszkowski:2014iqa ; Catalan:2015cna ; Aad:2015baa ; Beneke:2016jpw ; Arbey:2017eos ; Athron:2017yua ; Abel:2018ekz that have explored the observational status and prospects of detecting neutralino DM within the Minimal Supersymmetric Standard Model (MSSM). In our analysis we focus on the upcoming CTAAcharya:2017ttl and therefore on the heavy neutralino as DM in the nine-parameter version of the MSSM (p9MSSM) that will be defined below.

Our analysis improves previous works by: (i) deriving the projected CTA sensitivity via a state-of-the-art binned likelihood analysis to be used by the CTA Collaboration, (ii) using up-to-date experimental constraints and numerical tools that include, e.g., 13$\,\mathrm{TeV}$ LHC data and (iii) taking into account the Sommerfeld enhancement (SE) for all points in the scan, whereas previous works included it only as an estimate or only in some selected sectors of parameter space.

The paper is organized as follows. In Sec. 2 we provide an overview of the recent input in VHE gamma-ray results from the H.E.S.S. experiment in the context of searches for heavy DM. An update of the CTA sensitivity to DM searches in the GC region using the latest Monte Carlo simulations of the CTA instrument response functions is provided. In Sec. 3 we briefly describe the p9MSSM, scanning methodology and experimental constraints applied in the analysis. In Sec. 4 we compare the results of our scans with the reach of current and planned indirect and direct detection experiments. We stress the importance of CTA to provide coverage of one of the most interesting cases, the $\sim 1\,\mathrm{TeV}$ higgsino region, as emphasized in Kowalska:2015kaa (see alsoKrall:2017xij ; Kowalska:2018toh for a recent work and review) that otherwise would remain unexplored. Finally, we present our conclusions in Sec. 5.

2 Indirect detection with VHE gamma rays

2.1 Observation of the Galactic Center region with H.E.S.S.

The GC region is the brightest source of DM annihilation in gamma rays (for a review see e.g.vanEldik:2015qla ). However, it harbors numerous astrophysical sources that shine in very high energy ($E\gtrsim$ 100 GeV) gamma rays. Among them are H.E.S.S. J1745-290Aharonian:2009zk , a strong emission spatially coincident with the supermassive black hole Sagittarius A*, the supernova/pulsar wind nebula G09+01Aharonian:2005br , the supernova remnant H.E.S.S. J1745-303Aharonian:2008gw , and a diffuse emission extending along the Galactic planeAharonian:2006au ; Abramowski:2016mir . The rich observational dataset obtained from deep observations of the GC region by the H.E.S.S. phase-I instrument has been used to look for continuum and line signals from DM annihilations. Standard analyses of H.E.S.S.-I observations of the GC region provided about 250 hours of live time in the inner $1^{\circ}$ of the GC.

H.E.S.S. searches have been performed with 10 years of data of the 4-telescope array towards the GC for the continuumAbdallah:2016ygi and mono-energetic gamma lineAbdallah:2018qtu channels, respectively, using a 2-dimensional likelihood ratio test statistics to look for any possible excess over the measured background. In order to avoid modeling the complex standard astrophysical background in the GC region, a region of $\pm 0.3^{\circ}$ in Galactic latitude along the Galactic plane has been excluded from the dataset together with a disk of $0.4^{\circ}$ radius centered at the position of J1745-303. No excess in the signal region with respect to background was found and some of the strongest constraints on TeV DM were derived in various annihilation channelsAbdallah:2016ygi ; Abdallah:2018qtu .

2.2 Dark matter prospects with CTA in the Inner Galactic halo

The central region of the Milky Way is also a prime target for DM search with the planned CTAAcharya:2017ttl . CTA is envisaged as a two-site observatory to be built at the Paranal site (Chile) in the Southern hemisphere and at La Palma (Spain) in the Northern hemisphere. As the GC region can be observed under favorable and efficient conditions from the Southern hemisphere, the Chilean site of CTA is best suited to explore the GC region. The CTA observation strategy plans a survey of the GC region as a key-science observation program for DM searchesMoulin:2019oyc . A deep multi-year observation program is planned in the form of an extended and homogeneous survey of the inner several degrees of the GC. The CTA flux sensitivity is expected to improve by up to about one order of magnitude compared to H.E.S.S. and the energy resolution reaches 15% at about 100 GeV down to better than 5% in the TeV energy range. The performance of CTA used in this work is based on instrumental Monte Carlo simulations performed for the Southern array which comprises 4 Large-Size Telescopes (LSTs), 25 Mid-Size Telescopes (MSTs), and 70 Small-Size Telescopes (SSTs). See Ref.HASSAN201776 for further details. Following the methodology presented in Ref.Rinchiuso:2018ajn , in this work the sensitivity to DM annihilation signals for CTA observations of the GC region is computed using the latest publicly-available instrument response functions (IRFs) of the CTA Southern site at average zenith angle $20^{\circ}$CTAIRFs .

The main background for IACTs measurements consists of hadronic (proton and nuclei) cosmic rays (CRs) as well as electron and positron CRs, with a dominant contribution from the protons. In order to efficiently separate gamma rays from an overwhelming CR background, efficient discrimination techniques based of the shower image topology have been developedAharonian:2008zza . However, due to finite CR rejection of IACTs, a residual background consisting of misreconstructed CRs identified as gamma rays is unavoidable. The expected residual background determination for CTA has been performed through extensive Monte Carlo simulationsHASSAN201776 . Here, we use the so-called prod3b version of the instrument response functions that include the residual background determinations.

The region of interest for the DM signal extraction with CTA extends up to $\pm 5^{\circ}$ from the GC both in Galactic longitude and latitude. The overall region is split into squared pixels of side $0.5^{\circ}$. A homogeneous exposure of 500 hours is assumed over the entire field of view. The energy-differential residual background rate and acceptance are extracted from Ref.CTAIRFs , and the energy threshold is taken at 30$\,\mathrm{GeV}$. All observations are assumed to be taken at 20∘ zenith angle. The IRFs depends on the chosen analysis cuts. All simulations are based on the CTA-South site performance according to an event selection optimized for 50 hours of observation.

The above-mentioned IRFs are provided for on-axis measurement, i.e. for emission located near the center of the field of view (FoV). In case of emission distant from the center, the IRFs have been computed as a function of the off-axis angle and the CTA flux sensitivity has been computed accordingly. The radius of the FoV region in which the flux sensitivity is within a factor 2 of the one at the center is more than 3 degrees above several hundred$\,\mathrm{GeV}$CTAIRFs . A possible CTA GC survey can make use of a regular grid of pointing positions. Provided that the distance between two nearby pointings is close enough, an overall spatially homogeneous sensitivity can be obtained. At a few degree distance, the sensitivity reached from a single pointing position degrades significantly but is expected to be compensated by nearby pointings. An optimized and quantitative pointing position strategy for the GC survey with CTA to achieve the best possible sensitivity in the inner several degrees of the GC is much beyond the scope of this work. In what follows we will assume that an homogeneous flux sensitivity in the overall region of interest with a 500 hour flat exposure can be achieved provided that the overall adequate amount of observation time is granted to the GC survey to fulfill this goal.

2.3 Statistical method for sensitivity computation

A dedicated 3-dimensional likelihood ratio test statistics technique has been developed to exploit the spectral and spatial features of the expected DM signal with respect to the background. The spatial pixels are defined as squared pixels of $0.5^{\circ}$ between $\pm 5^{\circ}$ in both Galactic longitude and latitudes. 20 energy bins are taken logarithmically-spaced between energies from 10$\,\mathrm{GeV}$ to 100$\,\mathrm{TeV}$, following Ref.CTAIRFs .

The likelihood function for DM searches is defined as a product of the Poisson probabilities of event counting in the signal and background regions in the $i$-th energy bin, $j$-th Galactic longitude bin, and $k$-th Galactic latitude bin. It reads

[TABLE]

where the likelihood function follows a Poisson distribution given by $\textrm{Pois}(\lambda,n)=\lambda^{n}e^{-\lambda}/n!$ while $\alpha_{jk}$ corresponds to the ratio of the solid angle size of the background over the signal regions. The measured count numbers in the signal and background regions are $m_{ijk}$ and $n_{ijk}$, respectively. Following H.E.S.S.’ strategy for DM searches in the GC region, the OFF region (see below) measurements are taken in the same observational and instrumental conditions as for the signal measurement, which does not require any further acceptance correction for the background determination. In this case, $\alpha_{jk}$ is taken to be 1. In the context of a counting experiment using ON-OFF measurementsAbdallah:2018qtu , the signal is searched in an ON region where the measured number of events is $m_{ijk}$. The expected number of background events in the signal region, $b_{ijk}$, is determined from the measurement of number of events in control (OFF) regions, $n_{ijk}$, with no or little expected searched signal. $s_{ijk}$ is the expected signal in the signal region. In order to compute an expected sensitivity, no excess between the ON and the OFF regions is assumed, i.e. $m_{ijk}\equiv n_{ijk}$.where the likelihood function follows a Poisson distribution given by $\textrm{Pois}(\lambda,n)=\lambda^{n}e^{-\lambda}/n!$ while $\alpha_{jk}$ corresponds to the ratio of the solid angle size of the background over the signal regions. The measured count numbers in the signal and background regions are $m_{ijk}$ and $n_{ijk}$, respectively. Following H.E.S.S.’ strategy for DM searches in the GC region, the OFF region (see below) measurements are taken in the same observational and instrumental conditions as for the signal measurement, which does not require any further acceptance correction for the background determination. In this case, $\alpha_{jk}$ is taken to be 1. In the context of a counting experiment using ON-OFF measurementsAbdallah:2018qtu , the signal is searched in an ON region where the measured number of events is $m_{ijk}$. The expected number of background events in the signal region, $b_{ijk}$, is determined from the measurement of number of events in control (OFF) regions, $n_{ijk}$, with no or little expected searched signal. $s_{ijk}$ is the expected signal in the signal region. In order to compute an expected sensitivity, no excess between the ON and the OFF regions is assumed, i.e. $m_{ijk}\equiv n_{ijk}$.

The total likelihood function $\mathcal{L}$ is the product of the individual likelihood functions over each $ijk$ bin defined as ${\displaystyle\mathcal{L}=\prod_{ijk}\mathcal{L}_{ijk}}$. The log-likelihood ratio test statistics (LLRTS) is defined as:

[TABLE]

where the single and double carets indicate unconditional and conditional maximization, respectivelyCowan:2010js . For each mass, a LLRTS is computed and a LLRTS value equal to 2.71 for one degree of freedom corresponds to an one-sided upper limit at 95% C.L. on $\sigma v_{0}$. The expected sensitivity is computed for a 100% branching ratio in each of the channels $W^{+}W^{-}$, $ZZ$, $hh$, $Zh$, $c\bar{c}$, $b\bar{b}$, $t\bar{t}$, $e^{+}e^{-}$, $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$ and $\gamma\gamma$.

2.4 Expected signal and background events

The expected photon flux from pair-annihilation of DM particles of mass $m_{\rm DM}$ in a region of solid angle $\Delta\Omega$ in the sky can be expressed as

[TABLE]

where $\sigma v_{0}$ is the total annihilation cross section to all primary channels providing photons in the final sates, and ${\rm d}N_{\gamma}(E)/{\rm d}E$ is the photon spectrum per annihilation. $J(\Delta\Omega)$ is the so-called J-factor defined as the integral of the square of the DM density $\rho$ along the line of sight $s$ and over $\Delta\Omega$ by

[TABLE]

$s$ is the distance along the line of sight from the observer and is related to the radial distance $r$ in the coordinates centered at GC by r=\big{(}s^{2}+r_{\odot}^{2}-2\,r_{\odot}\,s\,\cos\theta\big{)}^{1/2}, where $\theta$ is the angle between the direction of observation and the GC, and $r_{\odot}$ = 8.5 kpc is the distance from the Sun to the GC. We consider a cuspy DM distribution at the GC for which suitable parametrizations are the Einasto1965TrAlm…5…87E and Navarro-Frenk-White (NFW)Navarro:1996gj profiles defined as

[TABLE]

where normalization $\rho_{s}$, scale radius $r_{s}$, and power index $\alpha_{s}$ are given in Table 1, following Ref.Pieri:2009je . The local DM density is taken to be $\rho_{\odot}=0.39\ \rm GeV\ cm^{-3}$Catena:2009mf . Since the DM density in the GC is rather uncertain we also consider a Cored Einasto profile with a core radius $r_{\rm c}$ such that $\rho_{\rm CE}(r<r_{\rm c})=\rho_{\rm E}(r_{\rm c})$ and $\rho_{\rm CE}(r\geq r_{\rm c})=\rho_{\rm E}(r)$. We note that the presence of possible DM substructures is known to play a subdominant role in DM searches towards the GC and is not therefore considered here (for a discussion see, e.g., Ref.Gaskins:2016cha and references therein).

The expected DM signal count number in the $ijk$-th bin writes as

[TABLE]

where $d\Phi^{\rm DM}_{\gamma}/dE$ is defined in Eq. (3), $A_{\rm eff}^{\gamma}(E)$ is the gamma-ray energy-dependent effective area, $T_{\rm obs}$ is the observation time, and $\mathcal{G}(m_{\rm DM}-E)$ is a Gaussian function centered at the DM mass $m_{\rm DM}$ of width taken as the CTA energy resolution in order to reproduce the effect of the energy resolution on the theoretical signal spectrum. The DM spectrum $dN_{\gamma}(E)/dE$ is taken from Ref.Cirelli:2010xx for continuum channels. The monoenergetic gamma-ray line is a Dirac delta function centered at $m_{\rm DM}$. The signal count rate in the $ijk$-th bin is integrated over the spatial pixel of solid angle size $\Delta\Omega_{jk}$ and energy bin width $\Delta E_{i}$.

The CR background count number is given by

[TABLE]

where $d\Gamma^{\rm CR}/dEd\Omega$ is the energy-differential residual background rate per steradian. The background modeling follows the Monte Carlo procedure outlined in the papers, Refs.APJ2012 ; HASSAN201776 .

A detailed modeling of the spectral and spatial extrapolation of the Galactic Diffuse Emission measured by Fermi-LAT in the TeV energy range is beyond the scope of this paper and we neglect it in our computation. A band of $\pm 0.3^{\circ}$ in Galactic latitudes is excluded from the ROIs as being populated by numerous standard astrophysical sources of VHE gamma rays. A 0.4∘ radius disk is removed at the position of HESS J1745-303, one of the brightest TeV gamma-ray sources in the overall ROI.

We present in Fig. 1a 111All sensitivity limits, including more annihilation channels, for all three halo profiles considered here, can be found in the supplementary material on the arXiv. Limits provided there also extend to 100$\,\mathrm{TeV}$. the projected CTA 95% C.L. sensitivity to DM annihilation as a function of DM mass $m_{\chi}$. For this figure the DM particle is assumed to annihilate into the specific SM final states described in the legend with 100% branching fraction. The exclusion lines are computed for each of the three different choices of the DM Galactic halo profile. The CTA 95% C.L. sensitivity to monochromatic $\gamma$-ray lines for the same three choices of halo profile is featured in Fig. 1b. Note that the current monochromatic $\gamma$-ray H.E.S.S. bound, Ref.Abdallah:2018qtu , is more constraining than the corresponding Fermi-LAT monochromatic bound, Ref.Ackermann:2015lka , in the mass range m\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}300\,\mathrm{GeV}.

3 The MSSM and details of the scan

Low energy-scale supersymmetry (SUSY) is the most thoroughly studied scenario for new physics that provides solutions to the problems of the SM - e.g. hierarchy problem, lack of DM candidate, unification of gauge interactions. Despite null results for any new physics signal at the LHC or direct and indirect detection experiments searching for DM, SUSY remains an attractive candidate for new physics, especially in light of the discovery at the LHC of a Higgs boson with mass not far above the $Z$ boson mass. Indeed, the experimental data have so far only excluded models based on optimistic expectations founded on purely theoretical, or aesthetic, arguments, like naturalness.

3.1 The p9MSSM

The simplest realization of SUSY that is also phenomenologically viable is the R-parity conserving MSSM, where the lightest supersymmetric particle (LSP) is stable and may be identified as a thermally-produced DM candidate. We will assume that the LSP is the lightest neutralino. However, with over 100 free soft-breaking parameters, it is almost impossible (nor actually even necessary) to study the MSSM in complete generality. Therefore, one has to study more constraining models with a specific high-scale mechanism for SUSY breaking (e.g. CMSSM/mSUGRA) or to consider a p(henomenological)MSSMDjouadi:1998di ; Cahill-Rowley:2014twa . The latter is based on the following assumptions: (i) CP conservation, (ii) Minimal Flavor Violation at the electroweak scale, (iii) degenerate first two generations of sfermion soft-mass parameters and (iv) negligible Yukawa couplings and trilinear couplings ($A$-terms) for the first two generations. In our numerical scan we consider the p9MSSM where in addition to the above, we set the gluino mass, the third-generation down-type right soft squark mass, and the first two generations of soft slepton masses at $20\,\mathrm{TeV}$, which decouples them (see Table 2). The p9MSSM provides a sufficiently generic parametrization and coverage of the DM properties of the MSSM with CP and R-parity conservation. It captures a rich electroweak scale phenomenology with multiple possibilities regarding its UV-completion, while being sufficient for our purpose of exploring heavy neutralinos as DM. Indeed, adding more MSSM parameters to the scan would not alter our results in any significant way.

3.2 p9MSSM scanning setup and constraints

We apply the projected sensitivity reach of CTA as calculated in Sec. 2 to the case of the MSSM parametrized by 9 free input parameters. The parameters that we scan over and their ranges are shown in Table 2. We employ the $\tt Multinest\_v.2.7$Feroz:2007kg ; Feroz:2008xx package for the scan, using flat priors. In order to ensure the best coverage of the parameter space of the model, several independent scans with $20,000$ live points each have been performed and the resulting points have been combined when presenting the results. The supersymmetric spectrum is calculated with SPheno v4.0.3Porod:2003um ; Porod:2011nf . We allow the bino mass $M_{1}$ and the $\mu$ parameter to assume negative values in order to accommodate blind spots in DM direct detectionEllis:2000ds ; Cheung:2012qy , which stem from the vanishing $h\chi\chi$ coupling for certain combinations of parameters (see alsoHan:2016qtc for a recent discussion). The remaining gaugino mass parameter $M_{2}$ is kept positive, starting from the a minimal value of 100$\,\mathrm{GeV}$, allowed by the LEP bounds on charginos. Most of the third generation sfermion masses are allowed to assume a broad range of values in between being almost mass degenerate with the lightest neutralino up to tens of TeV. The former regime allows for efficient co-annihilations to occur in the early Universe when the DM relic density is determined, while the latter, in case of squarks, can more easily lead to a correct value of the Higgs boson mass, $m_{h}$, thanks to an increase of the characteristic SUSY scale. As discussed above, the remaining sfermion mass parameters and the gluino mass $M_{3}$ are fixed at 20$\,\mathrm{TeV}$. They do not play any real role in a further discussion.

The SUSY mass parameters are defined at the scale of the geometrical average of the physical stop masses, $M_{\rm SUSY}=(m_{\tilde{t}_{1}}m_{\tilde{t}_{2}})^{1/2}$. The ratio of the Higgs doublets’ vevs, $\tan\beta$, and the top quark pole mass, $m_{t}$, which is treated here as a nuisance parameter, are defined at the electroweak symmetry breaking (EWSB) scale. We assume a Gaussian distribution for $m_{t}$, whose central value and experimental error are given inATLAS:2014wva , $m_{t}=(173.34\pm 0.76)\,\mathrm{GeV}$. Our numerical scans are driven by a global likelihood function, which incorporates a standard set of constraints described below.

Dark matter relic density

The constraint with the strongest impact on our numerical result is given by the measurement of the relic abundance of DM, as given by PlanckAghanim:2018eyx ,

[TABLE]

To calculate the relic density we employ $\tt micrOMEGAs\ v.5.0.6$Belanger:2001fz ; Belanger:2004yn supplemented by $\tt DarkSE$Hryczuk:2011tq . We additionally impose a $10\%$ theoretical uncertainty on the calculation to partially take into account the effects of, e.g., loop correctionsBaro:2007em ; Baro:2009na , variations in the renormalization scheme and scaleHarz:2016dql , and modifications to the QCD equations of stateHindmarsh:2005ix ; Laine:2006cp ; Drees:2015exa .

At the typical mass scale tested by CTA and H.E.S.S. ($\sim 1$ to a few TeV) the SE plays an important role and strongly affects both the calculation of the present-day neutralino annihilation cross section, $\sigma v_{0}$, and, to a lesser degree also the determination of the thermal neutralino relic densityHisano:2006nn ; Cirelli:2007xd ; Hryczuk:2010zi . An accurate treatment of the freeze-out process thus requires the incorporation of the SE coming from multiple exchanges of all the gauge bosons and of the SM Higgs and applied to all co-annihilation channels. At present, the only public code that gives the relic density with the SE included for a generic neutralino and all possible co-annihilation partners in the general MSSM is DarkSE - a package written for DarkSUSY v5Bringmann:2018lay .222DarkSE does not take into account some recent theoretical developments relative to the most proper way of implementing the SE computation. In particular, the code includes off-diagonal terms in the annihilation matrixBeneke:2012tg ; Beneke:2014gja ; Beneke:2014hja exclusively in the pure wino limit. In this respect it provides a less accurate determination than that of a new program that is currently being developedBeneke:inprep , which has been already used in several phenomenological studiesBeneke:2016ync ; Beneke:2016jpw . However, DarkSE also presents an additional functionality of having the SE implemented for sfermion co-annihilation, which is a necessary ingredient for the scan performed in this work.

A complete numerical treatment of the SE is very CPU-expensive and thus cannot be handled automatically in a scan. Therefore, we have adopted a two-step approach: 1) in the scan we use $\tt micrOMEGAs$ and include the SE by rescaling the result using a grid of the enhancements in the $M_{2}$-$\mu$ plane following the procedure ofRoszkowski:2014iqa ; 2) the final points are then post-processed with the accurate SE treatment using full $\tt DarkSE$ code. Sommerfeld enhancement is also included in the computation of the present-day $\sigma v_{0}$, as well as for $\sigma v_{\gamma\gamma}$ and $\sigma v_{Z\gamma}$.333It has been checked that the zero-velocity limit of these cross sections gives essentially the same result as when averaged over the Maxwellian velocity distribution of DM in the GC, with only minimal percent level differences in the close proximity of the SE resonance in the wino region. Ideally, one could use an approximate simplified treatments of the SE for the first step, as in, e.g., Slatyer:2009vg ; ElHedri:2016onc , but unfortunately these are known only for simple setups - there is no known method for estimating the SE for the relic density with co-annihilations, and there is also no simple functional dependence on either the input nor physical parameters.

Note that in our analysis we do not take into account possible bound-state formation of strongly interacting co-annihilating particles. This effect was noticed and first discussed for a simple toy model in a recent workHarz:2018csl and potentially can apply to the regions of the MSSM parameter space featuring one or more squarks almost degenerate in mass with the neutralino, particularly if the latter lies around the TeV scale. Implementing bound-state formation in our code would go far beyond the scope of this analysis. While this effect might modify the value of the predicted relic density for some points, these could only be sporadic cases with very strong co-annihilations with squarks.

Another effect that in principle could have some impact on the discussed limits is the modification of the end point of the energy spectrum of photons produced in the present-day DM annihilation due to soft and collinear gauge boson emission. Such processes, though formally of higher order, are enhanced by large Sudakov logarithms especially at energy scales much larger than the electroweak one. This has been noticed in the context of DM annihilation inCiafaloni:2010ti and explicitly seen in the wino annihilation computation at one-loopHryczuk:2011vi while finally approached with resummation techniques inBaumgart:2014vma ; Bauer:2014ula ; Ovanesyan:2014fwa ; Beneke:2018ssm ; Beneke:2019vhz . InBaumgart:2017nsr ; Baumgart:2018yed it was argued that for the neutralino annihilation the modification due to fully resummed exclusive cross section is most relevant in multi-TeV regime. However, more recentlyBeneke:2018ssm ; Beneke:2019vhz showed that for the nearly whole energy regime of interest for CTA, high precision calculations in fact would require resummation of Sudakov logarithms. Nevertheless, due to the fact that it is currently not possible to directly apply the framework ofBaumgart:2017nsr orBeneke:2018ssm to generic p9MSSM neutralinos, and since the expected corrections in the DM mass range of our interest are typically much lower than astrophysical uncertainties, we do not include this effect in our scan.

When performing the numerical scans, we study two commonly discussed cases:

the thermal relic density saturates Eq. (8), in which case we use a Gaussian distribution for the relic density, 2. 2.

the thermal relic density does not exceed the value given in Eq. (8), in which case we use a half-Gaussian distribution - with relic density imposed only as an upper bound.

In the former case, we assume that no deviations from the standard cosmological history of the Universe took place, as well as that the lightest neutralino is the only DM particle. In the latter case, the neutralino cannot be a single particle comprising the DM. We then present results of CTA sensitivity to underabundant neutralinos with local density rescaled by the square of the ratio of the neutralino density to the PlanckAghanim:2018eyx value.

On the other hand, the neutralino relic density can also be significantly affected by deviations from the standard cosmological history of the Universe, e.g., if neutralino freeze-out occurs during an extended reheating periodGiudice:2000ex ; Fornengo:2002db ; Gelmini:2006pq (see alsoRoszkowski:2014lga ; Drees:2018dsj for recent studies) or in presence of additional non-thermal production. In this case, the neutralino can be a single DM particle even though its standard freeze-out relic density does not saturate the Planck value. In order to accommodate such scenarios, we additionally present results for the aforementioned case 2 but without rescaling $\sigma v_{0}$.

Dark matter direct detection

The steady progress observed in recent years in direct detection (DD) searches for DM in underground liquid noble gas detectors has led to strong upper limits on the spin-independent cross section of the neutralino scattering off nucleons.

We include the most recent DD bounds here. For this we employ $\tt SuperIso\ Relic\ v4.0$Arbey:2018msw and the $\tt DDCalc\_v.2.0.0$ packageWorkgroup:2017lvb , assuming the Standard Halo Model (SHM) and the following values for the relevant astrophysical parameters: $\rho_{0}=0.39\leavevmode\nobreak\ \textrm{GeV}/\textrm{cm}^{3}$, $v_{\textrm{rot}}=220\leavevmode\nobreak\ \textrm{km}/\textrm{s}$, $v_{\textrm{esc}}=544\,\textrm{km}/\textrm{s}$. We note that slight modifications to the SHM that might be suggested by e.g. recent data release by the GAIA CollaborationBrown:2018dum , see e.g.Necib:2018iwb ; Evans:2018bqy for further discussion, would have minor impact on our results. The experimental limits that we take into account are the following: PandaX-2Cui:2017nnn , PICO-60Amole:2017dex , and the most recent results from the XENON1T collaborationAprile:2018dbl .

Almost universally in the parameter space of the MSSM the bounds on the spin-dependent cross section of the neutralino scattering off the proton or neutron cannot compete with the corresponding DD bounds on the spin-independent cross section. The current bounds on $\sigma^{\rm SD}_{p}$ for $m_{\chi}$ of our interest come from the searches by the IceCube Collaboration for neutrinos coming from the center of the Milky WayAartsen:2017ulx , the EarthAartsen:2016fep , and the SunAartsen:2016zhm . Since they are less constraining than the aforementioned DD bounds, and also the indirect detection bounds described below, we do not consider them here.

In the future, the neutralino scattering cross section on neutrons and protons can also be constrained by their interactions inside neutron stars and white dwarfsKrall:2017xij ; Baryakhtar:2017dbj . The corresponding limits, however, depend on additional astrophysical assumptions, as well as progress in observations and, therefore, they are not discussed further below.

Collider constraints

The TeV mass-range particle spectrum of the MSSM is very poorly constrained by direct SUSY searches at colliders (see, e.g.,Kowalska:2016ent ; Arbey:2017eos ; Athron:2017yua ; Athron:2018vxy and references therein), including the most recent data from the LHC. In our case, since we focus on the parameter space characterized by colored sparticles lying in the multi-TeV range, only very few points are affected by LHC bounds, with negligible impact on the results shown in Sec. 4. For completeness, we also take into account LEP and Tevatron limits on SUSY particlesTanabashi:2018oca .

Higgs physics

The Higgs mass determination and Higgs-sector LHC measurements in general can show their effect on the MSSM parameter space under probe in DM searches. Indirect constraints on the stop mass and mixing from the Higss mass measurement affect the extension of the regions potentially subject to stop co-annihilation; bounds on the mass and couplings of heavy Higgs bosons can end up influencing somewhat the shape of the funnel regions. In here, the Higgs sector is constrained with $\tt HiggsBounds-5.2.0beta$Bechtle:2008jh ; Bechtle:2015pma and $\tt HiggsSignals-2.2.1beta$Bechtle:2013xfa , while additional constraints from searches for heavy Higgs decays to $\tau^{+}\tau^{-}$ are implemented followingATLAS:2016fpj ; CMS:2016rjp .

Flavor physics

We calculate a few flavor observables with $\tt Superiso\ Relic\ v4.0$Arbey:2018msw . The parameter space of the MSSM is potentially sensitive, in particular, to the bounds from rare decays in $b\to sll$ processes and radiative decays like $b\to s\gamma$, which can constrain scan points characterized by large $\tan\beta$ values, and/or relatively light non-SM Higgs bosons, squarks, and charginos/neutralinos. We use the following experimental determinations:

[TABLE]

where, following, e.g., Ref.Workgroup:2017myk , in Eq. (9) we give the calculated averageMisiak:2017bgg of the determinations in Refs.Aubert:2007my ; Lees:2012wg ; Lees:2012ym ; Saito:2014das ; Belle:2016ufb , and in Eq. (10) we report the most recent LHCb measurement, based on 8$\,\mathrm{TeV}$ collision dataAaij:2017vad . We thus implicitly assume that Eq. (10) has superseded an older statistical combination of CMS and LHCb measurements with 7 and 8$\,\mathrm{TeV}$ dataCMS:2014xfa . Note that very recently the ATLAS Collaboration has presented a measurement of $\textrm{BR}\left(B_{s}^{0}\rightarrow\mu^{+}\mu^{-}\right)$, from a combination of data taken during their 8$\,\mathrm{TeV}$ and 13$\,\mathrm{TeV}$ runs, which agrees with Eq. (10): $\textrm{BR}\left(B_{s}^{0}\rightarrow\mu^{+}\mu^{-}\right)=(2.8_{-0.7}^{+0.8})\times 10^{-9}$Aaboud:2018mst .

We impose the bounds of Eqs. (9) and (10) at the 95% C.L., a posteriori on the points belonging to the $2\,\sigma$ region of the profile likelihood. This reduces the number of viable points in the scan by approximately 2%. Other potentially relevant flavor observables like $\textrm{BR}(B^{\pm}\to\tau\nu_{\tau})$ or the $B_{s}$ mass mixing measurement $\Delta M_{B_{s}}$ are not constraining at the mass scale relevant for this paper.

Note that we do not include constraints from observables that are currently showing a 2–3 $\sigma$ discrepancy with the SM, like the differential branching ratios and angular observables in $B^{0}\to K^{\ast 0}\mu^{+}\mu^{-}$Aaij:2014pli ; Aaij:2015oid , or the branching ratio measurements that have recently provided tantalizing hints of lepton flavor nonuniversalityLees:2013uzd ; Aaij:2015yra ; Hirose:2016wfn ; Aaij:2017vbb ; Aaij:2019wad . It is known that these anomalies cannot be explained consistently in the MSSM (see, e.g., Ref.Altmannshofer:2014rta ) and that, if confirmed to higher statistical significance with further release of data, will require new physics beyond the particle content of the MSSM. For analogous reasons, we do not apply to the parameter space the constraint from the measurementBennett:2006fi of the muon anomalous magnetic moment, which shows a $3.5\,\sigma$ discrepancy with the SM expectation, $\delta\left(g-2\right)_{\mu}=(27.4\pm 7.6)\times 10^{-10}$Davier:2016iru . It is well known that this value cannot be accommodated in the regions of the MSSM parameter space that feature a TeV-scale LSP, see, e.g., Ref.Kowalska:2015zja . In this case too, if the anomaly were to be confirmed by upcoming Fermilab dataGrange:2015fou , it will require a BSM explanation lying outside of the MSSM parameter space relevant for the current analysis. One should keep in mind, however, that is possible to extend the MSSM minimally by a U(1) gauge group, so that the all of the above-mentioned flavor anomalies become consistent with TeV-scale neutralinos with the exact same DM properties as in the vanilla MSSMDarme:2018hqg .

Dark matter indirect detection

The indirect detection constraints on neutralino DM, that are the main subject of this study, are not included in the likelihood function when performing initial numerical scans of the parameter space of the MSSM. Instead, we carefully study them by postprocessing the results obtained in these scans. This leads to a better understanding of their impact on the allowed parameter space.

The most constraining data for the TeV-scale mass range are currently provided by H.E.S.S. A more detailed description of DM ID limits from H.E.S.S. and future projections has been described in details in section 2.

When presenting the results below, we also take into account Fermi-LAT limits on DM-induced $\gamma$-rays that correspond to 6 years of data and observation of $28$ dSphsFermi-LAT:2016uux . These data are in principle most constraining in the MSSM for neutralinos of mixed gauge composition with a mass of a few hundred GeV, which are, however, already strongly bounded by the null DD results. They might also provide a complementary probe on the low-energy tail of spectra from the annihilation of winos including SE. We illustrate this below for a fixed annihilation final state into a $b\bar{b}$ pair. We have also verified numerically, following $\tt Superiso\ Relic\ v4.0$Arbey:2018msw , that taking into account a complete list of annihilation final states leads to similar results.

We note that $\sigma v$ can also be constrained by requiring that the CMB spectrum is not affected too much by the pre- and post-recombination energy injection from DM annihilationsPadmanabhan:2005es ; Galli:2009zc ; Slatyer:2009yq . However, for the heavy DM of our interest, this effect typically leads to less stringent bounds than null searches for DM annihilation signal in the GC by H.E.S.S. and in dSphs by Fermi-LAT (seeBeneke:2016jpw for recent discussion).

A recent determinationCuoco:2017iax of the bounds on the neutralino annihilation cross section from AMS-02 antiproton cosmic-ray (CR) dataAguilar:2016kjl has proven to be competitive with H.E.S.S. diffuse $\gamma$-ray searches in the $\sim{\rm TeV}$ mass range. We have verified that this is in general true also in the context of our scans using $\tt SuperIso\ Relic\ v4.0$Arbey:2018msw that employs a semi-analytic approach to solving the propagation equations followingBoudaud:2014qra . However, the limits that one derives from the AMS-02 data depend on the assumed CR propagation model and suffers from large astrophysical uncertainties (see, e.g., Refs.Evoli:2011id ; Cirelli:2013hv ). For this reason, we do not discuss them in details in the following section, which focuses on DM-induced $\gamma$-ray signal.

4 Results

4.1 CTA sensitivity to the p9MSSM parameter space

For each point in the scan, we compute the H.E.S.S. limit for the present-day annihilation cross section, $\sigma v_{0}$, and the corresponding projected sensitivity of CTA. We use the 95% C.L. bounds and projections for annihilation to pure channels (see Fig. 1a). In case of annihilation final states for which H.E.S.S. limits have not been reported by the collaboration, we employ the most relevant existing bounds. In particular, for $hh$ final state we use $ZZ$ limit, for final states with $c$ and $s$ quarks – $b\bar{b}$ limit, for the lightest quarks – $\tau^{+}\tau^{-}$ limit and for $e^{+}e^{-}$ we employ $\mu^{+}\mu^{-}$ limit. Instead, for CTA we derive bounds for a more complete set of annihilation final states, as discussed in section 2.4.

In order to verify whether a particular point in the p9MSSM parameter space is within current bounds and future sensitivities, we combine limits obtained for pure annihilation final states by taking their average weighted by the branching-ratios corresponding to those channels. In section 4.4 we show that this procedure is sufficient for our purpose, by comparing our results for several benchmark scenarios to a more detailed treatment in which photon spectra are carefully combined prior to obtaining the CTA limit.

For channels with non-SM particles in the final state, e.g., the neutral MSSM Higgs particles, $A_{0}$ and $H$, we employ the bounds computed for the SM Higgs; for the charged MSSM Higgs particle, $H^{\pm}$, we use the bounds derived for $W^{\pm}$. While the non-SM annihilation final states typically do not play a dominant role in our analysis, they might become important for selected points in the parameter space. For these points, we have verified our results against a more detailed procedure in which decays of the non-SM particles were taken into account employing $\tt HDECAY$Djouadi:1997yw ; Djouadi:2018xqq before generating the combined photon spectrum using Ref.Cirelli:2010xx .

In the plots below we only show points that belong to the 95% C.L. region of the global profile-likelihood, i.e. we select $\Delta\chi^{2}\leq 5.99$ from the best-fit point, where $\Delta\chi^{2}=-2\ln\left(\mathcal{L}/\mathcal{L}_{\max}\right)$.

4.2 Discussion of results

We present in Fig. 2 the scan points in the plane ($m_{\chi}$, $\sigma v_{0}$) of the present-day annihilation cross section of the neutralino versus its mass. The color code used in Fig. 2 refers to the gauge composition of the neutralino LSP, which, by construction, in the MSSM is never a 100% pure eigenstate.

How “pure” a certain mass eigentate is depends on the elements of the unitary matrix, $\boldsymbol{Z}$, diagonalizing the neutralino mass matrix after EWSB. In green we show the points with the LSP containing at least 90% of the pure bino gauge eigenstate (i.e., that is, in the basis of gauge eigenstes {bino, wino, down-type higgsino, up-type higgsino}, we require $|Z_{11}|^{2}>0.9$ for these points). In blue, the points for which it is for at least 90% a wino ($|Z_{12}|^{2}>0.9$). In cyan we show a mixture of these two gaugino states, with the additional constraint that the higgsino composition remain below 10%, $|Z_{13}|^{2}+|Z_{14}|^{2}<0.1$. In red we show neutralinos that are dominated for at least 90% by their higgsino fraction ($|Z_{13}|^{2}+|Z_{14}|^{2}>0.9$). We finally point out that only very few points characterized by a mixture of a gaugino and a large higgsino component appear, marked in gold, in the plot, as they are in strong tension with the latest bounds from direct detection searches.

We also note that due to our general focus on TeV-scale neutralino DM, which is of most relevance for H.E.S.S and CTA, our scanning procedure is not tuned to thoroughly explore the parameter space of the p9MSSM corresponding to light neutralinos with masses around the EWSB scale. For this reason, we do not show in our plots points corresponding to the region where $m_{\chi}\approx m_{h}/2$ where the correct neutralino DM relic density can be obtained thanks to efficient resonance annihilations via the Higgs boson exchange. We note, however, that the expected annihilation cross section for such light neutralino DM lies well below the reach of CTA. The same is also true for another instance of supersymmetry at the electroweak scale that has recently been discussed in the context of a collection of mild excesses present in the LHC dataAthron:2018vxy .

We show in Fig. 2 with a solid black line the current 95% C.L. upper bound on $\sigma v_{0}$ from 254 hours of observation of the GC at H.E.S.S., under the Einasto profile assumption, applied to the points of the p9MSSM. Importantly, when deriving these results, as well as CTA sensitivity lines discussed below, we take into account all the points obtained in the scan of the parameter space including the ones that violate some of other constraints and, therefore, are not shown in the plot. In particular, the presence of such otherwise excluded points allows us to determine the position of the H.E.S.S. limit in the region with a light neutralino and large $\sigma v_{0}$ which is virtually excluded by current bounds.

The latest observations exclude points whose neutralino is strongly dominated by the wino component (in blue, and some in cyan), for which the annihilation cross section in the present day has a large SECohen:2013ama ; Hryczuk:2014hpa . The plot updates Figure 5(a) of Ref.Roszkowski:2014iqa and is in agreement with e.g. Ref.Catalan:2015cna . Compositions of the neutralino very close to a pure wino state are in very strong tension with H.E.S.S. with continuum observations as well as monochromatic line searches.

The H.E.S.S. bound, on the other hand, does not bite into the $\sim 1\,\mathrm{TeV}$ (nearly pure) higgsino region of the parameter space, corresponding to the red points in Fig. 2, for which the SE is less pronounced. Upcoming increased statistics can tighten the bound but, realistically, batches of new data are at this point not expected to bring qualitative improvements to the current picture. It is CTA, with an effective area that is by about a factor 10 larger than that of H.E.S.S.’ at 1 TeV, that will be probing more deeply into the higgsino region of the parameter space.

We show with a dash-dotted black line our projection of the sensitivity of CTA in the p9MSSM in searches for DM-induced diffuse photon flux, with 500 hours of observation of the GC and under the Einasto profile assumption. In addition, we overlap gray triangles to the points that are within the sensitivity of the CTA $\gamma$-ray line search. As was described in Sec. 2, we factor in a detailed treatment of the statistical uncertainties, and the likelihood function is calculated with an improved design of the ROIs of the Galactic Plane with respect to previous analysesRoszkowski:2014iqa ; Lefranc:2015pza . The higgsino region of the parameter space is likely to be tested in its near entirety by CTA, and the same is true for points with bino-dominated neutralinos with annihilation cross section around the thermal freeze-out value (light dotted line).

We note that the actual CTA limit on the p9MSSM parameter space cannot be perfectly represented by a single line due to number of possible neutralino DM annihilation final states that need to be taken into account. However, we have verified that, for $m_{\chi}\gtrsim 1$ TeV the approximate limit that we present reproduces very well the true CTA sensitivity. For lower masses, the line shown in Fig. 2 corresponds to a conservative approach, i.e., all the points lying above the line are within the CTA reach. We follow a similar strategy to obtain the approximate H.E.S.S. limit shown as a solid line in Fig. 2.

The points that will remain untested, almost all characterized by a nearly pure bino-like composition of the LSP, are those for which the neutralino annihilation cross section is too small to yield the correct relic density, and thus either feature spectra with sparticles that co-annihilate in the early Universe with the LSP (near mass degeneracy between the bino-like neutralino and one or more sfermions), or spectra that include one or more Higgs bosons of mass within a few hundred GeV of $2m_{\chi}$, which provide a means for funnels, or resonant $s$-channel annihilation of the LSP in the early Universe due to the thermal broadening of the energy distribution. As is well known, the specifics of these spectra are very model-dependent. Moreover, their realization in explicit high-scale completions can encounter model building challenges and/or require some fine tuning of the initial parameters.

This is unlike in the case of (nearly pure) higgsinos and winos, which do fall inside the sensitivity of large IACTs, and for which the correct value of the relic density emerges naturally once the mass of the LSP is around either 1$\,\mathrm{TeV}$, or $\sim 2.5-3\,\mathrm{TeV}$, respectively, quite independently of the model details of the rest of the sparticle spectrum. Note, however, that for higgsino points with masses larger than about 1.3$\,\mathrm{TeV}$, shown outside of the CTA sensitivity in Fig. 2, one also relies on additional mechanisms like coannihilations with squarks in the early Universe to preserve the correct relic density. These points tend to feature lower present-day annihilation cross section than lighter higgsinos, and they are consequently more difficult to probe.

The projected sensitivity of CTA shown in Fig. 2 is obtained in the two limiting cases of Einasto and Cored Einasto DM halo profiles. A sensitivity line corresponding to the NFW profile can be easily obtained by multiplying the projected line for Einasto case by the factor of about 2.5 obtained from Fig. 1a.

In Fig. 2, we also show with the dot-dashed line the projected CTA sensitivity reach obtained for the Cored Einasto profile defined in Sec. 2.4. As can be seen, in this case CTA can still play an important role by probing the entire wino-like neutralino DM scenario which would otherwise remained not fully tested by the H.E.S.S. observations.

In Fig. 3 we show the p9MSSM points in the ($m_{\chi}$, $\sigma^{\rm SI}_{p}$) plane. The most recent XENON1T 90% C.L. upper limitAprile:2018dbl is shown by a solid violet line. The XENON1T results are included in the global likelihood function, and that explains the absence above the line of any point belonging to the $\sim 2\sigma$ region of the profile likelihood. The onset of the irreducible neutrino background is denoted by a solid black line. The color code is the same as in Fig. 2 and we additionally overlap violet triangles to points excluded by the H.E.S.S. bound on $\sigma v_{0}$. Black triangles are overlapped to points within our projection of the sensitivity of CTA in the Einasto profile.

The necessity of using both direct and indirect detection strategies to cover the most substantial portions of the parameter space of the MSSM with high-mass DM has been pointed out in the literature since early after the discovery of the Higgs boson at the LHCRoszkowski:2014iqa . We show the power of complementarity of direct and indirect detection in Fig. 4, where we project the points of the p9MSSM to the ($\sigma^{\rm SI}_{p}$, $\sigma v_{0}$) plane. The color code is the same as in Fig. 2 and Fig. 3.

The future reach of direct underground searches with noble liquids is bound to bite into the parameter space from right to left, until it reaches the irreducible neutrino background, shown here as a shaded region (recall that the value of $\sigma^{\rm SI}_{p}$ characteristic of the neutrino “floor” for direct DM searches depends on the DM mass, hence the boundary of the shaded area is jagged in Fig. 4). To guide the eye, we add a vertical dashed gray line denoting the neutrino background limit $\sigma^{\rm SI}_{p}\approx 4\cdot 10^{-12}\ \text{pb}$ taken at $m_{\chi}\approx 2\,\mathrm{TeV}$. Conversely, the sensitivity of IACTs gradually improves from the top down, providing a complementary means of testing the parameter space. The H.E.S.S. bound is denoted in the figure by a dashed black horizontal line while the projected sensitivity of CTA is denoted by a dashed double-dotted horizontal line.

4.3 Underabundant neutralinos

As discussed in Sec. 3.2, the neutralino can be a good DM candidate even when its thermally produced relic abundance is different from the total DM relic density in the Universe. It can then either be one of several DM components, or might even remain the only DM particle but in non-standard cosmological scenarios. In this subsection, we present the results of two scans corresponding to the cases in which the relic density constraint is imposed as an upper bound only, by means of a half Gaussian distribution. The corresponding results can be seen in Figs. (5a) and (5b) where only the points that belong to the 95% C.L. region of the global profile-likelihood are shown in the ($m_{\chi}$, $\sigma v_{0}$) plane.

In Fig. 5a we rescale $\sigma v_{0}$ by $(\Omega_{\chi}h^{2}/0.12)^{2}$ which corresponds to the case when neutralino DM can provide only a partial contribution to the total $\Omega_{\textrm{DM}}h^{2}$. Similarly, we rescale the DM DD cross section $\sigma^{\rm SI}_{p}$ by $\Omega_{\chi}h^{2}/0.12$ when imposing the corresponding constraints. As can be seen in the plot, underabundant higgsino-like and wino-like neutralinos with masses of order few hundred GeV are typically beyond the reach of CTA. There are, however, some higgsino-like points that can be probed by the CTA monochromatic photon ($\sigma v_{\gamma\gamma}+\frac{1}{2}\sigma v_{\gamma Z}$) search even though these points lie below the projected CTA sensitivity in searches for DM-induced diffuse photon spectrum (dash-dotted line in the plot). These points are denoted by gray triangles. The crucial impact of the monochromatic line search is even more pronounced for heavier neutralinos with masses $m_{\chi}\approx 1\,\mathrm{TeV}$. In particular, it is worth stressing that, e.g., an underabundant wino-like neutralino DM with $m_{\chi}\approx 1\,\mathrm{TeV}$ can be discovered by CTA in monochromatic-line searches with no corresponding signal in the diffuse spectrum searches. For even heavier, but still underabundant, wino-like neutralinos, CTA can provide a good way of indirectly detecting them in both types of searches.

In Fig. 5b we show the results that correspond to a scenario with the neutralino being the only DM particle and having its production in the early Universe supplemented by, e.g., some non-thermal contribution. Notably, this allows one to consider neutralino DM with significantly larger values of the annihilation cross section and, therefore, much better prospects for discovery in future indirect searches. In particular, in this scenario CTA could easily discover higgsino-like neutralino DM with the mass of order a few hundred GeV in both diffuse photon and monochromatic-line searches. As can be seen in the plot, the Fermi-LAT limitsFermi-LAT:2016uux bite into the low mass region of the parameter space, where IACTs lose sensitivity. This is illustrated in Fig. 5b by a dashed line for fixed annihilation final state into a $b\bar{b}$ pair, which well represents the position of the exclusion bound we would obtain when imposing Fermi-LAT as a constraint in the likelihood. This scenario is also independently constrained by DD searches of DM, which are taken into account in our scanning procedure.

4.4 Study of benchmark points

In the previous section we have computed the H.E.S.S. limits and CTA sensitivity in the p9MSSM by combining the bounds shown in Fig. 1a weighted by the branching fractions to the appropriate final states. We have already noted, however, that, in principle, a more robust procedure should be applied. It would involve summing over all weighted spectra of annihilation final states and then using up-to-date instrument response functions and background estimates to obtain the limit as described in detail in Sec. 2. The full procedure, on the other hand, has the disadvantage of being extremely time and CPU consuming. In this section, we test the simplified treatment against the more accurate one for some carefully selected representative benchmark scenarios.

For that purpose, we choose 7 benchmark points with different neutralino properties and diverse annihilation final states. The physical properties of these points are summarized in Table 3. Among these points, BM5$-$BM7 fail to provide the thermally produced relic abundance consistent with $\Omega_{\chi}h^{2}\approx 0.12$ in the standard freeze-out scenario, but could do this, e.g., in modifed cosmological scenarios. The last two rows of the table show the difference between the 95% C.L. CTA projected sensitivities computed with the simplified and full procedures.

The dependence of this difference on final states and specific branching ratios is shown in Fig. 6. We find good agreement, better than 10%, for typical points corresponding to higgsino-like, wino-like, mixed bino-wino and some pure bino-like neutralinos. The biggest discrepancy (up to $25\%$) occurs for BM6 and BM7 which are bino-like neutralinos that annihilate primarily to leptons (note the different shape of the limit for the $\tau^{+}\tau^{-}$ annihilation final state in Fig. 1a) but that also exhibit a significant branching fraction into hadronic final states. Such points are not found in Fig. 2$-$Fig. 5b as their thermal relic density would overclose the Universe. Moreover, their $\sigma v_{0}$ is orders of magnitude below the CTA sensitivity, hence they would be irrelevant for determining CTA prospects of detecting neutralino DM within the p9MSSM, even if their relic density was altered in the desired way by assuming a modified cosmological history.

However, it is interesting to note that this discrepancy is not due to the low statistics of the signal coming from neutralinos with a small value of the annihilation cross section. It actually persists if one multiplies $\sigma v_{0}$ by an appropriate factor that brings $\sigma v_{0}$ close to the projected CTA sensitivity reach. Therefore, it could potentially also affect some analyses performed for other models of new physics in which a particle DM candidate features mixed leptonic-hadronic final annihilation states, and could lead to a sizable discrepancy between the true reach of indirect detection experiments and the one determined by the simplified approach (or similar).

5 Conclusions

In this work we performed an updated and improved study of the reach of CTA in testing neutralino DM in minimal supersymmetric scenarios. The results were compared with the most recent bounds on $\sigma v_{0}$, as a function of DM mass, obtained by H.E.S.S. We conducted the analysis in the framework of the 9-parameter MSSM, or p9MSSM. We included the most recent constraints from DM direct detection searches, flavor physics, and Higgs searches, and constructed a state-of-the-art likelihood ratio test statistic approach to analyze the CTA sensitivity. The direct constraints on sparticle masses from the LHC are also included, although they are known to be of very limited impact for the parameter space leading to TeV-scale DM. Furthermore, on the theoretical side we refined the calculations of DM relic abundance and present-day annihilation cross section by taking into account the Sommerfeld enhancement effect for a completely generic mixed neutralino and its co-annihilation partners. In particular, for the first time sfermion co-annihilations were considered with Sommerfeld effect included in a scanning framework.

Having all these improvements implemented, we performed numerical scans of the p9MSSM parameter space focusing on a TeV scale neutralino DM. We find that, assuming the Einasto profile of DM halo in the Milky Way, H.E.S.S. has been able to nearly reach the so-called thermal WIMP value, while CTA will go below it by providing a further improvement of at least an order of magnitude. The results show that both H.E.S.S. and CTA are sensitive to several cases for which direct detection cross section will be below the so-called neutrino floor, with H.E.S.S. being sensitive to most of the wino region, while CTA also covering a large fraction of the 1$\,\mathrm{TeV}$ higgsino region. We additionally show the extent to which the CTA sensitivity will be further improved in the monochromatic photon search mode for both single-component and underabundant DM.

While we focused on the Einasto profile when presenting the results for the p9MSSM, we also studied two other DM profiles, namely the standard NFW profile and the version of the Einasto profile with a core with conservative radius $r_{c}=3$ kpc, for which we presented the most up-to-date CTA sensitivities in searches relevant for a number of fixed annihilation final states. These can be easily combined to derive actual results for any model of new physics predicting heavy WIMP DM. In particular, when applied to the p9MSSM, the aforementioned Cored Einasto profile leads to substantially weaker current bounds and future sensitivity reaches. In this case, the H.E.S.S. limits do not completely exclude the region of the parameter space with wino-like neutralino DM. Instead, CTA will be able to fully probe this important scenario.

Acknowledgements.

We would like to thank Alexander Pukhov for providing a refined micrOMEGAs module that calculates annihilation cross section to monochromatic photons. We would like to thank Martin Vollmann for useful discussions. This research has made use of the CTA instrument response functions provided by the CTA Consortium and Observatory, see http://www.cta-observatory.org/science/cta-performance/ (version prod3b-v1) for more details. K.J., L.Rosz. and S.T. are supported by the National Science Centre (NCN) research grant No. 2015/18/A/ST2/00748. L.Rosz. is also supported by the National Science Centre research grant No. 2016/22/M/ST9/00583 and by “AstroCeNT: Particle Astrophysics Science and Technology Centre” project that is carried out within the International Research Agendas programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund. E.M.S. is supported in part by the National Science Centre (Poland) under the research Grant No. 2017/26/D/ST2/00490. S.T. is supported in part by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC grant: ST/P000800/1. The use of the CIS computer cluster at the National Centre for Nuclear Research in Warsaw is gratefully acknowledged.