MSSM A-funnel and the Galactic Center Excess: Prospects for the LHC and Direct Detection Experiments

The pseudoscalar resonance or"A-funnel"in the Minimal Supersymmetric Standard Model~(MSSM) is a widely studied framework for explaining dark matter that can yield interesting indirect detection and collider signals. The well-known Galactic Center excess (GCE) at GeV energies in the gamma ray spectrum, consistent with annihilation of a $\lesssim 40$ GeV dark matter particle, has more recently been shown to be compatible with significantly heavier masses following reanalysis of the background. In this paper, we explore the LHC and direct detection implications of interpreting the GCE in this extended mass window within the MSSM A-funnel framework. We find that compatibility with relic density, signal strength, collider constraints, and Higgs data can be simultaneously achieved with appropriate parameter choices. The compatible regions give very sharp predictions of 200-600 GeV CP-odd/even Higgs bosons at low tan$\beta$ at the LHC and spin-independent cross sections $\approx 10^{-11}$ pb at direct detection experiments. Regardless of consistency with the GCE, this study serves as a useful template of the strong correlations between indirect, direct, and LHC signatures of the MSSM A-funnel region.


I. INTRODUCTION AND MOTIVATION
The Galactic Center (GC) of the Milky Way galaxy is the densest dark matter region in our vicinity and has long been earmarked as the most promising target for searches of dark matter (DM) signals. Intriguingly, recent years have seen a persistent and statistically significant excess in the gamma ray spectrum peaking at 2−5 GeV originating from the GC, above what is predicted from known sources and conventional astrophysics [1][2][3][4][5][6][7][8][9][10][11][12]. The signal was initially reported to be compatible with ∼ 40 (10) GeV dark matter annihilating into bb (τ τ ), with an annihilation cross section σv ∼ O(10 −26 ) cm 3 /s. Since this is approximately the annihilation cross section expected of a thermal relic, a dark matter interpretation of this excess presents itself as a very tantalizing possibility. This prospect has been explored by many authors in various contexts (see, for instance Refs. [9,10,13,14] and references therein), including the Minimal Supersymmetric Standard Model (MSSM) [15][16][17][18]. More recently, it has been shown that this excess might be attributable to unresolved point sources [19][20][21], although a conclusive verdict has not been reached.
Recently, the Fermi-LAT Collaboration has presented an analysis of the region around the GC with four different variants of foreground/background models, finding, for every variant, significant improvements in the agreement with data when an additional component centered at the GC with a peaked profile (NFW, NFW-contracted), i.e. a dark matter-like spectrum, was included in the fits [22,23] (see also Ref. [12] for an attempt at accounting for systematic uncertainties in the background). From a dark matter perspective, a recent study [15] found these additional components for the four choices of background models to be compatible with several annihilation channels (W W, ZZ, hh, tt) and significantly higher DM masses (165 GeV for bb, 310 GeV for tt) than previously thought possible. Similar conclusions were also reached in Refs. [17] and [18], which reported that a higher mass (175 − 200 GeV) dark matter annihilating into tt could give reasonable fits to the signal. This relaxation of the allowed range of dark matter masses compatible with the GC excess (GCE) has particularly interesting implications for MSSM dark matter, as it opens up the possibility of explaining the signal with the well-known pseudoscalar resonance or "A-funnel" mechanism, where the dark matter relic density is set by resonant s-channel annihilation through the pseudoscalar A, with m A ≈ 2m χ ( χ represents the lightest neutralino, which is the dark matter candidate). The pseudoscalar resonance has been studied in con-nection with the GCE outside the MSSM in Refs. [24][25][26]; however, realizing the mechanism in the MSSM is of particular interest given that the MSSM remains one of the most familiar and widely studied Beyond the Standard Model (BSM) theories. Previous fits to the GCE with m χ < ∼ 50 GeV did not allow for this possibility in the MSSM due to constraints on m A from direct LHC searches [27,28] (although this constraint can be circumvented in the the Next-to-Minimal Supersymmetric Model (NMSSM), allowing for an NMSSM explanation of the GCE [25,29]). This incompatibility is lifted if, as discussed in Ref. [15], m χ < ∼ 165 (310) GeV annihilates into bb (tt), allowing for m A large enough to evade collider constraints.
The aim of this paper is to explore whether, given this wider range of allowed masses, the MSSM pseudoscalar resonance can give reasonable fits to the GCE, consistent with stringent constraints from relic density, indirect/direct detection, collider search limits, and Higgs data. Since the mechanism requires a light (∼ 200 − 500 GeV) pseudoscalar, the SM-like nature of the 125 GeV Higgs boson is particularly constraining as the heavier CPeven Higgs is at the same mass as the pseudoscalar and can mix with the 125 GeV Higgs, resulting in deviations from SM-like properties inconsistent with measurements. For such light, non-decoupled heavier Higgs bosons, the Higgs sector needs to be "aligned" [30][31][32][33][34][35] to maintain SM-like properties for the 125 GeV mass eigenstate. As we will show in this paper, this can indeed be achieved while simultaneously satisfying all other DM requirements.
A successful realization of neutralino dark matter along with the GCE through the pseudoscalar resonance requires very precise choices of parameters in order to simultaneously achieve resonant annihilation, the Higgs mass, and alignment in the Higgs sector (this is also the reason why extensive scans in the MSSM parameter space [15][16][17][18] fail to uncover it as a viable explanation of the GCE). It is nevertheless worthwhile to pursue this direction for several reasons. First, the A-funnel is one of several "traditional" mechanisms in the MSSM that have been widely studied for a long time, and its compatibility with a possible DM signal is therefore of considerable interest. Second, while most scenarios put forward to explain the GCE could potentially be constrained by stringent spin-independent direct detection limits (indeed, avoiding these limits itself involves some nontrivial fine-tuning of parameters in supersymmetric models [36][37][38]), the A-funnel naturally gives small direct detection cross sections and is automatically safe from these bounds. Most importantly, the framework is eminently predictive, giving very specific predictions for heavy Higgs bosons that will be probed at the 13 TeV LHC and future colliders, as well as direct detection cross sections that may be probed by the next generation of experiments. Independent of these considerations, and independent of the applicability to the GCE, this study serves as a valuable template of the conditions necessary for the existence of a light pseudoscalar in the MSSM together with indirect detection signals of dark matter via the A-funnel.
The outline of the paper is as follows. Section II introduces the parameter space relevant for the study and discusses dark matter aspects such as the annihilation cross section and relic density. Section III is devoted to a discussion of various constraints from direct detection, indirect detection, collider constraints, Higgs data, and vacuum metastability. Section IV presents the details of our scans and the best fit regions to the GCE. Predictions for the 13 TeV LHC and future direct detection searches are presented in Section V. We summarize our results in Section VI. The Appendices contains additional details on the MSSM parameters and fits to the GCE.

