Prospects for dark matter searches in the pMSSM

We investigate the prospects for detection of neutralino dark matter in the 19-parameter phenomenological MSSM (pMSSM). We explore very wide ranges of the pMSSM parameters but pay particular attention to the higgsino-like neutralino at the ~ 1 TeV scale, which has been shown to be a well motivated solution in many constrained supersymmetric models, as well as to a wino-dominated solution with the mass in the range of 2-3 TeV. After summarising the present bounds on the parameter space from direct and indirect detection experiments, we focus on prospects for detection of the Cherenkov Telescope Array (CTA). To this end, we derive a realistic assessment of the sensitivity of CTA to photon fluxes from dark matter annihilation by means of a binned likelihood analysis for the Einasto and Navarro-Frenk-White halo profiles. We use the most up to date instrument response functions and background simulation model provided by the CTA Collaboration. We find that, with 500 hours of observation, under the Einasto profile CTA is bound to exclude at the 95% C.L. almost all of the ~ 1 TeV higgsino region of the pMSSM, effectively closing the window for heavy supersymmetric dark matter in many realistic models. CTA will be able to probe the vast majority of cases corresponding to a spin-independent scattering cross section below the reach of 1-tonne underground detector searches for dark matter, in fact even well below the irreducible neutrino background for direct detection. On the other hand, many points lying beyond the sensitivity of CTA will be within the reach of 1-tonne detectors, and some within collider reach. Altogether, CTA will provide a highly sensitive way of searching for dark matter that will be partially overlapping and partially complementary with 1-tonne detector and collider searches, thus being instrumental to effectively explore the nearly full parameter space of the pMSSM.


Introduction
The search for particles that comprise the dark matter (DM) in the Universe has in recent years made much progress. Alternative and complementary experimental strategies are employed ranging from direct detection of DM-nucleon scattering in underground laboratories, to indirect detection of DM through observation of the Standard Model (SM) products of annihilation in astrophysical phenomena (for a recent review of the large number of experiments dedicated to direct and indirect detection of DM see, e.g., the updated version of [1] in [2] and References therein), to DM direct production at colliders [3,4].
The most impressive advances in sensitivity have arguably been made in direct detection experiments, where improvement happened rapidly and led to the most recent null results by XENON100 [5] and LUX [6]. As a consequence, the upper bounds on the spin-independent DMnucleon elastic scattering cross section, σ SI p , have become increasingly constraining for many models of weakly interacting massive particles (WIMPs). On the other hand, interesting upper bounds on the cross section for WIMP production [7,8,9,10] have been placed at the LHC, which have become particularly constraining for many models of low-mass DM. At the same time, strong limits on the present-day DM annihilation cross section, σv, as a function of the WIMP mass, have been provided by γ-ray experiments. In particular, the most stringent ones for masses up to ∼ 1 TeV come from Fermi-LAT's data on dwarf Spheroidal Galaxies (dSphs) [11]. For larger masses the air Cherenkov radiation telescope H.E.S.S. produces the strongest limits from observation of the Galactic Center (GC) [12]. The strongest indirect limits on the spin-dependent DM-proton cross section, σ SD p , have instead been obtained at IceCube/DeepCore [13,14] and ANTARES [15] in observations of neutrinos from the Sun, and in monojet searches at the LHC [9] for some choices of interactions and mediator masses.
From a theoretical perspective, the most interesting solution to the DM puzzle arguably still comes from low scale supersymmetry (SUSY), as SUSY has the ability to solve many long-standing theoretical issues within one and the same elegant framework. In common scenarios where the lightest SUSY particle (LSP) is the lightest neutralino of the Minimal Supersymmetric Standard Model (MSSM) and makes up all of the DM in the Universe many of the experiments mentioned above have already started to exclude important parts of the parameter space. Some of the solutions [16,17,18] previously favoured by considerations of electro-weak (EW) naturalness, featuring the neutralino as a somewhat balanced admixture of higgsino and gaugino quantum states with a mass m χ in the range 80 to 200 GeV now show significant tension [19,20] with the limits from XENON100 and LUX. The same mixed solutions also start to show some tension [19] with the limits on σ SD p from IceCube and ANTARES. On the other hand, under the assumption that the lightest chargino and second lightest neutralino are not much heavier than the lightest neutralino, the regions of SUSY parameter space characterised by bino-like neutralinos with mass m χ 100 − 300 GeV (depending on whether light sleptons are present) are starting to be probed [19,21,22,23] by EW-ino searches at the LHC [24,25,26,27].
At the same time, it has been shown [28,29,30,31,20,32,33] in global fits of the Constrained MSSM (CMSSM) [34] and Non-Universal Higgs Model (NUHM) that one of the consequences of the discovery at the LHC of a ∼ 125 GeV Higgs boson [35,36] in agreement with the SM is that the now favoured parameter space quite naturally gives rise to DM candidates heavier than previously thought. Typically, the ∼ 1 TeV almost pure higgsino, whose existence as a solution for DM in the MSSM has been long known [37,38] but in the framework of unified SUSY was first pointed out in a pre-LHC study of the NUHM [39], is now substantially favoured in the above mentioned CMSSM and NUHM and in some non-universal scenarios characterised by reduced EW fine tuning [40,41,42]. On the other hand, in "split" SUSY scenarios with anomaly mediation [43,44], which also rose in popularity [45,46,47,48] after the Higgs discovery, the wino is the DM. To be a thermal relic that satisfies the relic density its mass should be even larger, in the range 2-3 TeV.
Incidentally, because of its relatively large cross sections for annihilations to SM particles, the wino LSP features excellent indirect detection prospects, which have been recently investigated in several papers [49,50,51]. The annihilation cross section is enhanced in this case by inclusion of the Sommerfeld enhancement [52,53], a nonperturbative effect which affects σv and also modifies the expectations for the relic density [54,55,56]. It is particularly substantial for the thermal wino, which has now consequently been excluded at the 95% C.L. by the absence of specific signatures in existing γ-ray observatories, particularly a monochromatic line in H.E.S.S. data [57]. (The limit is relaxed if one assumes a flat halo profile [51].) On the other hand, no such statement can be made for the other TeV-scale DM candidates of the MSSM, like the ∼ 1 TeV almost pure higgsino for which the Sommerfeld enhancement is much less important. Additionally, TeV-scale neutralinos lie probably outside of the direct reach of the LHC, and existing indirect detection experiments will not reach enough sensivity to test them either. The prospects for direct detection at 1-tonne detectors are, on the other hand, very good [19,20] but upon hypothetical direct detection of a 1 TeV higgsino complementary detection by some other means will be necessary to specify its properties.
Such complementarity will likely be provided by the Cherenkov Telescope Array (CTA) [58]. The CTA project will build the next generation air Cherenkov telescope observatory. Several sensitivity studies [59,60,61,62] have shown that for DM masses greater than ∼ 100 GeV CTA is expected to significantly exceed current limits for WIMP annihilation from the Cherenkov imaging telescopes H.E.S.S. [63], MAGIC [64], and VERITAS [65], and those from Fermi-LAT [11]. CTA may even probe cross sections below the "canonical" thermal relic value of 2.6 × 10 −26 cm 3 /s for some final states.
In a previous study [20] we showed that CTA has the potential to probe the ∼ 1 TeV higgsino region of the CMSSM and NUHM parameter space. As mentioned above, this is the region that our Bayesian analysis found to be favoured by the constraints in those models, encompassing approximately 70% of the 2σ credible region in the CMSSM and 90% if it in the NUHM. In this paper we extend the study of direct and indirect DM detection prospects to the more general 19 dimensional low-scale parametrisation called p19MSSM, or more commonly pMSSM [66]. This will allow us to cover a greater number of possibilities for heavy SUSY DM: nearly pure higgsinos, winos, and bino/higgsino/wino admixtures. We will pay particular attention to neutralinos around the TeV scale, which seem to be favoured after the discovery of the Higgs boson, and to the sensitivity of CTA, which is naturally poised to probe that particular region of the parameter space. We also review the present status of direct and indirect detection constraints on the pMSSM and compare the reach of CTA with that of 1-tonne direct detection experiments and of other detection methods.
The sensitivity of CTA to the pMSSM and its complementarity to other DM searches has also been recently analysed in Ref. [67]. Some of the conclusions of this paper overlap with that study, but we also present several elements not included in [67]: • We use the most up to date [68] instrument response functions and background estimates provided by the CTA Collaboration [69].
• The sensitivity of CTA is calculated from a binned likelihood function defined on the signal and background regions, similarly to what was done in Ref. [62]. (The details of our calculation are presented at the end of this paper, in Appendix A.) • We present results for CTA sensitivity under both the Einasto [70] and Navarro-Frenk-White (NFW) [71] DM halo profiles.
• Our analysis takes into account the Sommerfeld enhancement and the consequent limits on wino DM.  Table 1: Prior ranges for the pMSSM parameters, over which we perform our scans. All masses and trilinear couplings are given in TeV. * In order to avoid generating a large number of points strongly disfavoured by the LHC we impose an additional cutoff on the physical gluino mass, mg ≥ 0.75 TeV.
The paper is organised as follows: in Sec. 2 we specify the parameter and prior ranges of our scans and their distributions. We present there the set of experimental constraints applied to the likelihood. In Sec. 3 we describe the DM properties of the parameter space regions favoured by the experimental constraints, the dominant annihilation mechanisms, and we identify a few benchmark points for future detection of DM. In Sec. 4 we present a summary of the present status of direct and indirect bounds on neutralino DM and we expose future prospects for detection, particularly at CTA, for which we accurately calculate the sensitivity reach; we also present a comparison with present and future complementary experiments. We finally give our conclusions in Sec. 5. The details of our calculation of CTA's sensitivity including a treatment of alternative statistical approaches are given in Appendix A.

