Global analysis of the pMSSM in light of the Fermi GeV excess: prospects for the LHC Run-II and astroparticle experiments

We present a new global fit of the 19-dimensional phenomenological Minimal Supersymmetric Standard Model (pMSSM-19) that comply with all the latest experimental results from dark matter indirect, direct and accelerator dark matter searches. We show that the model provides a satisfactory explanation of the excess of gamma-rays from the Galactic centre observed by the Fermi~Large Area Telescope, assuming that it is produced by the annihilation of neutralinos in the Milky Way halo. We identify two regions that pass all the constraints: the first corresponds to neutralinos with a mass ~80-100 GeV annihilating into WW with a branching ratio of 95% ; the second to heavier neutralinos, with mass ~180-200 GeV annihilating into t tbar with a branching ratio of 87%. We show that neutralinos compatible with the Galactic centre GeV excess will soon be within the reach of LHC run-II -- notably through searches for charginos and neutralinos, squarks and light smuons -- and of Xenon1T, thanks to its unprecedented sensitivity to spin-dependent cross-section off neutrons.


Introduction
The existence of dark matter (DM) is now robustly established [1][2][3][4] and its cosmological abundance measured with high precision [5].Yet, the fundamental nature of the most abundant matter component in the Universe is unknown.According to the most promising theories, DM is a new fundamental particle.As a consequence, the search for DM is also a search for new physics beyond that of the well-known elementary particles, as laid down in the Standard Model (SM).Weakly interacting massive particles (WIMPs) are the leading DM candidates: they arise in many extensions of the Standard Model of particle physics, and naturally achieve the appropriate relic density through self-annihilation the early Universe.WIMPs can be searched for with three detection strategies: direct detection of the energy recoil of nuclei scattering off DM particles; indirect detection of the final stable products of DM annihilation or decay, as for example gamma rays; and accelerator searches for new particles, in particular at the Large Hadron Collider (LHC).
An excess in the γ-ray emission from the centre of our Galaxy has been discovered in data from the Large Area telescope (LAT), aboard the Fermi satellite [6][7][8][9][10][11][12][13][14][15][16].The nature of the so-called Fermi GeV excess remains a mystery.Several explanations have been put forward, the most exciting of which is perhaps DM annihilation in the halo of the Milky Way (see for example [13,14,17]).Astrophysical processes have also been suggested: the emission from a population of dim unresolved sources [18][19][20][21], and the inverse Compton emission from a new population of cosmic-ray, either from time-dependent events taking place in the past of the Galaxy [22][23][24] or from the high star formation activity in the inner Galaxy [25].
Very recently, two reanalyses of the gamma-ray emission from the inner Galaxy found strong evidence for the Fermi GeV excess being due to hundreds to thousands of dim or unresolved point sources [26,27], most likely millisecond pulsars [18].Even more recently, it was shown that part or all of the required millisecond pulsar population could originate from the disruption of globular clusters by tidal forces in the inner Galaxy [28].Directly detecting some of the millisecond pulsars in the inner Galaxy by radio observations is the next critical step for fully establishing this scenario [26].
It is however striking that the Fermi GeV excess spectrum and spatial distribution are well fitted by what is expected from DM annihilation.The excess could be the first nongravitational signal of DM particles.It is thus urgent to either corroborate or disprove the DM particle nature of the Fermi GeV excess in the framework of concrete models for physics beyond the SM.Supersymmetry (SUSY) is one of the best-motivated classes of renormalizable extensions of the SM, which can accommodate a stable DM particle together with new degrees of freedom that mediate interactions.The most generic R-parity conserving and phenomenologically viable supersymmetric model is the phenomenological Minimal Supersymmetric Standard Model (pMSSM) [29].In this work we address the following key question: can models of the pMSSM explain the observed properties of the Fermi GeV excess, while retaining consistency with other experimental data?And if so, what are the detection prospects for future direct detection and collider experiments?
Until now, MSSM scenarios [30] could not reproduce the Fermi GeV excess as observed by [14], as it was impossible to obtain in this framework a light neutralino (m χ ∼ 30 -40 GeV), as required to fit the Fermi GeV excess spectrum, which could also account for the cosmological DM as measured by Planck [5].However, [15] demonstrated that higher WIMP masses and annihilation channels different from b-quark pairs can give a good fit to the Fermi excess, owing to the freedom allowed by background model systematics.By fitting the GeV excess data of [15], it has been shown that viable solutions in the MSSM exist [31].In the context of the pMSSM, [32] demonstrated that a re-assessment of the theoretical uncertainties in the DM signal spectra opens up a new phenomenology at the LHC experiments.
Here we present the first systematic study of the pMSSM parameter space through global fitting techniques.This approach exhaustively covers all possible phenomenological signatures, allowing us a complete overview of the viable pMSSM interpretation of the Fermi GeV excess.
The paper is organized as follows: in section 2 we briefly present the well-known framework of the pMSSM and the parameters describing the model.In section 3 we describe the experimental set-up of the global fit and the implementation of the joint likelihood.We present the results of the parameter scan in section 4. We then discuss the implications for future direct and indirect detection experiments in section 5.1 and the prospects for detection at the LHC run II in section 5.2.Finally, we summarize our conclusions in section 6.