II. THE MSSM PSEUDOSCALAR RESONANCE: DARK MATTER ASPECTS
In R-parity conserving supersymmetric models, the lightest supersymmetric particle (LSP) is stable. If it is also neutral, it can be a dark matter candidate. In the MSSM, the LSP is often assumed to be the lightest of the neutralinos, the neutral superpartners of the gauge bosons and Higgs bosons (Bino, Wino and Higgsinos respectively). The Wino and the Higgsinos tend to annihilate too efficiently to explain the observed dark matter abundance. However, the Bino can yield the correct relic density via various mechanisms, including resonant annihilation via the pseudoscalar, and has long been regarded as the favored dark matter candidate.
We perform our study in the phenomenological MSSM (pMSSM) [39], which is defined in terms of 19 parameters, which are taken to be independent at the weak scale. Of these, our analysis will be entirely determined by the following seven parameters: • M 1 , the Bino mass parameter. The dark matter is mostly Bino, so this is also approximately the mass of the dark matter candidate m χ ≈ M 1 .
• µ parameter. This is the Higgsino mass, and controls the Higgsino fraction in the dark matter particle χ. As we will see later, the relic density, signal strength, and direct detection cross section all depend sensitively on this fraction.
• tan β, the ratio of the up-and down-type Higgs vacuum expectation values (vevs).
• m A , the heavy Higgs mass. This is the mass of the pseudoscalar that mediates the resonance (hence m A ≈ 2m χ ) as well as the mass of the heavier scalar, which feeds into Higgs phenomenology and expected direct detection cross-sections.
• m Q 3 , m u 3 , the left and right handed stop masses, which contribute significantly to the mass of the observed 125 GeV Higgs boson. In this paper we take the stop mass scale • A t , stop trilinear coupling. This determines the mixing in the stop sector and is again a relevant parameter for the mass of the observed Higgs boson.
All other masses, such as the other gaugino (wino and gluino) and sfermion masses, are assumed to be heavy and decoupled from the analysis.

A. Dark Matter Composition
The lightest neutralino in the MSSM is a combination of the Bino, Wino, and neutral Higgsinos: As mentioned above, we are mainly interested in the region of parameters where the lightest neutralino is predominantly a Bino, hence N 11 ∼ 1, N 12 = 0, and N 13 , N 14 1. In this regime, the Bino mass parameter M 1 and the neutralino components are approximately [25] Here, s θ , c θ denote sin θ, cos θ respectively and m χ is the dark matter mass.

