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

###### 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 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.

## 1 Introduction

Dark matter (DM) is the dominant component of matter in the universe but its nature is still 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 seeks 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 its direct detection 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 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 can 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 background. 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 H.E.S.S., MAGIC and VERITAS experiments are performing 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 few years.

The GC region is arguably the most likely astrophysical environment to detect DM annihilation signals in very high energy (VHE, ) 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 papers Fowlie: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 CTA Acharya: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 former works by: (i) obtaining 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 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 higgsino region (see Krall:2017xij (); Kowalska:2018toh () for a recent work and review) that otherwise would remain unexplored. Finally, we present our conclusions in Sec. 5.

## 2 Dark matter 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 ( 100 GeV) gamma rays. Among them are H.E.S.S. J1745-290 Aharonian:2009zk (), a strong emission spatially coincident with the supermassive black hole Sagittarius A*, the supernova/pulsar wind nebula G09+01 Aharonian:2005br (), the supernova remnant H.E.S.S. J1745-303 Aharonian:2008gw (), and a diffuse emission extending along the Galactic plane Aharonian: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 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 continuum Abdallah:2016ygi () and mono-energetic gamma line Abdallah: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 in Galactic latitude along the Galactic plane has been excluded from the dataset together with a disk of radius centered at the position of J1745-303. No excess in the signal region with respect to background was found and constraints among the strongest for TeV DM have been derived in various annihilation channels Abdallah: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 a DM search with the planned CTA Acharya: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 searches Moulin: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 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 (), 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 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 developed Aharonian:2008zza (). However, due to the 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 simulations HASSAN201776 (). Here, we use the so-called prod3b version of the instrument response functions including the residual background determinations.

The region of interest for the DM signal extraction with CTA extends up to from the GC both in Galactic longitude and latitude. The overall region is split into squared pixels of side . 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. 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 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 between in both Galactic longitude and latitudes. 20 energy bins are taken logarithmically-spaced between energies from 10 to 100, 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 -th energy bin, -th Galactic longitude bin, and -th Galactic latitude bin. It reads

(1) |

where the likelihood function follows a Poisson distribution given by while 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 and , 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, is taken to be 1. In the context of a counting experiment using ON-OFF measurements Abdallah:2018qtu (), the signal is searched in an ON region where the measured number of events is . The expected number of background events in the signal region, , is determined from the measurement of number of events in control (OFF) regions, , with no or little expected searched signal. 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. .

The total likelihood function is the product of the individual likelihood functions over each bin defined as . The log-likelihood ratio test statistics (LLRTS) is defined as:

(2) |

where the single and double carets indicate unconditional and conditional maximization, respectively Cowan: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 . The expected sensitivity is computed for a 100% branching ratio in each of the channels , , , , , , , , , and .

### 2.4 Expected signal and background events

The expected photon flux from pair-annihilation of DM particles of mass in a region of solid angle in the sky can be expressed as

(3) |

where is the total annihilation cross section to all primary channels providing photons in the final sates, and is the photon spectrum per annihilation. is the so-called J-factor defined as the integral of the square of the DM density along the line of sight and over by

(4) |

is the distance along the line of sight from the observer and is related to the radial distance in the coordinates centered at GC by , where is the angle between the direction of observation and the GC, and = 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 Einasto 1965TrAlm…5…87E () and Navarro-Frenk-White (NFW) Navarro:1996gj () profiles defined as

(5) |

where normalization , scale radius , and power index are given in Table 1. The local DM density is taken to be Catena:2009mf (). Since the DM density in the GC is far from being certain we also consider a Cored Einasto profile with a core radius such that and . 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.

Profiles | Einasto (E) | NFW | Cored Einasto (CE) |
---|---|---|---|

(GeVcm) | 0.079 | 0.307 | 0.079 |

(kpc) | 20.0 | 21.0 | 20.0 |

0.17 | 0.17 | ||

(kpc) | 3.0 |

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

(6) |

where is the gamma-ray energy-dependent effective area, is the observation time and is a Gaussian function centered at the DM mass 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 is taken from Ref. Cirelli:2010xx () for continuum channels. The monoenergetic gamma-ray line is a Dirac delta function centered at . The signal count rate in the -th bin is integrated over the spatial pixel of solid angle size and energy bin width .

The background count number is given by

(7) |

where is the energy-differential residual background rate per steradian. A careful modelling 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. A band of 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. (a)a ^{1}^{1}1All sensitivity limits, including more annihilation channels, for all three halo profiles considered, can be found in the supplementary material on the arXiv. Limits provided there also extend to 100.
the projected CTA 95% C.L. sensitivity to DM annihilation as a function of DM mass . 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 three different choices of the DM Galactic halo profile.
The CTA 95% C.L. sensitivity to monochromatic -ray lines for the same three choices of halo profile is featured in Fig. (b)b.

## 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 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 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)MSSM Djouadi:1998di (); Cahill-Rowley:2014twa (). The latter is based on the following assumptions 1) conservation, 2) Minimal Flavor Violation at the electroweak scale, 3) degenerate first two generations of sfermion soft-mass parameters and 4) negligible Yukawa couplings and trilinear couplings (-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 , which decouples them (see Table 2). The p9MSSM provides a sufficiently generic parametrization and coverage of the DM properties of the MSSM with 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 Feroz:2007kg (); Feroz:2008xx () package for the scan, using flat priors. The supersymmetric spectrum is calculated with SPheno v4.0.3 Porod:2003um (); Porod:2011nf (). We allow the bino mass and the parameter to assume negative values in order to accommodate blind spots in DM direct detection Ellis:2000ds (); Cheung:2012qy (), which stem from the vanishing coupling for certain combinations of parameters (see also Han:2016qtc () for a recent discussion). The remaining gaugino mass parameter is kept positive, starting from the a minimal value of 100, 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 set, while the latter, in case of squarks, can more easily lead to a correct value of the Higgs boson mass, , thanks to an increase of the characteristic SUSY scale. As discussed above, the remaining sfermion mass parameters and the gluino mass are fixed at 20. They do not play an important role in a further discussion.