Scanning methodology and experimental constraints
The pMSSM with 19 free parameters gives a generic coverage of the properties of the CP and R parity conserving MSSM. The parameters are defined at the scale of the geometrical average of the physical stop masses, M SUSY = (mt 1 mt 2 ) 1/2 , and we scan them in the ranges given in Table 1. In addition, we scan over the top quark pole mass, M t , treated here as a nuisance parameter. We assume a Gaussian distribution for M t , whose central value and experimental error are given in [72], M t = 173. 34 [41,20]. See [41,20]. See [41,20]. [ [86] is used to calculate BR B → X s γ , BR (B s → µ + µ − ), BR (B u → τ ν), and δ (g − 2) µ . M W , sin 2 θ eff , ∆M Bs are calculated using FeynHiggs. Dark matter observables Ω χ h 2 , σ SI p , σ SD p , and σv, the annihilation branching ratios, the photon fluxes for CTA and Fermi-LAT, the neutrino-induced muon flux for IceCube, and the positron flux are computed using micrOMEGAs v.3.5.5 [87].
The scans are subject to a set of constraints, applied through a global likelihood function L. The list of constraints, central values, theoretical and experimental uncertainties are presented in Table 2. We assume Gaussian distributions for the constraints, with the exception of those on the Higgs sector, which are imposed through HiggsSignals and HiggsBounds, and the constraints on σ SI p from LUX [6]. The LUX bound, which slightly improved on the limit from XENON100 [5], is included in the likelihood function following the procedure described in detail in [88,19,41,20]. Additionally, we impose 95% C.L. lower bounds from direct searches at LEP [2], smeared with 5% theoretical errors. The limits are given in Eq. (2) of Ref. [19], with the exception of the limit on the neutralino mass that has been replaced here by the LEP limit on the invisible Z width, Γ(Z → χχ) [2]. The points are gathered through several scans with both linear and log priors for the mass parameters.
We do not directly impose bounds on sparticle masses from direct SUSY searches at the LHC. As was explained in Sec. 1 we are particularly interested here on neutralinos predominantly in the range of a few hundred GeV to a few TeV, which are outside the reach of the LHC. The implementation of LHC searches in the likelihood function is for the pMSSM a numerically intensive task [19,95] that goes beyond our purpose here. We will explicitly mention in the text any situation in which a potential conflict with the limits from the LHC arises.
Although some of the scans are also driven by the constraint on the anomalous magnetic moment of the muon, δ (g − 2) µ [96,97], we present our results irrespective of whether the measurement of δ (g − 2) µ is satisfied, as the result is still somewhat controversial and favours the low-mass neutralino region, which is starting to show some tension with data from the LHC.
When appropriate, we present results for two cases. In the first, Ω χ h 2 0.12, we consider a Gaussian distribution for the relic density of DM, with experimental and theoretical uncertainties added in quadrature. In the second case, Ω χ h 2 0.12, the relic density is imposed as an upper bound only, by means of a half Gaussian likelihood.

Neutralino properties and benchmark points
In this and the following sections we show in the plots only the points the satisfy the constraints of Table 2 at the 95% C.L., i.e., we select ∆χ 2 ≤ 5.99 from the current best-fit point, where In Fig. 1(a) we show the distribution of our scan points in the (m χ , σv) plane for the case where the LSP saturates the relic abundance, Ω χ h 2 0.12 . We remind the reader that σv = σv | p→0 . The colour code gives the composition of the lightest neutralino. The equivalent distribution in the (m χ , σ SI p ) plane is shown in Fig. 1(b). The LUX bound on σ SI p is included in the likelihood; as a consequence in Fig. 1(b) almost no points lie above the 90% C.L. limit, shown here with a dashed red line for clarity.
As is well known, the neutralino mass and composition are determined by the relic density because it is a strong constraint with a relatively small uncertainty. The points of the elongated, almost vertical branches at m χ < 100 GeV belong to the Z-and h-resonance "regions" [98]. The neutralino mass is approximately half the mass of the Z boson or of the lightest Higgs, so that resonant annihilation in the early Universe leads to the correct relic density. The neutralino is predominantly bino-like with a small admixture of higgsino that does not exceed ∼ 40%. Because of their relatively low mass, neutralinos in this region will be the first among the SUSY particles to be tested at the LHC 14 TeV run, possibly even through direct DM production as in the monojet/monophoton searches, which provide limits that do not depend on the presence of light charginos or sleptons in the spectrum. On the other hand, because of their suppressed present-day annihilation cross section these points are in principle not very interesting for indirect detection.
As the neutralino mass increases, m χ > 100 GeV, predominantly bino-like LSPs (in green) satisfy the relic density through different well-understood mechanisms. From the left to the right, the point models are characterised by "bulk-like" annihilation to leptons [102,103], slepton/neutralino co-annihilation [104], or resonance with heavy A/H Higgs bosons [103]. 1 In the range 0.1 TeV m χ 1 TeV, neutralinos characterised by a mixed bino/higgsino composition satisfy the relic abundance partially through annihilation to gauge bosons, but also through off-shell Z exchange into top quarks, which can become dominant above the threshold since the other fermionic final states are suppressed by helicity conservation. The WIMPs in this region generally feature a large spin-independent cross section, so that if they have not been excluded by XENON100/LUX they are situated just below the 90% C.L. bound in Fig. 1(b). Some of the points surviving the bound show significant higgsino levels but a relative sign difference between µ and M 1 suppresses the coupling of the neutralino to the lightest Higgs boson. Additionally, cancellations between the heavy and light Higgs diagrams reduce σ SI p giving rise to known "blind spots" [105,106,107,108,109]. The position and depth of these spots depend on the relative sign of µ and M 1 and generally require that the mass of the heavy Higgs bosons are not much above the TeV scale. Note also that in the pMSSM there exists the additional possibility of augmenting neutralino annihilations in the early Universe with co-annihilations with a squark or gluino. This allows combinations of higgsino and bino components that can effectively evade the direct detection experiments. One can see some examples of this in the mixed points at m χ ∼ > 1 TeV. For indirect detection this scenario can feature a relatively large annihilation cross section but if co-annihilations are the dominant mechanism then this can be reduced below the projected sensitivity of CTA.
For masses equal and slightly larger than 1 TeV, almost pure higgsinos naturally satisfy the relic density and direct detection constraints simultaneously. Figure 1(a) shows that for these points (in red) the present-day annihilation cross section is large enough to be of interest to CTA. Pure higgsinos can lead to the correct relic density by quite a diversified set of annihilation mechanisms in the early Universe. Annihilation to gauge/Higgs bosons and chargino co-annihilation are the most natural options for masses on and above 1 TeV, but one can see in Fig. 1(a) some other higgsino points at 2-2.5 TeV, for which the above-mentioned coannihilations with stops or gluinos effectively decrease the relic abundance. Additionally, one can see a large set of higgsino points above the "canonical" thermal relic value for the cross section. There, σv is enhanced by a large resonant effect with the heavy Higgs bosons A/H that was thermally washed out in the early Universe (for an explanation see, e.g., Appendix B of Ref. [19]).
An almost pure wino LSP (in blue) can naturally satisfy the relic density for a quite wide range of mass values, 1.6 TeV m χ 2.8 TeV. Early Universe annihilation and co-annihilation of the winos characterised by the largest masses in this range receive a boost from the Sommerfeld enhancement [54,55,56]. We incorporated this effect as a correction factor to Ω χ h 2 , which we extrapolated from Fig. 5 of Ref. [56]. We also included the Sommerfeld correction to the present- day annihilation cross section by means of the rescaling factors to the different final states derived in [110]. The Sommerfeld enhancement gives rise to the peculiar "Eiffel Tower" shaped resonance that can be seen in Fig. 1(a) for wino and mixed wino/higgsino LSPs. Additionally, the spectrum also features an enhanced χχ → γγ cross section. Because of this, scenarios dominated by wino DM are in conflict with data from the current generation of γ-ray telescopes, namely they are excluded [49,50,51] at the 95% C.L. by the H.E.S.S. search for γ-ray lines [57], with the exception of cases where the halo profile is very flat. We will show in Sec. 4.1 that the same result applies to the wino-like points in our scans.
The bar chart in Fig. 2 shows the number of points with ∆χ 2 ≤ 5.99 organised by their predominant annihilation channel. We classify each point in terms of which final state has the largest branching fraction but for the majority of points annihilation will occur to several different final states. One can see that the largest number predominantly annihilates into W + W − , which is the preferred final state for both the higgsino and wino high-density regions in Fig. 1.
In Fig. 3(a) we show the distribution in the (m χ , σv) plane of the points in the first column of Fig. 2 (largest branching ratio to W + W − ). The colour code shows the actual value of the branching ratio for each point. The distributions of the points in the 3rd (bb), 4th (tt), and 5th (τ + τ − ) columns of Fig. 2 are shown in Figs. 3(b), 3(c), and 3(d), respectively. To allow a direct comparison, in the panels of Fig. 3 we indicate with a dashed black line the expected reach of CTA in the relative final state. We will discuss in detail in Sec. 4.2 and Appendix A the sensitivity of CTA. In Fig. 3 we also mark with orange stars the positions of a few benchmark points for indirect detection, which we present in Table 3 and whose characteristics we describe below.
In Fig. 3(a) one can see that the ∼ 1 TeV higgsino region must feature very substantial branching fractions to subdominant final states other than W + W − . In general, higgsinos can often have an almost equal branching fraction to ZZ, Zh or ZA, through t-channel exchange of the neutral, almost degenerate second-lightest neutralino. In Fig. 3(b) one can see that points annihilating into bb final states are very common and spread over the entire range of the parameter space. These are Figure 3: (a) Distribution in the (mχ, σv) of the points with ∆χ 2 ≤ 5.99 (see Table 2) annihilating primarily into typical of the A-resonance mechanism, when m A ≈ 2m χ , and they are characterised by moderateto-large tan β. Note that, for tan β 8, points characterised by resonance with the heavy Higgses prefer to annihilate to tt instead. One can see these points in Fig. 3(c) for masses m χ ∼ > 500 GeV (in the same part of the figure many points are also characterised by stop co-annihilation in the early Universe). The points at m χ < 500 GeV in Fig. 3(c) Table 3: Benchmark points for indirect detection. SUSY masses are given in GeV. "Decoupled" stands for any mass value above 5000 GeV. density through selectron (or smuon) co-annihilation, or annihilation to tops through an off-shell Z, for which s-wave production of light fermions is helicity-suppressed. In Fig. 3(d) one can see that many of the points that satisfy the relic density at the Z/hresonance, m χ 45 − 65 GeV, very predominantly annihilate in the present-day Universe to τ + τ − and their cross section is depleted with respect to the "canonical" thermal value. This is because, on the one hand the Boltzman distribution for WIMPs in this range is sharply peaked at momentum p > 0 so that one expects a drop in the cross section now, when p → 0; on the other hand, the lightest stau happens to be for these points light enough (mτ 1 100 − 500 GeV) to allow significant annihilation to taus through stau exchange in the t channel. Note that even if the points in this region are predominantly bino-like they need a non-negligible amount of mixing with the higgsino to efficiently annihilate in the early Universe. Thus, for many of these points the lightest chargino, χ ± 1 , and the second-lightest neutralino,χ 0 2 , have masses not much far above m χ , and they are therefore very likely to be already excluded by 3-lepton searches at the LHC [24,25,26,27]. However, there are several other points in this region for which m χ ± 1 ∼ > 450 GeV so that they are not currently excluded and remain viable scenarios. The 14 TeV run at the LHC should definitively probe these points.
We present four benchmark points for DM detection in Table 3, each one is marked by an orange star in the panels of Fig. 3. They are chosen such that they are not presently excluded but represent scenarios within the projected reach of CTA. Benchmark point 1 (BM1) belongs to the ∼ 1 TeV higgsino region. The dominant final states are W + W − and ZZ. The γγ and Zγ cross sections are suppressed compared to the continuous spectrum but remain relatively large due to the higgsino nature of the neutralino and the associated light chargino that appears in the loops. Benchmark point 2 (BM2) is a typical A-resonance point annihilating predominately to bb. It is the only benchmark point without sizeable higgsino-bino mixing. Due to this the γγ and Zγ cross sections are heavily suppressed. Benchmark point 3 (BM3) is a mixed bino-higgsino neutralino. The tt final state is produced non-resonantly and thus requires large mixing to increase the coupling of the neutralino to the Z boson. The chargino associated with the higgsino component is relatively light and induces the sub-dominant annihilation channels to W + W − and ZZ. Finally, benchmark point 4 (BM4) is another mixed bino-higgsino neturalino similar to BM3 although its dominant annihilation channel is now mediated by a light stau leading to a τ + τ − final state. For BM4 the neutralinos in the early Universe also annihilate into top quarks at threshold, but top production is suppressed for p → 0 because the neutralino is slightly too light.
It should be noted that both mixed bino-higgsino scenarios, BM3 and BM4, have a spinindependent scattering cross section close to the current limit from LUX, thus making these scenarios highly visible in future direct detection experiments. BM2 may be detectable at future 1-tonne direct detection experiments, while BM1 lies in a direct detection blind spot.

Prospects for dark matter detection
We present here the prospects for DM detection in the pMSSM. In Sec. 4.1 we summarise the current status of indirect-detection bounds on the parameter space of the pMSSM. We show the case where the neutralino saturates the relic abundance, Ω χ h 2 0.12, and the case Ω χ h 2 0.12. In Sec. 4.2 we present our estimate of the sensitivity reach of CTA. We give results for the Einasto and NFW DM halo profiles, obtained through the procedure described in detail in Appendix A.

Current indirect detection bounds on the pMSSM
It was shown, e.g., in Ref. [19] for the p9MSSM and in Ref. [67] for the p19MSSM that the bounds from the diffuse γ-ray spectrum from dSphs at Fermi-LAT [11] are too weak to affect the MSSM parameter space when the PLANCK value [89] on Ω χ h 2 is imposed. The limits can however become important in scenarios where the relic abundance is satisfied as an upper bound only. This opens up large regions of the parameter space, characterised by lighter higgsino and wino LSPs that annihilate away effectively in the early Universe and have a much suppressed relic abundance. These solutions could correspond to a scenario with two-component DM, in which case the sensitivity of Fermi-LAT should be rescaled according to the square of the density, R 2 , to account for the reduced presentday density of neutralinos (R = Ω χ /Ω Planck ). Alternatively, one can assume that the neutralino represents the entirety of the DM and the correct abundance is fixed by invoking some additional mechanism [111]. In this case one obtains much stronger limits from indirect detection than in the previous scenario because rescaling by R 2 is not necessary [112]. Only in this case does the limit from dSphs have significant impact on the parameter space. However, the limit from γ-ray searches at H.E.S.S. [57] can also affect a small part of the parameter space when one applies rescaling.
In Fig. 4(a) we show in red the points excluded at the 95% C.L. by the H.E.S.S. search under the assumption of rescaling by R 2 . In Fig. 4(b) we show the points excluded at the 95% C.L. by Fermi-LAT dSphs (in orange) and the above mentioned H.E.S.S. γ-ray line search (in red) projected to the (m χ , σv) plane in the case of no rescaling. In the absence of additional guidance from theory, the real limit is likely to lie somewhat in between the exclusions of Figs. 4(a) and 4(b).
Going back to the case where the neutralino saturates the relic density, the strongest indirect limits on the spin-dependent scattering cross section come from IceCube/DeepCore [13,14] and ANTARES [15], from observation of neutrinos from the Sun. It is well known [113] that σ SD p can easily be measured through the relation between the DM solar capture rate and the annihilation rate, which should give rise to a high-energy neutrino spectrum and flux.
In Fig. 5(a) we show the current limits from IceCube and ANTARES, provided for 100% branching ratios to W + W − , bb, and τ + τ − , compared to the points of the pMSSM. The colour code shows the dominant annihilation final state for each point, in the sense of Fig. 2.
As was mentioned in Sec. 3 the γ-ray line search at H.E.S.S. [57] is very effective in placing a strong bound on wino DM when Ω χ h 2 0.12. We show in Fig. 5(b) the current limit on the cross section under the assumption of the Einasto halo profile, compared to the pMSSM points. As was previously shown in [49,50,51], when one includes the Sommerfeld enhancement H.E.S.S. excludes at the 95% C.L. pure or almost pure wino DM, with the exception of cases where a flat DM profile is assumed. (We include the Sommerfeld enhancement here as a correction factor derived from the calculations of Ref. [110].) For the remainder of this paper we will not show in the plots the points excluded by H.E.S.S.
In Fig. 6 we show the DM annihilation contribution to the positron fraction for the benchmark points of Table 3. We compare our benchmark points to the measured positron fraction at Pamela [114] and AMS-02 [115]. We use micrOMEGAs to calculate the produced positron spectrum for each point and the positron flux at the Earth after propagation. We use the default values in micrOMEGAs for charged particle propagation. To compare with the positron fraction we use the parametrisation of the primary electron background and secondary electron and positron fluxes  Figure 5: (a) The current 90% C.L. limits on σ SD p from IceCube [13,14] and ANTARES [15]. The limits are presented for the W + W − , bb, and τ + τ − final states. Points are shown with ∆χ 2 ≤ 5.99 (see Table 2). The colour code shows the primary annihilation final state for each point. (b) The current 95% C.L. limits on σγγv from H.E.S.S. [57]. In red we show the points that are excluded by the experimental bound, indicated with a solid black line. Points are shown with ∆χ 2 ≤ 5.99 (see Table 2).
following [116,117]. The fluxes in GeV cm −2 s −1 sr −1 are given by, Figure 6: The DM annihilation contribution to the positron fraction for the benchmark points of Table 3 compared to the measured positron fraction at Pamela [114] and AMS-02 [115].  Table 2). (3) For all of the benchmarks the positron flux is many orders of magnitude too small to explain the anomalous positron fraction, which is also true for the rest of the points.