The theoretical framework: the pMSSM
We here study the pMSSM [29], in which the number of free parameters can be reduced to 19, given the present lack of experimental evidence for SUSY and no experimental indication that one requires the full freedom of a 22-dimensional pMSSM at present.
In this model, the lightest supersymmetric particle is the lightest neutralino, χ, a combination the neutral electroweak gauginos and higgsinos fields.The neutralino is one of the most well-motivated particle DM candidates since it is neutral, stable over cosmological timescales and can lead naturally to the correct DM relic abundance in the early Universe.In what follow, we assume that the neutralino is the particle DM candidate although we do not impose it to fully account for the DM relic abundance as measured by Planck but we allow for subdominant contributions to the DM content (cf.section 3).
We assume first and second generation mass universality, separately in the lepton and quark sectors (table 1).The trilinear couplings of the sfermions enter in the off-diagonal parts pMSSM parameters and priors Flat priors Log priors 173.2 ± 0.87 [33] (Gaussian prior) ρ 0 [GeV/cm 3 ] 0.4 ± 0.1 [34] (Gaussian prior) Table 1.pMSSM parameters and top mass value used in this paper and the prior range for the two prior choices adopted in our scans."Flat priors" are uniform on the parameter itself (within the ranges indicated), while "Log priors" are uniform in the log of the parameter (within the ranges indicated).
of the sfermion mass matrices.Since these entries are proportional to the Yukawa couplings of the respective fermions, we approximate the trilinear couplings associated with the first and second generation fermions to be zero, while the parameters A t , A b and A τ represent the third generation trilinear couplings.In our set-up, the Higgs sector is fully described by the ratio of the Higgs vacuum expectation values tan β, the higgsino mass parameter µ and the mass of the pseudoscalar Higgs m A , which are more directly related to the phenomenology of the model.This 19-dimensional realization of the pMSSM encapsulates all phenomenologically relevant features of the full model that are of interest for DM and collider experiments.The model parameters are displayed in table 1, along with their prior ranges.All of the input parameters are defined at the SUSY scale √ m t1 m t2 .