Parameter | Range |
---|---|

bino mass | |

wino mass | |

gluino mass | |

trilinear couplings | |

pseudoscalar mass | |

parameter | |

3rd gen. left soft squark mass | |

3rd gen. right up soft squark mass | |

3rd gen. right down soft squark mass | |

1st/2nd gen. soft squark masses | |

soft slepton masses | |

soft slepton masses | |

ratio of Higgs doublet VEVs | |

Nuisance parameter | Central value, error |

Top pole mass (GeV) | (173.34, 0.76) ATLAS:2014wva () |

The SUSY mass parameters are defined at the scale of the geometrical average of the physical stop masses, . The ratio of the Higgs doublets’ vevs, , and the top quark pole mass, , which is treated here as a nuisance parameter, are defined at the electroweak symmetry breaking (EWSB) scale. We assume a Gaussian distribution for , whose central value and experimental error are given in ATLAS:2014wva (), . 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 Planck Aghanim:2018eyx (),

(8) |

To calculate the relic density we employ Belanger:2001fz (); Belanger:2004yn () supplemented by Hryczuk:2011tq (). We additionally impose a theoretical uncertainty on the calculation to partially take into account the effects of, e.g., loop corrections Baro:2007em (); Baro:2009na (), variations in the renormalization scheme and scale Harz:2016dql (), and modifications to the QCD equations of state Hindmarsh:2005ix (); Laine:2006cp (); Drees:2015exa ().

At the typical mass scale tested by CTA and H.E.S.S. ( to a few TeV) the SE plays an important role and strongly affects both the calculation of the present-day neutralino annihilation cross section, , and, to a lesser degree also the determination of the thermal neutralino relic density Hisano: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 v5 Gondolo:2004sc ().^{2}^{2}2DarkSE 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 matrix Beneke: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 being currently developed Beneke:inprep (), which has been already used in several phenomenological studies Beneke: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 and include the SE by rescaling the result using a grid of the enhancements in the - plane following the procedure of Roszkowski:2014iqa (); 2) the final points are then post-processed with the accurate SE treatment using full code.
Sommerfeld enhancement is also included in the computation of the present-day , as well as for and .^{3}^{3}3It has been checked that the zero-velocity limit of these cross sections gives essentially the same result as 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, 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 work Harz: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 in Ciafaloni:2010ti () and explicitly seen in the wino annihilation computation at one-loop Hryczuk:2011vi () while finally approached with resummation techniques in Baumgart:2014vma (); Bauer:2014ula (); Ovanesyan:2014fwa (). In very recent works Baumgart:2017nsr (); Baumgart:2018yed (), however, it has been shown that for the neutralino annihilation the modification due to fully resummed exclusive cross section is most relevant in multi-TeV regime. Therefore, and also due to the fact that it is currently not possible to directly apply the framework of Baumgart:2017nsr () to generic p9MSSM neutralinos, we do not include this effect in our scan.

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

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 the only one 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 Planck Aghanim: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 the extended reheating period Giudice:2000ex (); Fornengo:2002db (); Gelmini:2006pq () (see also Roszkowski:2014lga (); Drees:2018dsj () for recent studies) or in presence of additional non-thermal production. In this case, the neutralino can be the only 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 .

#### 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. These are relevant for those TeV-scale neutralinos of predominant but not very pure higgsino and wino eigenstate composition.

We include the most recent DD bounds here. We employ Arbey:2018msw () and the package Workgroup:2017lvb (), assuming the Standard Halo Model (SHM) and the following values for the relevant astrophysical parameters: , , . We note that slight modifications to the SHM that might be suggested by e.g. recent data release by the GAIA Collaboration Brown: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-2 Cui:2017nnn (), PICO-60 Amole:2017dex (), and the most recent results from the XENON1T collaboration Aprile: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 for of our interest come from the searches by the IceCube Collaboration for neutrinos coming from the center of the Milky Way Aartsen:2017ulx (), the Earth Aartsen:2016fep (), and the Sun Aartsen: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 dwarfs Krall: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 particles Tanabashi: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 Bechtle:2008jh (); Bechtle:2015pma () and Bechtle:2013xfa (), while additional constraints from searches for heavy Higgs decays to are implemented following ATLAS:2016fpj (); CMS:2016rjp ().

#### Flavor physics

We calculate a few flavor observables with Arbey:2018msw (). The parameter space of the MSSM is potentially sensitive, in particular, to the bounds from rare decays in processes and radiative decays like , which can constrain scan points characterized by large values, and/or relatively light non-SM Higgs bosons, squarks, and charginos/neutralinos. We use the following experimental determinations:

(9) | |||||

(10) |

where, following, e.g., Ref. Workgroup:2017myk (), in Eq. (9) we give the calculated average Misiak: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 collision data Aaij:2017vad (). We thus implicitly assume that Eq. (10) has superseded an older statistical combination of CMS and LHCb measurements with 7 and 8 data CMS:2014xfa (). Note that very recently the ATLAS Collaboration has presented a measurement of , from a combination of data taken during their 8 and 13 runs, which agrees with Eq. (10): Aaboud:2018mst ().

We impose the bounds of Eqs. (9) and (10) at the 95% C.L., a posteriori on the points belonging to the region of the profile likelihood. This reduces the number of viable points in the scan by approximately 2%. Other potentially relevant flavor observables like or the mass mixing measurement 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 discrepancy with the SM, like the differential branching ratios and angular observables in Aaij:2014pli (); Aaij:2015oid (), or the branching ratio measurements that have recently provided tantalizing hints of lepton flavor nonuniversality Lees: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 measurement Bennett:2006fi () of the muon anomalous magnetic moment, which shows a discrepancy with the SM expectation, 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 data Grange: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 MSSM Darme: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 -rays that correspond to 6 years of data and observation of dSphs Fermi-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 pair. We have also verified numerically, following Arbey:2018msw (), that taking into account a complete list of annihilation final states leads to similar results.