B. Relic Density and Signal Strength
Both the relic density and the present day annihilation cross section are driven by the process χχ → ff with the pseudoscalar A in the s-channel (we are interested in the case where the fermion f is either b or t for compatibility with the GCE). When the process occurs close to resonance, it is well-known that the annihilation cross-section in the early universe (which sets the relic density at the time of freeze-out) is substantially different from that at present times (which sets the signal strength fitting the GCE) due to thermal broadening of the resonance during the former stage [40]. Thus, with appropriate parameter choices, one can scale the relic density and the present annihilation cross section independent of each other, thereby achieving better agreement with both measurements; this degree of freedom is not afforded in non-resonant scenarios, where these two quantities are strictly related to each other.
To understand this interplay, consider a simplified model describing a Majorana DM particle χ coupled to a pseudoscalar A through the interaction Lagrangian The entire parameter space of the model is then determined by m A , m χ , y aχχ and y af f . A crucial parameter in our analysis is the degeneracy parameter which characterizes the proximity to the resonant regime. We are interested in scenarios where δ ≈ 0.
The resonant annihilation cross-section at a given temperature T is [40] σv where x = m χ /T and Γ A is the decay width of A, This gives the relic abundance where x f is the value of x at freeze-out. This expression can be rewritten in a more illuminating form as [25] Ωh 2 ∼ 0.12 Likewise, the DM annihilation cross-section today is σv| v=0 3 2π Assuming that m A ∼ 2m χ so that the second term dominates in the denominator, one obtains (for 2m χ < m a ) [25] σv| v=0 ∼ 2 × 10 −26 cm 3 4m 2 Comparing Eq. 8 and Eq. 10, it is clear that the relic density and the current annihilation cross-section can be independently scaled with judicious choices of y af f and y aχχ δ/6. In terms of the fundamental MSSM parameters, these couplings are given by: where v = 174 GeV and g 1 is the SM U (1) Y gauge coupling. Note from the above that a non-vanishing y aχχ coupling requires a non-vanishing Higgsino component in χ. From the expressions for N 11 , N 13 , N 14 listed previously, we thus see that, for given values of m A and tan β, the desired relic density and an annihilation cross-section consistent with the GCE can be obtained simultaneously by appropriately choosing µ and δ (equivalently, m χ ).

III. CONSTRAINTS
As mentioned in Sec. I, the relevant A-funnel parameter space is constrained from several directions. Higgs phenomenology in our set-up is very directly linked to the GCE, hence LHC direct searches as well as the properties of the observed 125 GeV Higgs put stringent constraints on this scenario. Consistency with all collider observables can then create tension with constraints from requiring the stability of the electroweak vacuum. In addition, since the CP-even heavy Higgs H is expected to be approximately degenerate in mass with A, contributions to the spin-independent direct detection cross-section from H-exchange might be relevant. Finally, there are also several current and future indirect detection experiments that can probe the process of interest in this paper. In this section we detail the current status and future prospects in all of these different directions.

A. Collider and Higgs Sector Constraints
In the absence of CP-violation (which we assume in this paper), the physical spectrum GeV. One can define a "Higgs-basis", where a single field acquires all the vev: 1 Light m A /m H and heavily mixed stops (as usually needed for a 125 GeV Higgs in the MSSM) can also give large contributions to various flavor observables, for example B s → µ + µ − and B → X s γ. However, in this work we will mainly be interested in moderate to small value of tan β, hence there is no large enhancement of these effects. Moreover, the size of these contributions are heavily dependent on the signs of various contributions (see e.g. Ref. [42]), and consistency with all measured values could be obtained by tuning such cancellations.
where s β ≡ sin β, c β ≡ cos β, H SM = v, and H N SM = 0. The couplings of these states to the SM fields are: where V V, uu, dd refer to all vector, up-type and down-type states respectively, and g SM refers to the SM value of these couplings. Note that there is no coupling between the where , and α is the angle of rotation from the (H u , H d ) basis to the mass eigenstates. We want to identify the lightest CP-even mass eigenstate, h, with the recently observed 125 GeV scalar; given that all measurements suggest that its properties are SM-like, we also want to identify it as the SM-like field in the Higgs basis. That is, we require This requirement of vanishing mixing between the H N SM state and the 125 GeV Higgs, corresponding to κ h N SM ≈ 0, can be rewritten in terms of the fundamental parameters as [30,31] where M S is the geometric mean of the stop masses and Note that when the second Higgs becomes heavy (m H >> m h ), this relation is automatically satisfied; this is the familiar decoupling effect. Otherwise, one requires alignment without decoupling [30,31], brought about by an accidental cancellation in the fundamental parameters of the theory so as to satisfy Eq. 18. For small t β and M S ∼ O(1) TeV, large values of A t /M S are required to obtained an experimentally consistent Higgs mass whereas large values of (µA t )/M 2 S lead to close to alignment conditions [30,31]. The CMS and ATLAS collaborations present both the precision measurements of the 125 GeV Higgs and the searches for H → W W/ZZ as ratios to the expectations from a SM Higgs of the same mass. The predicted rate at the LHC for the decay of the mass eigenstate i = {h, H} into some final state XX as a ratio to the SM value is given by where  [44]. In our analysis we will take a conservative approach of assuming that observational consistency is obtained (that is, the Higgs sector is sufficiently aligned) for R h W W between 0.7 − 1.3. This range will narrow with additional data, and measurements at the level of 10% are expected at the high luminosity LHC [45,46].