The experimental setup
We implement experimental constraints with a joint likelihood function, whose logarithm takes the following form: where L GCE represents the Fermi GeV excess, L EW electroweak precision observables, L B(D) B and D physics constraints, L Ωχh 2 measurements of the cosmological DM relic density, L LUX (L IC ) direct (indirect) DM detection constraints and L Higgs (L SUSY ) Higgs (sparticles) searches at colliders.We discuss each component in turn: L GCE : For the Fermi GeV excess likelihood we follow the treatment in [15], to account astrophysical uncertainties.In addition, we include uncorrelated 10% uncertainties as DM modeling systematics, following ref.[32].We marginalize the likelihood function over the uncertainties in the Galactic centre J-value, which is assumed to follow a log-normal distribution with mean log 10 J/(GeV 2 cm −5 sr −1 ) = 23.29 and variance ∆ log 10 J/(GeV 2 cm −5 sr −1 ) = 0.37 [17].When the predicted neutralino relic density, Ω χ , is smaller than the Planck measurement, Ω DM , we follow [35] and adopt the so-called "scaling Ansatz", i.e. we assume that the local ratio of neutralino (ρ χ ) to total DM densities (ρ DM ) is equal to that for their cosmic abundances: L EW : This implements constraints from Z-pole measurements at LEP [36].We include the constraint on the effective electroweak mixing angle for leptons sin 2 θ eff , the total width of the Z-boson Γ Z , the hadronic pole cross-section σ 0 had , as well as the decay width ratios R 0 l , R 0 c and the asymmetry parameters A l , A b , A c and A 0,l F B , A 0,c F B .In addition, we also use the measurement of the mass of the W boson m W from the LEP experiment [36].We apply a Gaussian likelihood for all of these quantities, with mean and standard deviation as given in table II of [37].
The flavor observables related with B and D physics considered are BR( B → X s γ), ).For all we apply a Gaussian likelihood and for most of them we use the measurements following table II of [37].The experimental values assumed for BR(B s → µ + µ − ) and BR(B d → µ + µ − ) are (2.9 ± 0.8) × 10 −9 and (3.6 ± 1.55) × 10 −10 (theoretical uncertainties included [38]) [39] respectively.L Ωχh 2 : We include the Planck Cosmic Microwave Background data constraint on the DM relic abundance as an upper limit, to allow for the possibility that neutralinos are a subdominant DM component.We follow the formalism in the Appendix of [40], using as central value the result from Planck temperature and lensing data Ω χ h 2 = 0.1186±0.0031[41] with a (fixed) theoretical uncertainty, τ = 0.012, to account for the numerical uncertainties entering in the calculation of the relic density.
L LUX : For DM direct detection we include upper limits from the LUX experiment [42], as implemented in the LUXCalc code [43], including both the spin-independent and spindependent cross-sections in the event rate calculation.We adopt hadronic matrix elements determined by lattice QCD [44,45].
L IC : This implements conservative upper limits on the proton spin-dependent cross-section from the IceCube detector in its 79-string configuration [46] (IC-79).Comparable -if slightly weaker -limits have been set for the W W channel by Super-Kamiokande [47] and ANTARES [48].The most stringent constraint is for the case where WIMPs annihilate exclusively to W W pairs. Since the neutrino spectrum generated by Z bosons is similar to the W bosons we apply this constraint whenever the combined branching fraction to W W and ZZ is above 95%.In that case the likelihood is a step function smeared with half a Gaussian (as in Eq. (3.5) of [49]) to account for theoretical and experimental uncertainties that we set to be 50% of the predicted value.
L Higgs : The likelihood for the Higgs searches has two components: the first implements bounds obtained from Higgs searches at LEP, Tevatron and LHC via HiggsBounds [50], which returns whether a model is excluded or not at the 95% CL.The second component constrains the mass and the production times decay rates of the Higgs-like boson discovered by the LHC experiments ATLAS [51] and CMS [52].We use HiggsSignals [53] assuming a theoretical uncertainty in the lightest Higgs mass calculation of 2 GeV.
L SUSY : SUSY searches constraints at LEP and Tevatron follow the likelihood used [49].We have imposed the strict constraints from a large number of searches for supersymmetry at the LHC.The branching ratios of the sparticles have been calculated with SUSYHIT 1.5 [54].We have generated the hadronized event samples with Pythia 8.2 [55] and have employed the NNPDF 2.3 parton distribution functions [56].The generated events are passed on to CheckMATE [57], which is based on the fast detector simulation Delphes 3.10 [58].CheckMATE tests if the model point in question is excluded or not at 95% confidence level by comparing to current experimental searches at the LHC for supersymmetry in the relevant hadronic and leptonic final states with large missing transverse momentum.We assign a 0 log-likelihood L SUSY = 0 if the point passes all constraints, and exclude it if it fails any of them.
We have only included observables we consider robust in order to be conservative.For instance, we have dropped the electroweak precision observables R 0 b and A 0,b F B in the fit because it is unclear whether the large deviations of 2.5σ that are observed with respect to the SM predictions are due to unknown systematic uncertainties or to new physics.The experimental status of the magnetic anomaly of the muon, a µ = 1 2 (g − 2) µ remains unclear in the face of persistent discrepancies in the determination of the hadronic vacuum-polarization diagram using either e + e − or the hadronic τ −decay data.We do not include DM searches in dwarf spheroidal galaxies by the Fermi-LAT [59], as constraints are given only for 100% branching ratio into final states.
Finally, we refer to [37] for details about how the SUSY spectrum and observables are computed.
We use the MultiNest [60] algorithm as implemented in SuperBayeS-v2.0, to perform a global fit of the pMSSM parameter space, including all the data in eq.(3.1), excepting the SUSY searches at the LHC.This is because the LHC searches evaluation is computationally too expensive to be performed on-the-fly.Our scans were run using both log and flat priors to ensure a complete coverage of the parameter space, gathering ∼ 10 6 samples from ∼ 10 8 likelihood evaluations.Samples have been thinned by a factor of 10, focusing our search to regions of the parameter space that were not clearly ruled out by LHC-Run I constraints.This produced 10 5 representative samples to which the LHC SUSY searches have been applied.The ensuing ∼ 10 4 samples that pass LHC run I constraints are displayed in figure 1.