We note that 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 annihilations Padmanabhan: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 (see Beneke:2016jpw () for recent discussion).

A recent determination Cuoco:2017iax () of the bounds on the neutralino annihilation cross section from AMS-02 antiproton cosmic-ray (CR) data Aguilar:2016kjl () has proven to be competitive with H.E.S.S. diffuse -ray searches in the mass range. We have verified that this is in general true also in the context of our scans using Arbey:2018msw () that employs a semianalytic approach to solving the propagation equations following Boudaud: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 -ray signal.

## 4 Results

### 4.1 CTA sensitivity to the p9MSSM

For each point in the scan, we compute the H.E.S.S. limit for the present-day annihilation cross section, , and the corresponding projected sensitivity of CTA. We use the 95% C.L. bounds and projections for annihilation to pure channels (see Fig. (a)a). 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 final state we use limit, for final states with and quarks – limit, for the lightest quarks – limit and for we employ 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, and , we employ the bounds computed for the SM Higgs; for the charged MSSM Higgs particle, , we use the bounds derived for . 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 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 from the best-fit point, where .

### 4.2 Discussion of results

We present in Fig. 4 the scan points in the plane (, ) of the present-day annihilation cross section of the neutralino versus its mass. The color code used in Fig. 4 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, , 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., if we order gauge eigenstes as {bino, wino, down-type higgsino, up-type higgsino}, we must have for these points). In blue, the points for which it is for at least 90% a wino (). In cyan we show admixtures of these two gaugino states, with the additional constraint that the higgsino composition remain below 10%, . In red we show neutralinos that are dominated for at least 90% by their higgsino fraction (). We finally point out that only very few points characterized by admixtures of a gaugino and a large higgsino component appear, 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 scan 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 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 data Athron:2018vxy ().

We show in Fig. 4 with a solid black line the current 95% C.L. upper bound on from 254 h 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, this allows us to determine the position of the H.E.S.S. limit in the region with a light neutralino and large 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 SE Cohen: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 at 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 (nearly pure) higgsino region of the parameter space, corresponding to the red points in Fig. 4, 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.’s at 1 TeV, that will be probing 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 h 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 -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 analyses Roszkowski: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 TeV the approximate limit that we present reproduces very well the true CTA sensitivity. For lower masses, the line shown in Fig. 4 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. 4.

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 including one or more Higgs boson of mass within a few hundred GeV of , which provide a means for funnels, or resonant -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 the cases 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, or , respectively, quite independently of the model specifics of the rest of the sparticle spectrum.

The projected sensitivity of CTA shown in Fig. 4 is obtained in the two limiting cases of Einasto and Cored Einasto DM halo profiles. In order to infer the corresponding sensitivity in the case of the NFW profile, it is an easy task to multiply the projected line by rescaling factor obtained from Fig. (a)a, which reads approximately 1:2.5. In Fig. 4, 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. 5 we show the p9MSSM points in the (, ) plane. The most recent XENON1T 90% C.L. upper limit Aprile: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 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. 4 and we additionally overlap violet triangles to points excluded by the H.E.S.S. bound on . 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 LHC Roszkowski:2014iqa (). We show the power of complementarity of direct and indirect detection in Fig. 6, where we project the points of the p9MSSM to the (, ) plane. The color code is the same as in Fig. 4 and Fig. 5.

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 characteristic of the neutrino “floor” for direct DM searches depends on the DM mass, hence the onset of the shaded area is jagged in Fig. 6). Conversely, the sensitivity of IACTs gradually improves from the top down, providing a perpendicular 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. To guide the eye, we add a vertical dashed gray line denoting the neutrino background limit taken at .

### 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. ((a)a) and ((b)b) where only the points that belong to the 95% C.L. region of the global profile-likelihood are shown in the (, ) plane.

In Fig. (a)a we rescale by which corresponds to the case when neutralino DM can provide only a partial contribution to the total . Similarly, we rescale the DM DD cross section by 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 () 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 . In particular, it is worth stressing that, e.g., an underabundant wino-like neutralino DM with 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. (b)b 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 limits Fermi-LAT:2016uux () bite into the low mass region of the parameter space, where IACTs lose sensitivity. This is illustrated in Fig. (b)b by a dashed line for fixed annihilation final state into a 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. (a)a 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 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, BM5BM7 fail to provide the thermally produced relic abundance consistent with 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 diffence 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. 10. 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 ) occurs for BM6 and BM7 which are bino-like neutralinos that annihilate primarily to leptons (note the different shape of the limit for the annihilation final state in Fig. (a)a) but that also exhibit a significant branching fraction into hadronic final states. Such points are not found in Fig. 4Fig. (b)b as their thermal relic density would overclose the Universe. Moreover, their is by 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 by an appropriate factor that brings 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).

Benchmark | BM1 | BM2 | BM3 | BM4 | BM5 | BM6 | BM7 |

points | |||||||

1099 | 1765 | 1840 | 531 | 1516 | 2288 | 997 | |

Branching | 0.23 | 0.35 | 0.64 | 0.85 | 0.17 | 0.26 | 0.39 |

fractions | 0.23 | 0.29 | 0.14 | 0.14 | 0.16 | 0.22 | 0.37 |

0.21 | ZZ 0.24 | ZH 0.14 | 0.01 | 0.16 | 0.21 | 0.22 | |

ZZ 0.20 | 0.05 | 0.08 | ZH 0.16 | 0.14 | 0.01 | ||

Zh 0.06 | 0.04 | 0.16 | 0.13 | ||||

0.04 | Zh 0.03 | 0.157 | 0.03 | ||||