B. Vacuum Metastability
Another important constraint on these parameters comes from vacuum metastability.
Large values of the soft stop trilinear coupling A t , required for the Higgs mass and alignment (discussion above), can result in the appearance of charge-and color-breaking minima in the scalar potential of the MSSM. The condition for either these minima to be energetically unfavorable or the tunneling to these minima to have lifetimes longer than the age of the Universe leads to the approximate bound [47] where Hu + µ 2 , and r = m 2 u 3 /m 2 Q 3 . In our analysis we assume Minimization conditions of the Higgs potential give , hence the condition for vacuum metastability can be written as It is worth keeping in mind that this is only an approximate bound and depends on several assumptions (see Ref. [47] for details). However, consistency with the above provides a rough guide for the feasibility of the parameter region under investigation.

C. Direct Detection
Direct detection possibilities focusing on the A-funnel in the MSSM have been studied in Refs. [48][49][50]. The pseudoscalar A does not mediate spin-independent WIMP-nucleon scattering. Instead this cross section σ SI comes from light and heavy CP-even Higgs boson exchanges in the t-channel, facilitated by the Bino-Higgsino mixture of the LSP necessary to obtain the correct relic density. There are also contributions from tree level squark exchange in the s-channel and from gluon loops [51,52], but these are negligible when the sfermions are heavy. The cross section then depends only on M 1 , m A , tan β and µ.
For given values of m A and tan β, requiring the correct relic density and GCE leaves no free parameters, thereby fixing the direct detection cross section. This cross section in our region of interest can be written as approximately [25] σ SI where F u ∼ 0.15, F d ∼ 0.13 (the up and down type quark content respectively of the nucleon, proton or neutron), t β = tan β, m N is the mass of the nucleon, and m r = m N mχ m N +mχ is the reduced mass. For the correct dark matter relic density obtained via the A-funnel, this cross section is generally around 10 −11 pb [48][49][50]53], well below existing bounds from XENON100 [54] and LUX [55], which currently rule out σ SI > ∼ 5 × 10 −10 pb. Note that while the annihilation processes that determine the relic density as well as indirect detection signals are s-channel and therefore enhanced by the resonance, the direct detection cross-section is mediated by t-channel processes and does not receive this enhancement. Such small direct detection cross sections are therefore a generic feature of this region of parameter space.
Crucially, this cross section still lies above the neutrino background and is therefore within reach of future detectors, although detection will still be challenging.
As is well-known, an exception to this generic feature can occur for negative values of the µ parameter due to destructive interference between the light and heavy Higgs exchange contributions, giving cross sections several orders of magnitude below the neutrino background cross section [48,50]. Such blind spots can in general occur at any dark matter mass, but their appearance in the A-funnel framework is more strongly constrained as we also need m H ∼ m A ∼ 2 m χ . Approximating the up-and down-type quark content in the nucleus as roughly equal, this cancellation condition in the A-funnel region can be formulated as approximately [48]