Results
Our global fits identify two distinct viable solutions in the pMSSM parameter space (figure 1): the first exhibits a WIMP mass of ∼ 80 − 100 GeV, with the neutralino annihilating to W W Red crosses indicate the best-fit points in the two islands.The Fermi dwarfs limit for the W + W − channel [59] is plotted for reference only, and it has not been applied.To compare with Fermi dwarfs data, the annihilation crosssection needs to be rescaled by a factor f 2 χ , which would suppress the signal well below current limits.The spin-dependent and spin-independent cross-sections have been multiplied by f χ = Ω χ /Ω DM and g χ = ρ 0 /0.3 GeV/cm 3 to facilitate comparison with current and future limits (LUX [42], Xenon1T and a multi-ton liquid Xe detector with 500 t× yr exposure [61]).In the bottom right panel, we display the IC-79 limit [46] used in our analysis.
with a 95% branching ratio.The second solution has a larger neutralino mass, ∼ 180 − 200 GeV, and 87% tt annihilation final states.The overall best-fit point is in the W W region has −2 ln L ≡ χ 2 = 122.0.This is for a fit with 21 free parameters, and 125 Gaussian data points (we do not include limits as their χ 2 is normalized to 0 whenever the constraint is satisfied), so we adopt 104 degrees of freedom1 .Our best fit thus has a p-value of 0.11 versus a χ 2 = 127.6 and a p-value of 0.06 for the tt solution.
It is important to notice that including theoretical uncertainties in the GC fit is crucial in achieving reasonable p-values.The 10% theoretical uncertainty advocated in [32] is a reasonable reflection of the differences in the predicted spectra between current numerical codes.However, in absence of such an uncertainty, the quality of our global fit would degrade to 0.023 and 0.008 for the two best-fit points, respectively.
The contribution to the overall χ 2 for the two best-fit points from different types of observables are plotted in figure 2. The pulls have been normalized by the number of data points in each group, N , to facilitate a visual comparison.We notice that the χ 2 per data point is distributed fairly evenly across observables.There is a slight preponderance in the χ 2 contribution coming from the GC Fermi fit (with N = 24 bins), which is exacerbated if one neglects the 10% theoretical uncertainty in the DM spectra (dashed bars).The contribution to the pull from the LUX likelihood comes almost exclusively from the SD neutron crosssection limit, as the SI constraint is easily satisfied by our best-fit points.
Let us now analyze in more detail both type of solutions.
In the W W region, model points providing a better fit have a neutralino mostly binolike (∼ 80 -90%) with a similar fraction of both wino and higgsino.Besides we find points in which neutralinos can be dominantly higgsinos with a bino fraction as small as ∼ 10%.Those provide a worse fit though, because a large higgsino composition, basically, implies a large annihilation cross section which drops the relic density below the Planck limit leading to a tension with the Fermi GeV excess due to the scaling Ansatz given in eq.(3.2).
argument we adopt is meant to be representative of what one would get in the simplest scenario.
Analogously to Ref. [32], also in the present work, we find a WW solution with a binohiggsino neutralino composition of about 84 -92 GeV.However, in our work, this solution has a slightly worse χ 2 than in Ref. [32].The reason for this is that we adopt slightly different values of the form-factors for the computation of the SD cross-section.Since such a WW solution is now just below the WW IceCube limits, it is punished a bit in the likelihood.
Typically, since third generation squarks enter into the SI cross section at loop level through LSP-gluon effective interactions [62], their contribution can be comparable to treelevel effective interactions mediated by squarks of the first two families when they are light.In this type of solutions, because current constraints allow sbottoms of a few hundred GeV, their contribution can be sizable and, indeed, cancels out the Higgs exchange contributions.This effect allows to relax the tension with LUX data when Higgs exchange contributions are large.
In terms of the impact of LHC Run-I data, we notice that first and second generation squarks as well as gluinos are decoupled.The stops and the heavier sbottom mass eigenstates are also kinematically inaccessible at the LHC Run-I energies.However, the lighter sbottom eigenstate with a mass around 400 GeV can be produced at the LHC at a considerable rate but the sbottoms evade detection from third generation searches due to complicated cascade decays.The sparticles in the electroweak sector are relatively light.There is a large mass splitting between the SU(2) doublet and singlet sleptons.The SU(2) doublet sleptons are light with masses around 250 GeV and narrowly evade detection in searches for two lepton and large missing transverse momentum final states.The lightest wino-like chargino and the second lightest neutralino escape detection since they are almost mass degenerate with the bino-like neutralino.The production rate of the higgsino eigenstates is too small to yield an observable signal at the LHC run I.
Results about the annihilation cross-section are shown in the top left panel of figure 1.One can see that the points with a better fit exhibit σv ∼ 10 −26 cm 3 /s, consistent with the results found in [15].We also show the Fermi dwarfs limit for the W + W − channel [59], but we emphasize that this limit has not been applied in the fit.In order to compare with the constraint coming from Fermi dwarfs observations, the annihilation cross-section needs to be rescaled by a factor f 2 χ (to account for the possibility of sub-dominant dark matter relic density, which translates according to our Ansatz in a correspondingly reduced local density in the dwarfs), which would suppress the signal well below current and future limits from dwarfs.
Regarding DM direct detection, the top right panel of figure 1 shows the spin-independent (SI) cross-section versus neutralino mass plane.In order to facilitate the visual comparison of our pMSSM models with existing and future limits, we have rescaled the theoretical crosssection by a factor f χ (to account for models where the neutralino does not make up all of the cosmological relic density) and a factor g χ ≡ ρ 0 /0.3 GeV/cm 3 .This accounts for the fact that the local density we have used to predict the number of counts for LUX, ρ 0 , is a nuisance parameter which is generally different from the value assumed by the LUX collaboration in deriving their limit [42], namely 0.3 GeV/cm 3 .
The points that appear above the nominal 95% exclusion limit from LUX cannot be excluded because of a combination of effects: (i) our LUX likelihood is slightly less stringent than what has been published by the LUX collaboration (and depicted in figure 1), (ii) our likelihood function for LUX allows for values of the cross section above the 95% limit to be included (albeit penalized by a smaller likelihood value), and (iii) the global likelihood func-tion can -to an extent -compensate for a poor fit to LUX data by gaining an improvement from other data sets.
Since squarks are heavy, the contribution coming from the exchange of a CP-even Higgs dominates and therefore the SI cross-section scales as ∝ |N 11 N 13/14 | 2 .Because, as noticed above, the higgsino fraction in this region is not negligible, it can be large.Indeed, there are model points above the LUX limit allowed by the scaling Ansatz we apply to the local dark matter density and to some extend due to the fact that we vary the local dark matter density.Another interesting feature is that the SI cross-section spans down to ∼ 10 −14 pb.That is possible because the heavy CP-even Higgs contribution can be sizable and cancellations with the lightest Higgs channel might occur [63].
In the bottom panels of figure 1 we display the spin dependent (SD) cross-section for scattering off neutrons/protons (left/right panels).While in the case of SI interactions the contributions for proton and neutron are comparable, the SD cross-sections may differ significantly.However, we find a tight correlation between the SD cross-sections for scattering off neutrons and those off protons in our model points.The SD cross-section is dominated by the exchange of a Z boson and therefore the SD cross-section is largely determined by the higgsino content of the neutralino and likewise for the SI cross-section it can be sizable in this region as it can be seen in both panels.
In terms of dark matter direct detection experiments, at present, LUX data represents the strongest constraint on the SD-neutron scattering cross-section because the Xenon contains neutron-odd isotopes therefore we overlay the LUX constraint properly rescaled as for the SI case.For the SD-proton scattering cross-section, IC-79 represents the strongest current constraint for the particular case when the neutralinos annihilate to a W + W − final state so we show the IC-79 90% CL upper limit [46] , and we have rescaled the value of the SD cross-section by a factor f χ .This assumes that equilibrium between capture and annihilation is reached in the Sun, which is a good assumption for the bulk of the models shown here.In fact, in the region where the neutralino annihilates mainly to W W and to ZZ to a lesser extend, the IC-79 limits apply here and disfavors a large number of model points.
One can see that there are model points above LUX and IC exclusion lines.In those, the higgsino component of the neutralino is dominant over the gaugino one leading to a large SD cross-section for both neutrons and protons.Those points still provide a reasonable fit to the data due to two effects, first because neutralinos with a large higgsino component yields to a relic density sensible below the Planck bound and therefore the scaling Ansatz applies and second because the local dark matter density is a nuisance parameter in our analysis.Beside there is a small fraction of points above the nominal IC-79 W W limit with a branching ratio to W W and ZZ low enough to evade the stringent IC-79 bound.
In the tt region, points providing a better fit have a neutralino dominantly bino-like (∼ 90%) with a ∼ 10% of higgsino.Those points have the characteristic that the neutralinos annihilate to top quark pairs via an exchange of a right-handed stop which is relatively heavy (∼ 1 TeV).This is possible because on top of the non-helicity suppression the neutralinostop-top coupling component, which is proportional to the top quark Yukawa coupling, is sizeable due to the non-vanishing higgssino fraction of the neutralino.We also found points in which the right-handed stops are light (∼ 300 GeV) being almost bino-like (also found in [32]).However those provide a worse fit to the BR( B → X s γ) data because higgssinostop loops have a sizable positive contribution which leads to values above the experimental constraint.
Annihilation into t t through a t-channel stop exchange requires, for bino-like neutralinos, stops with masses of a few hundred GeV.On the other hand, relatively light stops are not able to lift the tree-level Higgs mass to fit a 125 GeV Higgs.In order to enhance the DM annihilation cross sections to match the Planck measurement, the neutralino coannihilates with sneutrinos which have to be about the same mass as the DM particle (∼ 200 GeV).This induces the splitting in the left/right sleptons spectrum.
These benchmark points are characterized by left handed slepton next-to-lightest supersymmetric particle with masses above 200 GeV.The higgsinos have masses around 260 GeV.Since the resulting mass differences between the bino-like LSP and the sleptons/higgsinos are small, the final states are rather soft and thus the detection is suppressed in events with dilepton or trilepton final states and large missing transverse energy momentum.The production rate of the sbottom is quite suppressed and hence avoids detection.The spectrum of the remaining supersymmetric particles is decoupled.
The phenomenology of the model points in the tt region in terms of DM detection is similar to the W W one.The main difference is in that the IC-79 limits does not apply and therefore larger SD cross-sections are possible.The main difference is that the higgsino composition is not as large as in the WW type of solutions and therefore the SD cross sections for neutron interactions are below the LUX current sensitivity.
Lastly, in figure 3, we show the spectrum of the Fermi GeV excess together with the systematic uncertainties associated with the galactic diffuse emission modeling [15].We compare the data with the spectra of the pMSSM model points giving the best global pvalue in the two regions identified in Figure 1.
It is apparent from Figure 3 that the best-fit DM spectra are systematically offset from the mean values of the Galactic center excess spectrum (gray dots and boxes), by about 1 to 2 sigma, and do not provide a good fit to the data at first sight.However, since the systematic astrophysical uncertainties, indicated by the gray boxes (±1σ) are correlated, this still provides an acceptable fit to the data.To illustrate this point, we show with black dots and error bars the excess spectrum where we moved all data points systematically down, according to the freedom allowed by the covariance (the error bars show now statistical errors only).Together with the 10% uncorrelated systematic modeling uncertainty that we adopted for the DM signal, this provides a reasonable fit to the data, with p-values, whereas without the DM signal modeling uncertainties, the p-values would be prohibitively small (see Figure 2 above).