0.10 | 0.16 | 0.13 | 0.13 | - | - | - | |

Main | , | , | slepton | ||||

mechanism | coann. | -funnel | coann. | t-channel | coann. | t-channel | coann. |

## 5 Conclusions

In this work we have performed an updated and improved study of the reach of CTA in testing neutralino DM in minimal supersymmetric scenarios. The results are compared with the most recent bounds on , as a function of DM mass, obtained by H.E.S.S. We perform the analysis in the framework of the 9-parameter MSSM, or p9MSSM. We have 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, sfermion co-annihilations for the first time 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 higgsino region. We additionally show to what extent the CTA sensitivity will be further improved in the monochromatic photon search mode for both single-component and underabundant DM.

While we focus 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 radius kpc, for which we have 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 predicitng 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.

We would like to thank Alexander Pukhov for providing a refined micrOMEGAs module that calculates annihilation cross section to monochromatic photons. 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. and L.Ro. are supported by the National Science Centre (NCN) research grant No. 2015/18/A/ST2/00748. L.Ro. is also supported by the National Science Centre research grant No. 2016/22/M/ST9/00583 and by the project âAstroCeNT: Particle Astrophysics Science and Technology Centreâ carried out within the International Research Agendas programme of the Foundation for Polish Science 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/L000520/1. The use of the CIS computer cluster at the National Centre for Nuclear Research in Warsaw is gratefully acknowledged.

## References