D. Indirect Detection
Currently the strongest bounds on the annihilation cross section are given by the Fermi/LAT analysis of 6 years of data on 15 known dwarf galaxies [56]. For 100 − 300 GeV dark matter, which is our region of interest, this analysis constrains the annihilation cross-section to be less than ∼ a few×10 −26 cm 3 /s. The cross section required to explain the GCE is also in this region over this mass range (see [15]), hence the dwarf constraints are in some tension with a DM interpretation of the GCE. However, the large uncertainties in the dark matter distribution (J-factor) in these dwarf galaxies leave room for compatibility (see Fig. 8 in Ref. [56]). For instance, the 95% C.L. annihilation cross-section exclusion limit for a 100 GeV WIMP annihilating to bb is 2.  [59]. Bounds similar to those from the Fermi dwarf observations are also found by the Planck satellite from CMB measurements [60].
Likewise, since DM of interest in this paper annihilates primarily through hadronic channels (bb and tt), this is expected to generate a significant flux of antiprotons. There already exists some tension between models that explain the GCE and derived constraints from antiproton bounds on dark matter annihilation [61][62][63]. However, calculation of the antiproton flux suffers from significant uncertainties related to the propagation model in the galaxy (see [63][64][65][66] and references therein), and the GCE can be made compatible with the measured antiproton flux for conservative choices of propagation model parameters.
Bounds on the dark matter annihilating cross-section into quarks are also obtained by neutrino experiments like IceCube. The most current results from the IceCube-79 experiment exclude σv ≥ 2 × 10 −22 cm 3 /s into bb at 90% confidence level [67]. This lower limit is ∼ 10 4 larger than the cross-section required for the GCE [15] and thus irrelevant.
Therefore, no indirect detection results robustly rule out a DM interpretation of the GCE at present, although future measurements, particularly from Fermi-LAT observation of dwarfs, AMS-02 antiproton results, and the CMB could have interesting implications.

IV. NUMERICAL RESULTS
Building on the parameter space and constraints described in the previous sections, we present the fits to the GCE excess in this section. We used the following tools for our numerical analysis: the neutralino relic abundance and annihilation crosssection was calculated with Micromegas-4.1.7 [68], the MSSM particle spectra were computed using SuSpect-2.41 [69], and the Higgs phenomenology was obtained with FeynHiggs-2.11.0 [70][71][72][73][74].
For the gamma ray spectrum corresponding to the signal, we follow the approach employed in Ref. [15] and consider two of the four spectra presented in Fig. 13 of Ref.
[23] 2 , 2 The first version of our paper used the spectra presented in Ref. [22], and Ref. [23] is the corresponding which were derived by fitting the excess over various choices of background as exponentially cut off power laws (see Ref. [15,22] for further details). The four spectra are referred to as spectra (a)-(d) in Ref. [15], and just as they do, we pick spectra (b) and (d) for our analysis; spectrum (a) is very similar to what has been studied for light (m χ < ∼ 40 GeV) DM in previous papers and not amenable to the MSSM, whereas spectrum (c) is very similar to spectrum (d) and does not yield any new insight.
Spectrum (b) corresponds to a fit with OB stars as cosmic ray (CR) sources and a tuned index for pion production within the solar circle (see [22,23]); the analysis in Ref. [15] found it to be well fit by 75 − 95 GeV DM annihilating into bb or < ∼ 200 GeV DM annihilating into tt. Annihilation into gauge or Higgs bosons were also found to give good fits, but these are irrelevant for our analysis since they are always subdominant channels in the MSSM pseudoscalar resonance scenario. Note that spectrum (b) is also in agreement with other studies performed in Refs. [17] and [18], which also found that 175−200 GeV DM annihilating into tt could be compatible with the GCE. Likewise, spectrum (d) corresponds to a fit with OB stars as cosmic ray (CR) sources but with only the intensity of pion production tuned (using pulsars instead of OB stars gives a very similar spectrum); Ref. [15] found it to correspond to higher mass DM, with 130 − 165 GeV DM annihilating into bb or 250 − 310 GeV DM annihilating into tt giving good fits.
In this section, we will perform fits to the two spectra (b) and (d) with the idea of gaining intuition about the range of possibilities that the GCE allows for the MSSM pseudoscalar resonance. We note that the continuous region spanning spectra (b) and (d) could also plausibly explain the GCE for some reasonable background, but do not pursue this direction any further.