Sensitivity of CTA to the pMSSM
The sensitivity of CTA is obtained by testing the background-only hypotesis with a likelihood function. The construction of the likelihood is detailed in Appendix A. We use the most up to date [68] instrument response functions and background estimates provided by the CTA Collaboration [69]. We assume an observation time of approximately 500 hours, and we provide results under the Einasto and NFW DM halo profiles, for which we calculate the J-factors.
We start with the case where the lightest neutralino makes up all of the DM, Ω χ h 2 0.12. In Fig. 7(a) we show the 95% C.L. sensitivity of CTA in the (m χ , σv) plane. We show points within the reach of CTA under the assumption of the NFW profile in red, while the points excluded under the Einasto profile but not the NFW are shown in orange. Green points lie beyond the reach of CTA for either profile. We repeat that we do not show in the plots the points excluded by H.E.S.S. (see Fig. 5(b)), which feature a nearly pure wino LSP or an admixture of wino and higgsino. Our improved analysis shows that the reach of CTA for SUSY is even more optimistic than our previous study had anticipated [20] and under the Einasto profile CTA is bound to exclude at the 2σ level the ∼ 1 TeV higgsino region of the pMSSM in almost its entirety, effectively closing the window for heavy SUSY DM in many realistic models.  Figure 8: The reach of CTA in the pMSSM for the case Ωχh 2 0.12. The colour code is the same as in Fig. 7(a).
(a) The 95% C.L. reach of CTA in the (mχ, R 2 · σv) plane for the case with rescaling of the local density. (b) The 95% C.L. reach of CTA in the (mχ, σv) plane for the case without rescaling.
The point is emphasised in Fig. 7(b), where we show the boost factor, b F , required for each point to be within the 95% C.L. sensitivity of CTA for the Einasto profile. Points below b F = 1 can be excluded without a boost factor. One can see that, in agreement with Fig. 7(a), the majority of points in the ∼ 1 TeV higgsino region do not require any boost factor.
Incidentally, we want to underline here the power of the binned likelihood function especially in relation to the points that show a significant cross section to γγ. This is what happens for the points in red and orange at m χ < 300 GeV in Fig. 7(a), which are projected to be excluded under one or the other profile assumption. Their γ-ray flux effectively gives rise to an excess in one or a few particular bins, large enough to kill the likelihood function irrespectively of all other bins.
We now move on to describe the reach of CTA in the case where the neutralino does not saturate the relic density, Ω χ h 2 0.12. Figure 8(a) shows the expected reach of CTA in the (m χ , R 2 · σv) plane in the case of rescaling σv by the ratio (squared) of the neutralino density to the PLANCK value. We do not show in the plot the points excluded by Fermi-LAT dSphs and the γ-ray line search at H.E.S.S., see Fig. 4. One can see that CTA will probe regions of the parameter space beyond the limit presently set by H.E.S.S., especially if one assumes the Einasto halo profile. A few points (in red) are excluded by CTA also under assumption of the NFW profile but in that case much of the parameter space lies beyond the reach of CTA due to the rescaling factor R 2 .
As was the case in Fig. 4, when one assumes no rescaling the sensitivity of CTA can exclude a large fraction of the parameter space under both profile assumptions, as one can see in Fig. 8(b) where the CTA reach is shown in the (m χ , σv) plane.