Implications for direct and indirect dark matter searches
Generally, spin-dependent and spin-independent scattering cross sections are driven by the higgsino content of the neutralino.Therefore, as explained above, the sizable higgsino fraction of points in the two best-fit regions imply large SI cross-sections that makes the direct detection prospects promising, although the SI cross-section range spans down to ∼ 10 −14 pb due to cancellations with the lightest Higgs.In figure 1, top right panel, we also display the projected sensitivity limit (defined as the 90% CL exclusion limit) for the Xenon1T experiment and an hypothetical liquid Xe detector with 500 t× yr exposure [61].The latter  experiment essentially saturates the ultimate detection floor set by coherent neutrino scattering [64].Xenon1T data will be crucial in discovering or firmly ruling out models belonging to the tt island and will probe a significant fraction of the parameter space preferred by the first type of solutions.In figure 1, bottom left panel, we overlay the projected 90% exclusion limit for Xenon1T [61].We find that Xenon1T will be able to prove the entirety of the SD neutron scattering cross-section parameter space favoured by our models. 2A multi-ton experiment with 500 t× yr exposure would reach sensitivities a couple of order of magnitudes smaller than the smallest SD neutron cross-sections found in our scan.
The bottom right panel of figure 1 shows that our best-fit points easily evade the constraint set by IC-79 on the SD proton cross-section in the implementation we adopted in this paper.However, an event-level implementation of the likelihood (including the events' energies [66]) would increase the constraining power of the IC-79 limit, to the point that some of the surviving models could be probed [67].
Finally, as for indirect detection, the preferred parameter space is mostly out of reach even for the future 10-yr Fermi analysis of dwarf spheroidal galaxies [68], where an further improvement of current sensitivities by a factor 2 -3 can be expected.Note that points with a large annihilation cross section as shown in the top right panel of figure 1 usually correspond to suppressed relic densities, making these points hard to detect due to the f 2 χ factor suppression of their signal.

Prospects for detection at the LHC run II
The models yielding the best p-values show several interesting properties relevant for LHC searches.In the following we briefly discuss the discovery prospects depending on the produced SUSY particles.Sparticle mass distributions and MSSM parameters relevant for LHC searches are shown in figure 4.
Light squarks: In the top left panel of Fig. 4 we plot the lightest first/second generation squark mass vs the neutralino mass.The first and second generation squarks have masses > 1400 GeV, i.e. above the usual run I constraints.About 70% of the models have squark masses below 2000 GeV3 .The upcoming run II searches for light squarks will exclude many scenarios.In almost all models the left-and right-handed squarks have different branching ratios and tend to decay to the heavy neutralino and chargino states.Cascade decays including W, Z and Higgs bosons are common.
Stop: In the top right panel of Fig. 4 we plot the stop mass vs the neutralino mass.Some models have light stops with masses down to 200 -300 GeV decaying to chargino and a b-jet.
The neutralino has a mass of around 95 GeV.These models are not excluded by current LHC searches [69].Another interesting region also found in [32] has a stop mass of around 200 -220 GeV and a mass of the lightest neutralino around 180 GeV.A slight excess in the ATLAS data prevents exclusion with run-1 data in this region [69].These solutions will likely be tested with early run II analyses.Other solutions yield much heavier stop masses decaying predominantly into the heavier neutralino and chargino states.Dedicated searches for such decays are important.
Sbottom: In the central left panel of Fig. 4 we plot the sbottom mass vs the neutralino mass.
Several model points have a sbottom mass as low as 400 GeV.The points are not excluded in our procedure due to multi-step cascade decays involving heavy neutralinos.Typically, the lighter sbottom state typically has masses around 400 GeV while the lightest neutralino is mostly bino-like with a mass around 90 GeV.The second lightest winolike neutralino lies around 107 GeV.The two heavier neutralinos are higgsino dominated states and have masses around 250 GeV.
If the sbottoms predominantly decayed into a bottom quark and the lightest neutralino, these benchmark points would clearly be excluded.However, the bottom squark decays into all four neutralino as well as the lighter chargino mass eigenstates with comparable rates.The decay modes of the neutralino and chargino mass eigenstates are relatively complex and hence the limits from simplified sbottom searches do not apply.For instance, the second lightest neutralino eigenstate has large hadronic three body decay modes into the lightest neutralino via off shell Z bosons.Moreover, the second lightest neutralino can radiatively decay into a photon and the lightest neutralino.The corresponding lighter chargino eigentstate has relatively large leptonic three body decay modes via off shell W bosons.Finally, the two heaviest higgsino dominated neutralino mass eigenstates mainly decay into electroweak gauge bosons and the lighter electroweakino states.As a consequence, many events have leptons in the final state which are vetoed in the searches for direct sbottom production.In addition, the higher final state multiplicity tends to soften the net missing transverse momentum distribution compared to scenarios with direct sbottom decays into a bottom quark and the LSP.As a result fewer events pass the selection cuts of the relevant sbottom searches.We explicitly tested those light sbottom scenarios against experimental searches at the LHC with the computer tool CheckMATE and confirmed that those points were allowed.
As seen also for the light squarks about 70% of the models would be excluded with a sbottom limit of ≈ 1000 GeV.
Gluinos: In the central right panel of Fig. 4 we plot the sbottom mass vs the neutralino mass.Gluinos have masses > 1600 GeV.An upcoming early run II exclusion on gluinos with a mass up to 2500 GeV would exclude about 15% of these models.Run II searches for first and second generation squarks will thus likely be more constraining.
Sleptons: In the bottom left panel of Fig. 4 we plot the lightest smuon mass vs the neutralino mass.The lightest smuons found in the best fit models have masses < 400 GeV in about 60% of the best fit points.This makes searches for smuons in run II very sensitive to these solutions.
Chargino/neutralino: In the bottom right panel of Fig. 4 we plot the chargino mass vs the neutralino mass.Due to the GeV excess likelihood several neutralinos and charginos are typically light, the lightest having a mass fixed to 80 − 100 GeV for the W W solutions and 180 − 200 GeV for the tt solutions.The higgsino component in the tt and W W solutions typically involves that µ is only slightly larger than these mass scales, leading to 2 more neutralinos and the light chargino at masses around 100 GeV (W W ) or 200 -300 GeV (tt).These states are often mass compressed with the lightest neutralino which makes the solutions evade LHC chargino/neutralino searches so far.Dedicated chargino/neutralino searches will have sensitivity to most models, e.g. by a mono-jet and soft lepton search as proposed also in ref. [32].The Wino mass scale is quite unconstrained and lies between 100 GeV and 1.5 TeV.The Wino will decay to lighter states yielding final states with Z, W or Higgs bosons.
Heavy Higgs: About 50% of the best models have m A < 800 GeV making searches for heavy Higgs bosons very sensitive.Several chargino and neutralino states are light and have a large coupling to A/H/H ± .Consequently heavy Higgs decays to charginos and neutralinos can have huge branching ratios up to 30% competing with top and bottom decays.Dedicated searches for heavy Higgs bosons decaying into final states with W/Z/h with missing transverse momentum would help to constrain these scenarios.

Conclusions
In this paper we addressed the issue of finding model points in the pMSSM that might explain simultaneously the large set of independent data we gathered from astrophysics, cosmology and high-energy particle physics.We showed that no tension exists between currently available particle physics constraints and the interpretation of the Fermi GeV excess in terms of dark matter annihilation in the framework of the pMSSM.Furthermore, we found evidence for two regions that are able to explain the gamma-ray data, while being consistent with other various experimental constraints: (a) a first region where the neutralino is mostly bino-like and the dominant annihilation channel today is 95% into W bosons pairs and (b) a second region where the the annihilation into top-quark pairs dominates and the neutralino is again mainly bino-like.We showed that these models are very appealing since they will be soon in the reach of the next generation of direct detection experiments -Xenon1T will probe the entirety of the best-fit regions thanks to its sensitivity to spin-dependent neutron crosssection -and of the LHC run II, in particular through searches for charginos and neutralinos, squarks and light smuons.

Figure 1 .
Figure1.2D map of the p-values of our fit, showing the annihilation (top left), the spin-independent (top right) and the spin-dependent (bottom panel, left for neutrons scattering and right for protons) cross-sections vs neutralino mass.The color-bar represents the p-value from the global fit.The yellow overlay highlights points that are within 2σ of the Planck relic density.Red crosses indicate the best-fit points in the two islands.The Fermi dwarfs limit for the W + W − channel[59] is plotted for reference only, and it has not been applied.To compare with Fermi dwarfs data, the annihilation crosssection needs to be rescaled by a factor f 2 χ , which would suppress the signal well below current limits.The spin-dependent and spin-independent cross-sections have been multiplied by f χ = Ω χ /Ω DM and g χ = ρ 0 /0.3 GeV/cm 3 to facilitate comparison with current and future limits (LUX[42], Xenon1T and a multi-ton liquid Xe detector with 500 t× yr exposure[61]).In the bottom right panel, we display the IC-79 limit[46] used in our analysis.

Figure 2 .
Figure 2. Contribution to the overall χ 2 for the two best-fit points, grouped by type of observable (see section 3 for details).The pulls have been normalized by the number of data points in each group, N , to facilitate a visual comparison.The dashed bars show the Galactic Centre Fermi likelihood contribution when the 10% theoretical uncertainty is neglected, which would degrade the p-values to 0.023 (for the m χ = 88.3GeV solution) and 0.008 (for the m χ = 188 GeV model).

Figure 3 .
Figure 3. Spectral energy distribution of the Fermi-GeV excess data.Grey dots show Fermi GeV excess mean values w.r.to background model variations, together with associated systematic uncertainties (grey boxes).Black dots represent the excess for a variation of the Galactic diffuse emission contribution (within its systematic uncertainty).The solid lines show the prediction for the pMSSM models that give the best fits in the two regions of figure1and have dominant annihilation channel into W bosons or top quarks (magenta lines in the left and right panel, respectively).The green lines indicate the adopted 10% systematic uncertainties in the spectra.Furthermore, A denotes the required boost-factor with respect to a generalized Navarro-Frenk-White DM density profile with inner slope γ=1.26.The p denotes the p-value of the global fit (including all data).

Figure 4 .
Figure 4. 2D map of the p-values of our fit, showing the lightest squark mass of the first and second generation (top left), the lightest stop (top right), the lightest sbottom (middle left), the gluino (middle right), the lightest smuon (bottom left) and the lightest chargino (bottom right) vs the neutralino mass.The yellow overlay shows points within 2σ of the Planck relic density value.