A. Fit Procedure
The astrophysical information regarding the distribution of dark matter is encoded in the publication that recently appeared; we have chosen the spectra from Ref. [23] that correspond most closely to the spectra we used in the first version.
where ∆Ω is the region of interest (ROI), l.o.s. stands for line of sight, and ρ is the dark matter density.J can. = 2.0 × 10 23 GeV 2 /cm 5 is the canonical value of the J-factor obtained from evaluating the integral with an NFW profile. Following the analyses in Ref. [15], we parametrize the uncertainty in the dark matter density profile with the factor J , which is allowed to vary between [0.14, 4].
The gamma-ray spectrum is computed for the following MSSM parameters: • The pseudoscalar mass is allowed to vary over 200 GeV ≤ m A ≤ 700 GeV. Below 200 GeV, we find that the Higgs sector cannot be sufficiently aligned while remaining consistent with bounds from H/A → τ + τ − from the 8 TeV LHC run. We terminate the scan at 700 GeV since good fits to the GCE (either spectrum (b) or (d)) are not expected for m χ ≥ 310 GeV.
• tan β is scanned over the range 4 ≤ tan β ≤ 10. Below tan β = 4, extremely heavy • For given values of m A and tan β, we next scan over δ (equivalently, m χ as shown in Eq. 4) and µ for points such that the relic density constraint is satisfied: the neutralino makes up all of dark matter (0.08 ≤ Ωh 2 ≤ 0.16); and the annihilation cross section σv is within the 2σ best-fit annihilation crosssection contours from Ref. [15].
We scan over δ ∈ [0, 0.1] in order to stay close to resonance, and over µ ∈ [0.7, 10] TeV in order to obtain a mostly bino DM. We take the branching ratio to W W normalized to the SM value R h W W to be a measure of alignment and select (for each m A , tan β, µ, m χ combination) the combination of M S and A t that gives R h W W closest to 1 while maintaining 122 ≤ m h ≤ 128 GeV.
• All other MSSM input parameters (gaugino/wino masses, trilinear couplings, slepton/squark masses) are set to 5 TeV so that they decouple from this analysis.
The goodness of fit is obtained by performing a χ 2 analysis between the gamma-ray spectrum obtained from Micromegas and the GCE (Fermi spectra (b) and (d)). For a given MSSM point, the χ 2 is calculated as: where the subscript k runs over the 20 energy bins of the Fermi/LAT measurement [22], dN/dE is the gamma-ray spectrum obtained from Micromegas, the subscript obs denotes the spectrum consistent with the Fermi excess (i.e. spectrum (b) or (d)), σ k denotes the statistical uncertainty [15], and J is the value of J ∈ [0.14, 4] that minimizes the χ 2 value. The χ 2 analysis includes statistical errors, but neglects possible systematic errors from modeling backgrounds near the Galactic Center.

B. Fit Results
The fits resulting from the above procedure are presented in Fig. 1 as contours of χ 2 in the m A -tan β plane for Fermi spectrum (b) and (d). The pink crosses in each panel denote the points with the best fit to the corresponding spectrum; the gamma-ray spectra of these best fit points are presented in Fig. 2 along with the MSSM parameters 3 In Fig. 1 we also include, in solid black lines, the 1-σ and 2-σ bounds from A/H → τ + τ − searches at the 8 TeV LHC [41]; points that lie above these curves in the shaded region are inconsistent with these bounds. These τ τ searches, however, lose sensitivity at low tan β, hence light pseudoscalars can mediate DM annihilations capable of explaining the GCE in this region. The dashed blue corresponding to "OB stars index scaled" and "OB stars intensity scaled" spectra from Fig. 13 of Ref. [23] (see Ref. [15,22,23]  This generically requires close to resonance conditions 2m χ ≈ m A for consistency with both the GCE and relic density. We found that the χ 2 value did not change significantly between distinct values of (µ, δ, A t , and M S ) for the same m A , tan β. This is expected, since the fit quality is driven by the shape of the spectrum, which is controlled mainly by tan β via the branching ratios, and the position of the peak, which is controlled by m A (≈ 2m χ ). Although the fit should also depend on the signal strength, which is controlled by µ and δ via the annihilation cross section and relic density, the freedom in choosing J ∈ [0.14, 4], which essentially rescales the signal strength, smears out this dependence. In our region of interest, we find that δ < ∼ 0.04 while M S , A t , and µ all take multi-TeV values; we present contour plots of these parameters in Fig. 8 in Appendix A. The condition for vacuum metastability, Eq. 22, is also found to be satisfied in most parts of the parameter space allowed by the 8 TeV LHC A/H → τ + τ − bounds (see Fig. 9 in Appendix A).
From the left panel of Fig. 1, the best fit regions to Fermi spectrum (b) appear to be separated into two distinct islands. The m A < ∼ 250 GeV region has relatively low χ 2 for all values of tan β. In this region, annihilation into top quark pairs is kinematically forbidden, so the dominant annihilation channels is bb for all values of tan β. Recall that an approximately 100 GeV DM particle annihilating into bb can fit the GCE [15]; this region reflects this behavior. However, we see that this region is incompatible with the 8 TeV LHC A/H → τ + τ − bounds and/or the Higgs data (that is, R h W W < ∼ 0.7 in this region, signaling that the heavier CP-even scalar is so light that alignment does not work well). A second island opens up at 350 GeV < ∼ m A < ∼ 450 GeV, when annihilation into tt becomes kinematically feasible, and tan β < ∼ 6. This is consistent with Ref. [15] finding a ∼ 200 GeV DM annihilating into tt providing a good fit to spectrum (b). Note that the best fit point occurs at the lowest allowed value of tan β(=4) in our scan, where the coupling of A to top quarks is the largest. The fit deteriorates as tan β gets larger, as the branching ratio into bb gets larger due to the tan β enhancement of the Abb coupling. This region is also compatible with Higgs data as R h W W > ∼ 0.7, and safe from the current A/H → τ + τ − bounds. Beyond this island, the fit deteriorates rapidly as m A and/or tan β are increased.
Similar patterns are observed for the fit to spectrum (d). A small region of good fit exists at m A ∼ 300 GeV and low tan β, safe from the A/H → τ + τ − bounds and borderline compatible with Higgs data. Again, DM in this region annihilates dominantly to bb since tt is kinematically forbidden, and this observation is compatible with Ref. [15], where DM with mass 130 − 165 GeV annihilating into bb was found to give good fits to the spectrum.
A second region with better fits is again observed for larger m A once decay into tt opens up. This regions roughly spans 450 GeV < ∼ m A < ∼ 600 GeV and tan β < ∼ 8, and appears to correspond to the 250 − 310 GeV DM annihilating into tt region reported in Ref. [15] as a good fit to spectrum (d). Similarly to spectrum (b), the best fit occurs for small values of tan β : tan β ∼ 4.0. This suggests that a DM candidate that annihilates significantly into tt with BR(χχ → tt) = 0.66 at the best fit point) provides the best fit to spectrum (d). This can be confirmed by comparing the shape of the spectrum in Fig. 2, right panel, which fits the shape of Fermi spectrum (d) quite well. Finally, the fit deteriorates for larger m A and tan β values and we do not expect any good fits beyond the region shown in the plot.