Complementarity of CTA with other experiments
In Fig. 9(a) we show the projected reach of CTA in the (m χ , σ SI p ) plane to compare it with direct detection experiments, which are sensitive to the spin-independent neutralino-proton cross section, σ SI p . We use the same colour code as in Fig. 7(a) for the sensitivity of CTA. The projected sensitivity (a) (b) Figure 9: In all panels the colour code is the same as in Fig. 7(a). Points are shown with ∆χ 2 ≤ 5.99 (see Table 2). of XENON-1T [118], indicative of the generic reach of 1-tonne detectors like, e.g., DEAP-3600 [119] and LZ [120], is shown as a dashed grey line. The dot-dashed magenta line shows the onset of the irreducible background due to atmospheric and diffuse supernova neutrinos [99,100,101]. Since there is considerable overlap between the regions within the reach of CTA and those out of reach, it should be noted that the points are overlaid from the most constrained to the least constrained. Thus, the plots represent a conservative estimate of the reach of CTA in these planes. Figure 9(a) shows that over most of the pMSSM parameter space the reach of CTA is orthogonal to that of detectors directly measuring σ SI p . CTA will be able to probe the vast majority of the points that lie well beyond the reach of 1-tonne detectors and even reach the region where the irreducible neutrino background will strongly curb sensitivity advances for direct detection experiments.
In Fig. 9(b) we show the sensitivity of CTA compared to the projected 5-year sensitivity of the IceCube/DeepCore 86-string configuration, shown as a light shading in the plot. For clarity, we also separate the shaded points from the normally coloured points with a dashed grey line. The sensitivity has been here separately estimated for the IceCube upward, IceCube contained, and DeepCore contained events, following the procedure described in detail, e.g., in Refs. [121,122]. We take the effective area for IceCube and the effctive volume for DeepCore as parametrised in [123,121,124]. The light shaded points in Fig. 9(b) fall within the reach of at least one of these three estimates.
To highlight the idea of complementarity, we show in Fig. 10(a) the reach of CTA compared to the one of 1-tonne detectors in the (σ SI p , σv) plane. The projected sensitivity of CTA appears as horizontal bands that follow the same colour code as in Fig. 7(a), while the points within the projected reach of 1-tonne detectors are shown as lighter shaded. The dashed grey line vertically  Figure 10: In all panels the colour code is the same as in Fig. 7(a). Points are shown with ∆χ 2 ≤ 5.99 (see Table 2). (a) The sensitivity of CTA to the pMSSM points in the (σ SI p , σv) plane. Lighter shaded points are within the projected sensitivity of 1-tonne detectors. The dashed grey line gives an approximate reference value for future direct detection reach in σ SI p . (b) The sensitivity of CTA to the pMSSM points in the (σ SD p , σv) plane. Lighter shaded points are within the projected 5-year sensitivity of IceCube/DeepCore. The dashed grey line gives an approximate reference value for future reach in σ SD p .
separates the points within the sensitivity of 1-tonne detectors (lighter shading) from the ones below sensitivity (regular colouring). In Fig. 10(b) we present the equivalent picture in the (σ SD p , σv) plane, where the light-shaded region indicates the points within 5-year sensitivity at IceCube/DeepCore, and the dashed grey line separates the points within sensitivity from those outside the reach.
As was mentioned in Secs. 2 and 3, a detailed analysis of the LHC reach in the pMSSM is beyond the scope of this paper because it interests regions of the parameter space characterised by a low mass neutralino, orthogonal in some sense to the regions in which CTA is most sensitive. However, for completeness we show in Fig. 11 the reach of CTA compared to the present limits on sparticle masses obtained in simplified models at the LHC. These are meant to be seen as a generic indication of the present reach of the LHC and the reader must keep in mind that in complex models a detailed analysis can often produce limits that can be either weaker or stronger than the ones under simplified model assumptions [95].
In Fig. 11(a) we show the reach of CTA in the (mg, m χ ) plane. The solid black line shows the 95% C.L. bound obtained at ATLAS for the simplified model of [125]. The limit obtained at CMS is very similar [128]. The reach of CTA in the (mt 1 , m χ ), (mb 1 , m χ ), and (m χ ± 1 , m χ ) planes is shown in Figs. 11(b), 11(c), and 11(d), respectively. With the exception of the coannihilation bands, shown in Fig. 11 as regions of mass degeneracy (e.g. m χ ≈ mg in Fig. 11(a)), the reach of CTA is largely independent of the sparticle spectrum, as was to be expected. Improvements in the limits on the gluino and squark masses are not expected to have any effect on the sensitivity of CTA. Indeed, CTA remains sensitive to spectra where the gluinos and squarks lie well beyond the reach of present and future colliders.
Note that in Fig. 11(d) we have indicated with a solid black line the bound from 3 lepton EW-ino searches obtained under the assumption that a light slepton of the first two generations is present Figure 11: In all panels the colour code is the same as in Fig. 7(a). (a) Sensitivity of CTA in the (mg, mχ) plane, the thick black line shows the current ATLAS limit in the simplified model of Ref. [125]. (b) Sensitivity of CTA in the (mt 1 , mχ) plane, the thick black line shows the current ATLAS limit in the simplified model of Ref. [126]. (c) Sensitivity of CTA in the (mb 1 , mχ) plane, the thick black line shows the current ATLAS limit in the simplified model of Ref. [127]. (d) Sensitivity of CTA in the (m χ ± 1 , mχ) plane, the thick black line shows the current ATLAS limit in the simplified model of Ref. [25]. and has a mass in between m χ and m χ ± 1 or mχ0 2 [25,26]. As it stands the limits seems to exclude almost entirely the Z/h-resonance region of the parameter space, shown in the lower left-hand side of the plot, at m χ < 100 GeV. One must remember that for the points of the pMSSM characterised by heavier sleptons the limit is actually weaker than shown here [25,27], thus allowing part of the Z/h-resonance region to survive the bounds.