- (1) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, The waning of the WIMP? A review of models, searches, and constraints, Eur. Phys. J. C78 (2018), no. 3 203, [arXiv:1703.07364].
- (2) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, WIMP dark matter candidates and searches – current status and future prospects, Rept. Prog. Phys. 81 (2018), no. 6 066201, [arXiv:1707.06277].
- (3) A. Fowlie, K. Kowalska, L. Roszkowski, E. M. Sessolo, and Y.-L. S. Tsai, Dark matter and collider signatures of the MSSM, Phys. Rev. D88 (2013) 055012, [arXiv:1306.1567].
- (4) M. Cahill-Rowley, J. L. Hewett, A. Ismail, and T. G. Rizzo, Lessons and prospects from the pMSSM after LHC Run I, Phys. Rev. D91 (2015), no. 5 055002, [arXiv:1407.4130].
- (5) L. Roszkowski, E. M. Sessolo, and A. J. Williams, Prospects for dark matter searches in the pMSSM, JHEP 02 (2015) 014, [arXiv:1411.5214].
- (6) M. E. Cabrera-Catalan, S. Ando, C. Weniger, and F. Zandanel, Indirect and direct detection prospect for TeV dark matter in the nine parameter MSSM, Phys. Rev. D92 (2015), no. 3 035018, [arXiv:1503.00599].
- (7) ATLAS Collaboration, G. Aad et al., Summary of the ATLAS experimentâs sensitivity to supersymmetry after LHC Run 1 â interpreted in the phenomenological MSSM, JHEP 10 (2015) 134, [arXiv:1508.06608].
- (8) M. Beneke, A. Bharucha, A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, The last refuge of mixed wino-Higgsino dark matter, JHEP 01 (2017) 002, [arXiv:1611.00804].
- (9) A. Arbey, M. Boudaud, F. Mahmoudi, and G. Robbins, Robustness of dark matter constraints and interplay with collider searches for New Physics, JHEP 11 (2017) 132, [arXiv:1707.00426].
- (10) GAMBIT Collaboration, P. Athron et al., A global fit of the MSSM with GAMBIT, Eur. Phys. J. C77 (2017), no. 12 879, [arXiv:1705.07917].
- (11) S. Abel, D. G. Cerdeño, and S. Robles, The Power of Genetic Algorithms: what remains of the pMSSM?, arXiv:1805.03615.
- (12) Cherenkov Telescope Array Consortium Collaboration, B. S. Acharya et al., Science with the Cherenkov Telescope Array, arXiv:1709.07997.
- (13) R. Krall and M. Reece, Last Electroweak WIMP Standing: Pseudo-Dirac Higgsino Status and Compact Stars as Future Probes, Chin. Phys. C42 (2018), no. 4 043105, [arXiv:1705.04843].
- (14) K. Kowalska and E. M. Sessolo, The discreet charm of higgsino dark matter - a pocket review, Adv. High Energy Phys. 2018 (2018) 6828560, [arXiv:1802.04097].
- (15) C. van Eldik, Gamma rays from the Galactic Centre region: a review, Astropart. Phys. 71 (2015) 45–70, [arXiv:1505.06055].
- (16) H.E.S.S. Collaboration, F. Aharonian et al., Spectrum and variability of the Galactic Center VHE gamma-ray source HESS J1745-290, Astron. Astrophys. 503 (2009) 817, [arXiv:0906.1247].
- (17) H.E.S.S. Collaboration, F. Aharonian et al., Very high energy gamma rays from the composite SNR G0.9+0.1, Astron. Astrophys. 432 (2005) L25–L29, [astro-ph/0501265].
- (18) H.E.S.S. Collaboration, F. Aharonian, Exploring a SNR/Molecular Cloud Association Within HESS J1745-303, Astron. Astrophys. 483 (2008) 509–517, [arXiv:0803.2844].
- (19) H.E.S.S. Collaboration, F. Aharonian et al., Discovery of very-high-energy gamma-rays from the galactic centre ridge, Nature 439 (2006) 695–698, [astro-ph/0603021].
- (20) H.E.S.S. Collaboration, A. Abramowski et al., Acceleration of petaelectronvolt protons in the Galactic Centre, Nature 531 (2016) 476, [arXiv:1603.07730].
- (21) H.E.S.S. Collaboration, H. Abdallah et al., Search for dark matter annihilations towards the inner Galactic halo from 10 years of observations with H.E.S.S, Phys. Rev. Lett. 117 (2016), no. 11 111301, [arXiv:1607.08142].
- (22) H.E.S.S. Collaboration, H. Abdallah et al., Search for -Ray Line Signals from Dark Matter Annihilations in the Inner Galactic Halo from 10 Years of Observations with H.E.S.S., Phys. Rev. Lett. 120 (2018), no. 20 201101, [arXiv:1805.05741].
- (23) E. Moulin, J. Carr, J. Gaskins, M. Doro, C. Farnier, M. Wood, and H. Zechlin, Dark Matter Programme, in Science with the Cherenkov Telescope Array, pp. 45–81. 2019.
- (24) T. Hassan, L. Arrabito, K. BernlÃ¶hr, J. Bregeon, J. Cortina, P. Cumani, F. D. Pierro, D. Falceta-Goncalves, R. Lang, J. Hinton, T. Jogler, G. Maier, A. Moralejo, A. Morselli, C. T. Peixoto, and M. Wood, Monte carlo performance studies for the site selection of the cherenkov telescope array, Astropart. Phys. 93 (2017) 76.
- (25) L. Rinchiuso, N. L. Rodd, I. Moult, E. Moulin, M. Baumgart, T. Cohen, T. R. Slatyer, I. W. Stewart, and V. Vaidya, Hunting for Heavy Winos in the Galactic Center, arXiv:1808.04388.
- (26) “CTA’s expected baseline performance.” https://www.cta-observatory.org/science/cta-performance/.
- (27) F. Aharonian, J. Buckley, T. Kifune, and G. Sinnis, High energy astrophysics with ground-based gamma ray detectors, Rept. Prog. Phys. 71 (2008) 096901.
- (28) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics, Eur. Phys. J. C71 (2011) 1554, [arXiv:1007.1727]. [Erratum: Eur. Phys. J.C73,2501(2013)].
- (29) J. Einasto, On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters, Trudy Astrofizicheskogo Instituta Alma-Ata 5 (1965) 87–100.
- (30) J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal density profile from hierarchical clustering, Astrophys. J. 490 (1997) 493–508, [astro-ph/9611107].
- (31) R. Catena and P. Ullio, A novel determination of the local dark matter density, JCAP 1008 (2010) 004, [arXiv:0907.0018].
- (32) M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 1103 (2011) 051, [arXiv:1012.4515]. [Erratum: JCAP1210,E01(2012)].
- (33) MSSM Working Group Collaboration, A. Djouadi et al., The Minimal supersymmetric standard model: Group summary report, in GDR (Groupement De Recherche) - Supersymetrie Montpellier, France, April 15-17, 1998, 1998. hep-ph/9901246.
- (34) F. Feroz and M. P. Hobson, Multimodal nested sampling: an efficient and robust alternative to MCMC methods for astronomical data analysis, Mon. Not. Roy. Astron. Soc. 384 (2008) 449, [arXiv:0704.3704].
- (35) F. Feroz, M. P. Hobson, and M. Bridges, MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics, Mon. Not. Roy. Astron. Soc. 398 (2009) 1601–1614, [arXiv:0809.3437].
- (36) W. Porod, SPheno, a program for calculating supersymmetric spectra, SUSY particle decays and SUSY particle production at e+ e- colliders, Comput. Phys. Commun. 153 (2003) 275–315, [hep-ph/0301101].
- (37) W. Porod and F. Staub, SPheno 3.1: Extensions including flavour, CP-phases and models beyond the MSSM, Comput. Phys. Commun. 183 (2012) 2458–2469, [arXiv:1104.1573].
- (38) J. R. Ellis, A. Ferstl, and K. A. Olive, Reevaluation of the elastic scattering of supersymmetric dark matter, Phys. Lett. B481 (2000) 304–314, [hep-ph/0001005].
- (39) C. Cheung, L. J. Hall, D. Pinner, and J. T. Ruderman, Prospects and Blind Spots for Neutralino Dark Matter, JHEP 05 (2013) 100, [arXiv:1211.4873].
- (40) T. Han, F. Kling, S. Su, and Y. Wu, Unblinding the dark matter blind spots, JHEP 02 (2017) 057, [arXiv:1612.02387].
- (41) ATLAS, CDF, CMS, D0 Collaboration, First combination of Tevatron and LHC measurements of the top-quark mass, arXiv:1403.4427.
- (42) Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, arXiv:1807.06209.
- (43) G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov, MicrOMEGAs: A Program for calculating the relic density in the MSSM, Comput. Phys. Commun. 149 (2002) 103–120, [hep-ph/0112278].
- (44) G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov, micrOMEGAs: Version 1.3, Comput. Phys. Commun. 174 (2006) 577–604, [hep-ph/0405253].
- (45) A. Hryczuk, The Sommerfeld enhancement for scalar particles and application to sfermion co-annihilation regions, Phys. Lett. B699 (2011) 271–275, [arXiv:1102.4295].
- (46) N. Baro, F. Boudjema, and A. Semenov, Full one-loop corrections to the relic density in the MSSM: A Few examples, Phys. Lett. B660 (2008) 550–560, [arXiv:0710.1821].
- (47) N. Baro, F. Boudjema, G. Chalons, and S. Hao, Relic density at one-loop with gauge boson pair production, Phys. Rev. D81 (2010) 015005, [arXiv:0910.3293].
- (48) J. Harz, B. Herrmann, M. Klasen, K. Kovarik, and P. Steppeler, Theoretical uncertainty of the supersymmetric dark matter relic density from scheme and scale variations, Phys. Rev. D93 (2016), no. 11 114023, [arXiv:1602.08103].
- (49) M. Hindmarsh and O. Philipsen, WIMP dark matter and the QCD equation of state, Phys. Rev. D71 (2005) 087302, [hep-ph/0501232].
- (50) M. Laine and Y. Schröder, Quark mass thresholds in QCD thermodynamics, Phys. Rev. D73 (2006) 085009, [hep-ph/0603048].
- (51) M. Drees, F. Hajkarim, and E. R. Schmitz, The Effects of QCD Equation of State on the Relic Density of WIMP Dark Matter, JCAP 1506 (2015), no. 06 025, [arXiv:1503.03513].
- (52) J. Hisano, S. Matsumoto, M. Nagai, O. Saito, and M. Senami, Non-perturbative effect on thermal relic abundance of dark matter, Phys. Lett. B646 (2007) 34–38, [hep-ph/0610249].
- (53) M. Cirelli, A. Strumia, and M. Tamburini, Cosmology and Astrophysics of Minimal Dark Matter, Nucl. Phys. B787 (2007) 152–175, [arXiv:0706.4071].
- (54) A. Hryczuk, R. Iengo, and P. Ullio, Relic densities including Sommerfeld enhancements in the MSSM, JHEP 03 (2011) 069, [arXiv:1010.2172].
- (55) P. Gondolo, J. Edsjö, P. Ullio, L. Bergstrom, M. Schelke, and E. A. Baltz, DarkSUSY: Computing supersymmetric dark matter properties numerically, JCAP 0407 (2004) 008, [astro-ph/0406204].
- (56) M. Beneke, C. Hellmann, and P. Ruiz-Femenia, Non-relativistic pair annihilation of nearly mass degenerate neutralinos and charginos I. General framework and S-wave annihilation, JHEP 03 (2013) 148, [arXiv:1210.7928]. [Erratum: JHEP10,224(2013)].
- (57) M. Beneke, C. Hellmann, and P. Ruiz-Femenia, Non-relativistic pair annihilation of nearly mass degenerate neutralinos and charginos III. Computation of the Sommerfeld enhancements, JHEP 05 (2015) 115, [arXiv:1411.6924].
- (58) M. Beneke, C. Hellmann, and P. Ruiz-Femenia, Heavy neutralino relic abundance with Sommerfeld enhancements - a study of pMSSM scenarios, JHEP 03 (2015) 162, [arXiv:1411.6930].
- (59) M. Beneke, A. Bharucha, F. Dighera, C. Hellmann, A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, In preparation, .
- (60) M. Beneke, A. Bharucha, F. Dighera, C. Hellmann, A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, Relic density of wino-like dark matter in the MSSM, JHEP 03 (2016) 119, [arXiv:1601.04718].
- (61) T. R. Slatyer, The Sommerfeld enhancement for dark matter with an excited state, JCAP 1002 (2010) 028, [arXiv:0910.5713].
- (62) S. El Hedri, A. Kaminska, and M. de Vries, A Sommerfeld Toolbox for Colored Dark Sectors, Eur. Phys. J. C77 (2017), no. 9 622, [arXiv:1612.02825].
- (63) J. Harz and K. Petraki, Radiative bound-state formation in unbroken perturbative non-Abelian theories and implications for dark matter, JHEP 07 (2018) 096, [arXiv:1805.01200].
- (64) P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia, and A. Urbano, Weak Corrections are Relevant for Dark Matter Indirect Detection, JCAP 1103 (2011) 019, [arXiv:1009.0224].
- (65) A. Hryczuk and R. Iengo, The one-loop and Sommerfeld electroweak corrections to the Wino dark matter annihilation, JHEP 01 (2012) 163, [arXiv:1111.2916]. [Erratum: JHEP06,137(2012)].
- (66) M. Baumgart, I. Z. Rothstein, and V. Vaidya, Calculating the Annihilation Rate of Weakly Interacting Massive Particles, Phys. Rev. Lett. 114 (2015) 211301, [arXiv:1409.4415].
- (67) M. Bauer, T. Cohen, R. J. Hill, and M. P. Solon, Soft Collinear Effective Theory for Heavy WIMP Annihilation, JHEP 01 (2015) 099, [arXiv:1409.7392].
- (68) G. Ovanesyan, T. R. Slatyer, and I. W. Stewart, Heavy Dark Matter Annihilation from Effective Field Theory, Phys. Rev. Lett. 114 (2015), no. 21 211302, [arXiv:1409.8294].
- (69) M. Baumgart, T. Cohen, I. Moult, N. L. Rodd, T. R. Slatyer, M. P. Solon, I. W. Stewart, and V. Vaidya, Resummed Photon Spectra for WIMP Annihilation, JHEP 03 (2018) 117, [arXiv:1712.07656].
- (70) M. Baumgart, T. Cohen, E. Moulin, I. Moult, L. Rinchiuso, N. L. Rodd, T. R. Slatyer, I. W. Stewart, and V. Vaidya, Precision Photon Spectra for Wino Annihilation, arXiv:1808.08956.
- (71) G. F. Giudice, E. W. Kolb, and A. Riotto, Largest temperature of the radiation era and its cosmological implications, Phys. Rev. D64 (2001) 023508, [hep-ph/0005123].
- (72) N. Fornengo, A. Riotto, and S. Scopel, Supersymmetric dark matter and the reheating temperature of the universe, Phys. Rev. D67 (2003) 023514, [hep-ph/0208072].
- (73) G. Gelmini, P. Gondolo, A. Soldatenko, and C. E. Yaguna, The Effect of a late decaying scalar on the neutralino relic density, Phys. Rev. D74 (2006) 083514, [hep-ph/0605016].
- (74) L. Roszkowski, S. Trojanowski, and K. Turzyński, Neutralino and gravitino dark matter with low reheating temperature, JHEP 11 (2014) 146, [arXiv:1406.0012].
- (75) M. Drees and F. Hajkarim, Neutralino Dark Matter in Scenarios with Early Matter Domination, arXiv:1808.05706.
- (76) A. Arbey, F. Mahmoudi, and G. Robbins, SuperIso Relic v4: A program for calculating dark matter and flavour physics observables in Supersymmetry, arXiv:1806.11489.
- (77) The GAMBIT Dark Matter Workgroup Collaboration, T. Bringmann et al., DarkBit: A GAMBIT module for computing dark matter observables and likelihoods, Eur. Phys. J. C77 (2017), no. 12 831, [arXiv:1705.07920].
- (78) Gaia Collaboration, A. G. A. Brown et al., Gaia Data Release 2, Astron. Astrophys. 616 (2018) A1, [arXiv:1804.09365].
- (79) L. Necib, M. Lisanti, and V. Belokurov, Dark Matter in Disequilibrium: The Local Velocity Distribution from SDSS-Gaia, arXiv:1807.02519.
- (80) N. W. Evans, C. A. J. O’Hare, and C. McCabe, Refinement of the standard halo model for dark matter searches in light of the Gaia Sausage, Phys. Rev. D99 (2019), no. 2 023012, [arXiv:1810.11468].
- (81) PandaX-II Collaboration, X. Cui et al., Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment, Phys. Rev. Lett. 119 (2017), no. 18 181302, [arXiv:1708.06917].
- (82) PICO Collaboration, C. Amole et al., Dark Matter Search Results from the PICO-60 CF Bubble Chamber, Phys. Rev. Lett. 118 (2017), no. 25 251301, [arXiv:1702.07666].
- (83) XENON Collaboration, E. Aprile et al., Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121 (2018), no. 11 111302, [arXiv:1805.12562].
- (84) IceCube Collaboration, M. G. Aartsen et al., Search for Neutrinos from Dark Matter Self-Annihilations in the center of the Milky Way with 3 years of IceCube/DeepCore, Eur. Phys. J. C77 (2017), no. 9 627, [arXiv:1705.08103].
- (85) IceCube Collaboration, M. G. Aartsen et al., First search for dark matter annihilations in the Earth with the IceCube Detector, Eur. Phys. J. C77 (2017), no. 2 82, [arXiv:1609.01492].
- (86) IceCube Collaboration, M. G. Aartsen et al., Search for annihilating dark matter in the Sun with 3 years of IceCube data, Eur. Phys. J. C77 (2017), no. 3 146, [arXiv:1612.05949].
- (87) M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, Dark Kinetic Heating of Neutron Stars and An Infrared Window On WIMPs, SIMPs, and Pure Higgsinos, Phys. Rev. Lett. 119 (2017), no. 13 131801, [arXiv:1704.01577].
- (88) K. Kowalska, Phenomenological MSSM in light of new 13 TeV LHC data, Eur. Phys. J. C76 (2016), no. 12 684, [arXiv:1608.02489].
- (89) GAMBIT Collaboration, P. Athron et al., Combined collider constraints on neutralinos and charginos, arXiv:1809.02097.
- (90) Particle Data Group Collaboration, M. Tanabashi et al., Review of Particle Physics, Phys. Rev. D98 (2018), no. 3 030001.
- (91) P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein, and K. E. Williams, HiggsBounds: Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron, Comput. Phys. Commun. 181 (2010) 138–167, [arXiv:0811.4169].
- (92) P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak, and G. Weiglein, Applying Exclusion Likelihoods from LHC Searches to Extended Higgs Sectors, Eur. Phys. J. C75 (2015), no. 9 421, [arXiv:1507.06706].
- (93) P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak, and G. Weiglein, : Confronting arbitrary Higgs sectors with measurements at the Tevatron and the LHC, Eur. Phys. J. C74 (2014), no. 2 2711, [arXiv:1305.1933].
- (94) ATLAS Collaboration, T. A. collaboration, Search for Minimal Supersymmetric Standard Model Higgs Bosons in the final state in up to 13.3 fb of pp collisions at = 13 TeV with the ATLAS Detector, .
- (95) CMS Collaboration, C. Collaboration, Search for a neutral MSSM Higgs boson decaying into with of data at , .
- (96) The GAMBIT Flavour Workgroup Collaboration, F. U. Bernlochner et al., FlavBit: A GAMBIT module for computing flavour observables and likelihoods, Eur. Phys. J. C77 (2017), no. 11 786, [arXiv:1705.07933].
- (97) M. Misiak and M. Steinhauser, Weak radiative decays of the B meson and bounds on in the Two-Higgs-Doublet Model, Eur. Phys. J. C77 (2017), no. 3 201, [arXiv:1702.04571].
- (98) BaBar Collaboration, B. Aubert et al., Measurement of the branching fraction and photon energy spectrum using the recoil method, Phys. Rev. D77 (2008) 051103, [arXiv:0711.4889].
- (99) BaBar Collaboration, J. P. Lees et al., Exclusive Measurements of Transition Rate and Photon Energy Spectrum, Phys. Rev. D86 (2012) 052012, [arXiv:1207.2520].
- (100) BaBar Collaboration, J. P. Lees et al., Precision Measurement of the Photon Energy Spectrum, Branching Fraction, and Direct CP Asymmetry , Phys. Rev. Lett. 109 (2012) 191801, [arXiv:1207.2690].
- (101) Belle Collaboration, T. Saito et al., Measurement of the Branching Fraction with a Sum of Exclusive Decays, Phys. Rev. D91 (2015), no. 5 052004, [arXiv:1411.7198].
- (102) Belle Collaboration, A. Abdesselam et al., Measurement of the inclusive branching fraction, photon energy spectrum and HQE parameters, in Proceedings, 38th International Conference on High Energy Physics (ICHEP 2016): Chicago, IL, USA, August 3-10, 2016, 2016. arXiv:1608.02344.
- (103) LHCb Collaboration, R. Aaij et al., Measurement of the branching fraction and effective lifetime and search for decays, Phys. Rev. Lett. 118 (2017), no. 19 191801, [arXiv:1703.05747].
- (104) CMS, LHCb Collaboration, V. Khachatryan et al., Observation of the rare decay from the combined analysis of CMS and LHCb data, Nature 522 (2015) 68–72, [arXiv:1411.4413].
- (105) ATLAS Collaboration, M. Aaboud et al., Study of the rare decays of and mesons into muon pairs using data collected during 2015 and 2016 with the ATLAS detector, Submitted to: JHEP (2018) [arXiv:1812.03017].
- (106) LHCb Collaboration, R. Aaij et al., Differential branching fractions and isospin asymmetries of decays, JHEP 06 (2014) 133, [arXiv:1403.8044].
- (107) LHCb Collaboration, R. Aaij et al., Angular analysis of the decay using 3 fb of integrated luminosity, JHEP 02 (2016) 104, [arXiv:1512.04442].
- (108) BaBar Collaboration, J. P. Lees et al., Measurement of an Excess of Decays and Implications for Charged Higgs Bosons, Phys. Rev. D88 (2013), no. 7 072012, [arXiv:1303.0571].
- (109) LHCb Collaboration, R. Aaij et al., Measurement of the ratio of branching fractions , Phys. Rev. Lett. 115 (2015), no. 11 111803, [arXiv:1506.08614]. [Erratum: Phys. Rev. Lett.115,no.15,159901(2015)].
- (110) Belle Collaboration, S. Hirose et al., Measurement of the lepton polarization and in the decay , Phys. Rev. Lett. 118 (2017), no. 21 211801, [arXiv:1612.00529].
- (111) LHCb Collaboration, R. Aaij et al., Test of lepton universality with decays, JHEP 08 (2017) 055, [arXiv:1705.05802].
- (112) LHCb Collaboration, R. Aaij et al., Search for lepton-universality violation in decays, arXiv:1903.09252.
- (113) W. Altmannshofer and D. M. Straub, New physics in transitions after LHC run 1, Eur. Phys. J. C75 (2015), no. 8 382, [arXiv:1411.3161].
- (114) Muon g-2 Collaboration, G. W. Bennett et al., Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D73 (2006) 072003, [hep-ex/0602035].
- (115) M. Davier, Update of the Hadronic Vacuum Polarisation Contribution to the muon g-2, Nucl. Part. Phys. Proc. 287-288 (2017) 70–75, [arXiv:1612.02743].
- (116) K. Kowalska, L. Roszkowski, E. M. Sessolo, and A. J. Williams, GUT-inspired SUSY and the muon anomaly: prospects for LHC 14 TeV, JHEP 06 (2015) 020, [arXiv:1503.08219].
- (117) Muon g-2 Collaboration, J. Grange et al., Muon (g-2) Technical Design Report, arXiv:1501.06858.
- (118) L. Darmé, K. Kowalska, L. Roszkowski, and E. M. Sessolo, Flavor anomalies and dark matter in SUSY with an extra U(1), JHEP 10 (2018) 052, [arXiv:1806.06036].
- (119) DES, Fermi-LAT Collaboration, A. Albert et al., Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT, Astrophys. J. 834 (2017), no. 2 110, [arXiv:1611.03184].
- (120) N. Padmanabhan and D. P. Finkbeiner, Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects, Phys. Rev. D72 (2005) 023508, [astro-ph/0503486].
- (121) S. Galli, F. Iocco, G. Bertone, and A. Melchiorri, CMB constraints on Dark Matter models with large annihilation cross-section, Phys. Rev. D80 (2009) 023505, [arXiv:0905.0003].
- (122) T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, CMB Constraints on WIMP Annihilation: Energy Absorption During the Recombination Epoch, Phys. Rev. D80 (2009) 043526, [arXiv:0906.1197].
- (123) A. Cuoco, J. Heisig, M. Korsmeier, and M. KrÃ¤mer, Constraining heavy dark matter with cosmic-ray antiprotons, JCAP 1804 (2018), no. 04 004, [arXiv:1711.05274].
- (124) AMS Collaboration, M. Aguilar et al., Antiproton Flux, Antiproton-to-Proton Flux Ratio, and Properties of Elementary Particle Fluxes in Primary Cosmic Rays Measured with the Alpha Magnetic Spectrometer on the International Space Station, Phys. Rev. Lett. 117 (2016), no. 9 091103.
- (125) M. Boudaud, M. Cirelli, G. Giesen, and P. Salati, A fussy revisitation of antiprotons as a tool for Dark Matter searches, JCAP 1505 (2015), no. 05 013, [arXiv:1412.5696].
- (126) C. Evoli, I. Cholis, D. Grasso, L. Maccione, and P. Ullio, Antiprotons from dark matter annihilation in the Galaxy: astrophysical uncertainties, Phys. Rev. D85 (2012) 123511, [arXiv:1108.0664].
- (127) M. Cirelli and G. Giesen, Antiprotons from Dark Matter: Current constraints and future sensitivities, JCAP 1304 (2013) 015, [arXiv:1301.7079].
- (128) A. Djouadi, J. Kalinowski, and M. Spira, HDECAY: A Program for Higgs boson decays in the standard model and its supersymmetric extension, Comput. Phys. Commun. 108 (1998) 56–74, [hep-ph/9704448].
- (129) A. Djouadi, J. Kalinowski, M. Muehlleitner, and M. Spira, HDECAY: Twenty years after, Comput. Phys. Commun. 238 (2019) 214–231, [arXiv:1801.09506].
- (130) T. Cohen, M. Lisanti, A. Pierce, and T. R. Slatyer, Wino Dark Matter Under Siege, JCAP 1310 (2013) 061, [arXiv:1307.4082].
- (131) A. Hryczuk, I. Cholis, R. Iengo, M. Tavakoli, and P. Ullio, Indirect Detection Analysis: Wino Dark Matter Case Study, JCAP 1407 (2014) 031, [arXiv:1401.6212].
- (132) V. Lefranc et al. Phys. Rev. D 91 (2015), no. 12 122003, [arXiv:1502.05064].