Fit to a Modified Spectrum
So far, we performed fits to spectra (b) and (d) as defined in Ref. [15], corresponding to the "OB stars index scaled" and "OB stars intensity scaled" spectra from Fig. 13 of Ref. [23], which were obtained by modeling the excess with an NFW profile with a single power law with an exponential cutoff. This mimics what is expected of a dark matter source, Red (blue) contour regions denote the best (worst) fits. Black and blue contours are as in Fig. 1.
and serves the purpose of demonstrating how the preferred theory parameter space changes for two different choices of interstellar emission models of the background (matching the philosophy in Ref. [15]). However, Ref. [23] also finds significantly better fits to the excess if more freedom is allowed in the fit -in particular, if the spectrum of the NFW profile is modeled with a power-law that is allowed to vary per energy band over the 1 -100 GeV range; the resulting spectra for various choices of interstellar emission models are presented in Fig. 18 of Ref. [23]. In order to study how the MSSM fit is affected if the latter is used, we performed a similar fit (as described above) to the "pulsars index-scaled" spectrum from Fig. 18 of Ref. [23]; the result is shown in Fig. 3. We find that the overall fit quality worsens due to the tail of the spectrum, but the best fit regions in the MSSM parameter space still closely match those from the fit for spectrum (b) (see Fig. 1 (left)); consequently, the theoretical implications from fitting to spectrum (b) (discussed below) will also apply in this case.

A. LHC Prospects
There are several projections for the 14 TeV LHC provided by the CMS and ATLAS collaborations for heavy Higgs searches in [75,76]. The top row shows that both BR(A → τ τ ) and BR(A → Zh) are a few percent throughout the parameter region of interest, with the former always comparable to or larger (in some cases, by more than an order of magnitude). We can understand this behavior by noting that due to the close to alignment conditions, the AZh coupling is very suppressed.
Hence, despite the tan β enhancement of the gluon fusion production of A, we find that the rates for A → Zh are at least 2 orders of magnitudes smaller than the current exclusion limits [78,79] and therefore unlikely to be probed even at the high luminosity LHC [75,76].
Due to the absence of any other relevant decay modes, the decays to down-type fermions will still be the dominant decay modes and offer the best prospects for discovery of the pseudoscalar.
For the heavier CP-even Higgs H, in addition to the τ + τ − channel, there are nonnegligible branching ratios into W W or hh despite being suppressed due to alignment (recall that, close to alignment, H ≈ H N SM ). These branching ratios are largest at low tan β below the top mass threshold, whereas Br(H → τ + τ − ) is larger at higher tan β. Note again that in the low tan β region, the main production of H is via gluon fusion, which is enhanced     [80], hence dedicated searches at the LHC could probe the GCE best-fit regions, particularly for m A < ∼ 350 GeV, where R H W W can be within a factor of 10 of the current exclusion limit [75,76].
For H/A heavier than about 350 GeV and low values of tan β ( < ∼ 7), both the CP-odd and even Higgs bosons preferentially decay to top quark pairs. However, due to the large SM tt background, this is a very challenging signature for the LHC [81,82]; nevertheless, stronger sensitivity is expected at a 100 TeV collider [82]. The standard τ + τ − searches can probe regions with larger values of tan β.
It should be kept in mind that, in addition to these searches for heavier Higgs bosons, the good fit regions at low m A < ∼ 350 GeV also predict deviations in R h W W (see Eq. 20 for definition) at the 10% level or more, hence such deviations from SM-like properties of the 125 GeV Higgs could be a stark signal of this scenario. All of the above search modes as well as the precision measurements of the 125 GeV Higgs are expected to improve substantially in sensitivity with the higher luminosity and energy of the 13 TeV LHC [45,46].