Summary and conclusions
In this paper we have investigated the prospects for detection of neutralino DM in the framework of the pMSSM. We have focused here particularly on models for which the LSP is a neutralino at about the TeV scale, outside the reach of the 14 TeV run at LHC. TeV-scale neutralinos have been shown to be promising candidates for DM in many well-motivated SUSY scenarios after the Higgs boson discovery.
In order to satisfy the constraints from the relic density neutralinos in this mass range must preferentially be nearly pure higgsinos, or winos, or a mixture of the two. However, as is known from many studies recently appeared in the literature, thermal-relic winos are subject to large Sommerfeld enhancement of their annihilation cross section, so that they are excluded at the 2σ level by the H.E.S.S. search for γ-ray lines from the GC under most choices of the DM halo profile. Consequently the most likely candidate for heavy SUSY DM that is not excluded by the present constraints is the ∼ 1 TeV nearly pure higgsino.
We showed in this paper that the air Cherenkov radiation telescope array of imminent construction CTA has the sensitivity reach to nearly exclude at the 95% C.L. the ∼ 1 TeV higgsino region of the pMSSM with 500 hours of observation. We derived the sensitivity of CTA to diffuse and γ-ray line photons arising from neutralino DM annihilation by constructing an energy-binned likelihood function and adopted the most up to date estimates of the detector response functions and modelling of the background provided by the CTA collaboration. We obtained results for the Einasto and NFW profile assumptions. We applied our results to the parameter space of the pMSSM, but also presented the limits for single annihilation final state channels.
We showed that under the NFW profile assumption CTA will be able to probe 70% of the points belonging to the ∼ 1 TeV higgsino region in our scans, specifically those for which the annihilation cross section is enhanced by one of the following factors: resonance with a heavy Higgs boson, A/H; non-negligible admixture with the wino; or a significant annihilation to a monochromatic γγ line. Assumption of the Einasto profile allows us to formulate even more optimistic projections, as we found that in this case CTA will also probe the remaining part of the ∼ 1 TeV higgsino region, with the exclusion of some points characterised by heavy squark or gluino co-annihilation, and will probe additionally a significant part of the parameter space leading to bino and mixed bino/higgsino DM.
We also found that, in complementarity with other direct and indirect detection experiments, CTA will significantly probe the favoured parameter region of the pMSSM, far beyond the reach of 1-tonne underground detectors alone. We showed that many of the points well within our calculated sensitivity of CTA lie below the onset of the irreducible atmospheric and diffuse supernova neutrino background for direct detection.
Finally we found that by combining different experiments detection prospects for SUSY DM are very good also in scenarios where the neutralino is not the only particle comprising the DM in the Universe. The existing bounds and projections strongly depend in this case on the adopted prescription for rescaling the local DM abundance. We found that even under the most conservative assumption, that the local density is rescaled with the square of the ratio of the neutralino density to the PLANCK value, CTA will be able to exclude a large region of the parameter space, significantly outperforming the reach of current Fermi-LAT dSphs and H.E.S.S. γ-ray line searches.
Alexander Pukhov for clarifying some issues regarding micrOMEGAs. This work has been funded in part by the Welcome Programme of the Foundation for Polish Science. L.R. is also supported in part by a STFC consortium grant of Lancaster, Manchester, and Sheffield Universities. The use of the CIS computer cluster at the National Centre for Nuclear Research is gratefully acknowledged. L.R. would like to thank the Theory Division of CERN for hospitality during the final stages of the project.

A Derivation of the CTA reach
In this appendix we present the calculation of CTA's sensitivity to the WIMP annihilation cross section. We compare the sensitivities that can be obtained through different likelihood functions, and we refer to some other recent work on the subject that can be found in the literature [59,60,61,62].
We use the Ring Method as defined in [59,60,61,62]. Two regions are identified in the plane of the galactic coordinates l and b, as shown in Fig. 12. The "signal", or ON, region is based on a circle of angular radius ∆ cut = 1.36 • around the GC. The "background", or OFF, region is based on a ring centred at the offset coordinate b off = 1.42 • , with an inner angular radius of r 1 = 0.55 • and an outer radius of r 2 = 2.88 • , from which the ON region is subtracted. The strip of sky characterised by |b| < 0.3 • about the GC and the region of the sky within the inner radius r 1 do not belong to either the ON or OFF regions.
We proceed to create a binned likelihood function. For each energy bin, i, the expected number of counts from DM annihilation is calculated: where A eff is the effective area of the detector, δ(E) is the energy resolution, dN γ /dE is the annihilation spectrum, and J is the J-factor for either the ON or OFF region. For A eff and δ(E) we use the most up to date instrument response functions provided by the CTA Collaboration [68]. We take an observation time t obs = 500 h . As is well known, the J-factor is defined by integration along the line of sight (l.o.s.) of the DM halo profile squared, ρ 2 (r): where in the reference frame centred at the observer the distance from the GC depends on the azimuthal angle θ. Using the Einasto DM density profile [70], we obtain J Ein ON = 7.44 × 10 21 GeV 2 /cm 5 and J Ein OFF = 1.21 × 10 22 GeV 2 /cm 5 for the J-factors of the ON and OFF regions, respectively, in agreement with the results of [62]. The corresponding values for the NFW [71] profile are J NFW ON = 3.89 × 10 21 GeV 2 /cm 5 and J NFW OFF = 5.78 × 10 21 GeV 2 /cm 5 . To gauge the constraining power of the binned likelihood we compare the projected sensitivity for DM annihilating to bb for three different methods.
1. We first calculate the binned Poisson likelihood defined by In Eq. (6) µ ij is the expected number of photons in bin ij, with the index i running over the energy bins and j = 1, 2 for the ON and OFF region, respectively. It is given by the sum of the background Figure 12: In blue, the ON region: the angular radius is ∆cut = 1.36 • . In pink, the OFF region: the offset from the GC is b off = 1.42 • ; the inner radius is r1 = 0.55 • and the outer radius is r2 = 2.88 • . The strip of sky at |b| < 0.3 • is subtracted from the ON and OFF regions. and the expected count from DM annihilation, Eq. (4). n ij is the number of photons counted in bin ij. We find the 95% C.L. limit by setting n ij equal to the number of background-only photons and calculating L by increasing the annihilation cross section from the best fit value (σv= 0) until the difference in −2 ln L from the best fit value is 2.71 (one-sided 95% C.L.). We use an updated [68] background estimate provided by the CTA Collaboration [69].
2. In the limit of a single energy bin one can construct a test statistics for the background-only hypothesis [129] − 2 ln L = 2 N ON ln 1 + α α where N ON and N OFF are the number of counts in the ON and OFF regions, respectively, and α = ∆Ω ON /∆Ω OFF is the ratio of the ON and OFF solid angles. We calculate the projected sensitivity by setting N ON and N OFF equal to the expected number of counts N j = i µ ij . We calculate α = 0.2457. We will compare the sensitivity obtained with Eq. (7), hereafter referred to as the Li and Ma method, with the binned likelihood method, Eq. (6). 3. We finally compare the binned Poisson likelihood, Eq. (6), to the the binned Skellam distribution [130], on which the results of [61] are based. Note that the results of Ref. [61] were used in our previous study [20] to estimate the sensitivity of CTA to the parameter space of the CMSSM and the NUHM.
If one calculates for each energy bin i the difference in the measured number of events θ diff,i = n i,ON − α n i,OFF , the probability of observing a set of particular values {θ diff,i } is given by In Fig. 13 we show the limits obtained for the bb final state and Einasto profile using these three methods. For low masses the limits obtained using the Skellam distribution and the Li and Ma method are comparable. For larger masses the Li and Ma method loses sensitivity compared to the other two methods as they benefit from the data being binned by energy. The binned likelihood method produces the strongest limits and we adopt this method for producing our projected sensitivities. We find limits that are stronger than the estimates of [61] and [62] but comparable to [60].
It should be noted that we have not included uncertainties in the DM distribution, systematic uncertainties in the detector response, or included the diffuse γ-ray background. The uncertainty in the DM distribution enters into the calculation of the limits through the J-factors in Eq. (4). We account for this by presenting limits for two different realisations of the DM distribution, the NFW and Einasto profiles. The systematic uncertainty due to the finite energy resolution of the experiment already appears in Eq. (4). However further systematic effects can be present such as varying acceptance across the field of view or uncertainties in the effective area.
Finally, the diffuse astrophysical γ-ray background around the galactic centre measured by H.E.S.S. [131] and Fermi-LAT [132] presents a challenge to this type of ON/OFF analysis since this background will be larger in the ON (signal) region than the OFF (background) region mimicking the searched for signal. Thus the sensitivity presented here is somewhat optimistic and would be reduced due to the diffuse background and systematic uncertainties, however the treatment of Ref. [62] suggests that a morphological analysis could partially mitigate these effects.  Figure 14 shows the derived 95% C.L. limits for some specific final states including the most common primary annihilation channels found in the pMSSM using the binned likelihood of Eq. (6). The limits obtained can probe values of σv below the "canonical" thermal relic value for all of the final states. We compare these to the actual final states found by the scan in Fig. 3, Sec. 3.