B. Direct Detection
Our predictions for spin-independent direct detection experiments are plotted in Fig. 7, which shows DM masses and spin-independent DM-nucleon (proton) direct detection cross sections compatible with the GCE (Fermi spectrum (b) in blue, spectrum (d) in red). We only show points with χ 2 ≤ 50 that are compatible with both the 2σ A/H → τ + τ − 8 TeV LHC constraints and 0.7 ≤ R h W W ≤ 1.3. As discussed in Section III C, we see that DM via the pseudoscalar resonance corresponds to generic cross sections of O(10 −11 )pb, and these are comfortably safe from the existing Xenon100 [54] and LUX [55] bounds. A major fraction of the predicted parameter space can be probed with the next generation of direct detection experiments such as Xenon1T and LZ [83]. We note that almost all points predicted from our fit lie above the neutrino floor and therefore a signal can in principle be detected. The green cross and star correspond to the best fit points from Fig. 1 for spectrum (b) and (d) respectively.

VI. SUMMARY
To conclude, we summarize the main findings of this paper: • Recent reanalysis of GC background has found that the GCE could be consistent with annihilation of DM with much higher masses [15,17,18,22]. This allows the GCE to be explained by the MSSM pseudoscalar resonance or "A-funnel". We fit to two different dark matter spectra, Fermi spectrum (b) and (d) from [15,22], and find that reasonable fits can be obtained while maintaining consistency with stringent constraints from collider searches, Higgs data, and direct and indirect detection.
• For spectrum (b), the best fit region corresponds to 350 GeV < ∼ m A < ∼ 450 GeV and tanβ < ∼ 6. This region can be probed with searches for H → W W and tt resonance searches. m A < ∼ 250 GeV also gives reasonable fits but is incompatible with Higgs data.
• For spectrum (d), there are two regions with reasonable fits to the GCE: 450 GeV < ∼ m A < ∼ 600 GeV at tanβ < ∼ 8, and m A ∼ 300 and tanβ < ∼ 5.5. The former region can yield signals at the LHC in the A/H → τ τ or tt resonance searches at the LHC. The latter region can also be probed with the same channels, and should also lead to measurements of deviations of the 125 GeV Higgs couplings from SM-like values.
• The best fit regions for both spectra (b) and (d) predict spin-independent direct detection cross sections of O(10 −11 )pb for a 110 GeV < ∼ m χ < ∼ 350 GeV neutralino. The entire region lies above the neutrino background, and the majority of the region is within reach of Xenon1T and LZ (see Fig. 7).
This exercise therefore leads to very sharp predictions for the next round of the LHC and direct detection experiments. Although the best fits obtained in this paper are noticeably worse than the best fit dark matter scenarios discussed elsewhere in literature, this highly predictive framework, coupled with the wide popularity of the MSSM, makes these results noteworthy. Even if the GCE turns out to be incompatible with the MSSM pseudoscalar resonance and is ultimately explained by some other (dark matter or astrophysical) phe-nomenon, this study still serves as a valuable template for the interplay between existing collider and Higgs constraints and the indirect, direct, and collider signatures of the A-funnel region with a light pseudoscalar in the MSSM.
Appendix A: Parameters and Vacuum Metastability Fig. 8 presents contour plots of the scanned parameters in the m A -tan β plane. The approximate check for vacuum metastability from Eq. 22 is shown in Fig. 9. It is seen that the desired condition is satisfied (corresponding to the plotted ratio being less than 1) in most of the parameter space not ruled out by the 8 TeV LHC A/H → τ + τ − bound.  FIG. 9: Vacuum metastability requires this ratio to be approximately less than 1 [47], so we see that most of our points are compatible with vacuum metastability bounds.