Towards a Supersymmetric Description of the Fermi Galactic Center Excess

We attempt to build a model that describes the {\it Fermi} galactic gamma-ray excess (FGCE) within a UV-complete Supersymmetric framework; we find this to be highly non-trivial. At the very least a successful Supersymmetric explanation must have several important ingredients in order to fit the data and satisfy other theoretical and experimental constraints. Under the assumption that a {\it single} annihilation mediator is responsible for both the observed relic density as well as the FGCE, we show that the requirements are not easily satisfied in many TeV-scale SUSY models, but can be met with some model building effort in the general NMSSM with $\sim 10$ parameters beyond the MSSM. We find that the data selects a particular region of the parameter space with a mostly singlino lightest Supersymmetric particle and a relatively light CP-odd Higgs boson that acts as the mediator for dark matter annihilation. We study the predictions for various observables within this parameter space, and find that searches for this light CP-odd state at the LHC, as well as searches for the direct detection of dark matter, are likely to be quite challenging. It is possible that a signature could be observed in the flavor sector; however, indirect detection remains the best probe of this scenario.


Introduction and Background
Observations indicate that 85% of the matter in the universe is non-luminous and point towards the existence of cold dark matter particles. Determining the identify of dark matter (DM) is one of the most pressing issues before us today. An important method of searching for dark matter is indirect detection, which searches for DM annihilating into a pair of Standard Model (SM) particles, producing an excess of gamma rays or antiparticles. The gamma rays can appear as a distinct line at an energy corresponding to the mass of the annihilating dark matter particles, or as a continuous energy spectrum created by radiation from scattering, hadronization, and decay of the SM final states. Several anomalies with these features have recently appeared in data sets [1][2][3][4]. However, it is unfortunately a fact of life that indirect signals of dark matter can be ambiguous, and a dark matter explanation for many of the observed excesses to-date competes with more mundane astrophysical sources and instrumental backgrounds [5].
The Fermi Gamma-Ray Space Telescope has performed an all-sky survey of gamma rays [6] and has searched for dark matter signals from the galactic center, as well as from faint dwarf galaxies. As data accumulates, the latter pursuit has proven to provide relatively background-free constraints on dark matter candidates [7]. However, an excess of gamma rays from the galactic center (GC) has become increasingly intriguing as its characterization has improved [3,4]. The spectrum and angular distributions of this excess are well described by annihilating dark matter, with an energy spectrum peaking at 1-3 GeV, and an approximate spherical distribution about the center of the Milky Way. The flux displays the characteristics of a NFW halo profile of ρ ∼ r −γ with γ = 1.1 − 1.3. A fit to the spectrum is consistent with 30-40 GeV dark matter particles annihilating dominantly into bb final states (a lighter 7-10 GeV particle annihilating to taus is also possible, but somewhat disfavored), with a cross section of σv = 1 − 2 × 10 −26 cm 3 s −1 . This value of the cross section is tantalizingly consistent with that expected for a Weakly Interacting Massive Particle (WIMP) thermal relic, which is a particularly interesting DM candidate since it can account for the observed DM relic density through the simple mechanism of thermal freeze-out. Accordingly, several ideas have recently been put forward to explore a particle physics explanation of this excess in the center of our galaxy [8][9][10][11].
Given the significance of the galactic center excess observed by Fermi, and the apparent match of its characteristics to WIMP dark matter expectations, it is important to explore whether or not it can be accommodated by a UV-complete theoretical framework, such as Supersymmetry (SUSY). SUSY with R-Parity conservation is the most appealing and widely examined scenario for physics beyond the SM and naturally incorporates a DM candidate, the lightest Supersymmetric particle (LSP). Within the framework of conventional thermal WIMP relics, the LSP is generally identified with the lightest neutralino, as is the case in the Minimal Supersymmetric SM (MSSM), although other possibilities exist in more extensive structures. The purpose of this paper is to determine whether the present-day pair annihilation of a neutralino LSP in the galactic center can be the source of the observed Fermi excess. Our results show that this exercise is not trivial, and that specific requirements are placed on model building efforts which attempt to realize this possibility within the SUSY framework. In particular, it appears that fermionic dark matter annihilating via a relatively light pseudoscalar mediator is the best option, as recently discussed by others [10,11].
As part of our model building efforts, we first examined a number of SUSY scenarios that failed to reproduce the GC excess. The 19-parameter p(henomenological) MSSM [12] does not have enough freedom to explain the GC excess, assuming a thermal history for the LSP. A model with Dirac gauginos [13] requires t-channel exchange of a light charged sparticle, in conflict with LEP, TRISTAN and LHC bounds. Extending both the matter and gauge content of SUSY within a GUT context, such as E 6 [14], leads to an annihilation cross section that is far too small. As each of these cases teaches a useful lesson, we briefly summarize our findings for each case in Section 2.
After investigating these options, we examine the Next-to-MSSM (NMSSM), which incorporates an additional gauge singlet superfield [15]. This model is attractive in its own right as it naturally generates the µ-term in the superpotential and easily accommodates the mass of the observed Higgs boson [15]. The extra superfield results in an additional neutralino -the singlino -and two extra Higgs fields, one CP-odd and one CP-even. The extra field content allows a wider range of dark matter properties and annihilation channels [16]. The most common form of this model incorporates a Z 3 symmetry that eliminates terms with a mass dimension in the superpotential. With this Z 3 symmetry imposed, we find that pseudoscalar exchange alone is not sufficient to describe the observations. A sliver of parameter space with significant annihilation through both the Z boson and light pseudoscalar a channels has recently been shown to accommodate [11] the Fermi galactic center excess (FGCE). We summarize our results for the Z 3 option in Section 3.
We next expand our quest to one last SUSY model -the general NMSSM, which incorporates five additional parameters over the Z 3 case. We find that the FGCE is well-described and selects a well-defined region of the parameter space within this model. This region is found to be consistent with other measurements and we study some of its predictions for future experiments. Our results are discussed in Section 4.
In Section 5 we present our summary, where we conclude that the most robust SUSY scenario describing the FGCE is the general NMSSM. Within this scenario, we find that the data selects a specific region of parameter space and that the signatures, such as in Higgs physics at the LHC and in direct detection, are small and will be difficult, if not impossible, to detect. It is possible that indirect detection thus provides the best probe of this scenario.

Initial Model Building Attempts -Lessons Learned
Here, we provide a brief overview of the set of well-motivated SUSY models that we investigated in light of the FGCE. All of these scenarios fail to describe the observed excess while retaining consistency with other measurements. We believe that it is instructive to understand why each case falls short of explaining the excess and will incorporate these lessons in constructing the successful scenario presented below.

pMSSM
The pMSSM [12] is a natural place to begin the search for a plausible SUSY scenario that can explain the FGCE since it contains the minimal particle content, yet enjoys a reasonable amount of freedom in its parameter space. Here, we consider a version of the pMSSM with a Majorana neutralino LSP (χ 0 1 ) as part of a 19-dimensional subspace of the full ∼100parameter R-Parity conserving MSSM. A set of experimentally-motivated constraints have been imposed: (a) The soft-breaking parameters are taken to be real so that no new sources of CP violation are present, (b) MFV holds at the TeV scale so that flavor physics is essentially controlled by the CKM matrix, (c) the first and second generation sfermions are assumed to be degenerate for each sparticle type, and (d) the A-terms and Yukawa couplings of these two generations are assumed to be negligible. In recent work [17], we have generated two sets of pMSSM models (i.e., unique points in this parameter space) with a neutralino LSP. These samples are consistent with collider, flavor, cosmological and precision measurements. We employ the "low-FT" model set, which has a large sample of LSPs with masses between 30 and 40 GeV, for the present study. Other than increasing the number of light LSPs, the fine-tuning requirements placed on this model set have no effect on the properties of LSPs within this mass range.
It is challenging for a ∼ 30−40 GeV neutralino LSP to generate the observed dark matter thermal relic density in the pMSSM. Pure winos or Higgsinos of this mass annihilate very efficiently and yield extremely small values for the relic density; they are also excluded by LEP data since the LSP will necessarily be accompanied by a corresponding nearby charged state, i.e., a chargino. Pure binos require a t-channel exchange of a light sfermion, which is frequently in tension with LEP or LHC bounds on light charged sparticles. A mostlybino LSP with a mass of ∼ 35 GeV can annihilate through an s-channel Higgs (or H/A) boson if it contains a significant Higgsino component. However, the Higgsino component also controls the coupling to the Z boson, so that Z exchange is dominant in the mass region of interest (large couplings of the LSP to the Higgs are also constrained by null results from direct detection experiments). The LSP coupling to the Z boson is constrained by limits on the spin-dependent dark matter direct detection cross section, and by measurements of the invisible width of the Z (∆Γ inv 2 MeV). Interestingly, these combined constraints exclude a mixed bino-Higgsino LSP with a mass below ∼ 31 GeV.
We now examine the annihilation cross section for a LSP with a bino-Higgsino admixture. The top panel in Fig. 1 displays the present-day value of < σv > for the LSPs in the low-FT model set. All of these models predict the correct LSP abundance, and we have selected only models with LSP masses in the interesting range. The panel shows the estimate for this thermal cross section including only the Z exchange contribution, as well as the more complete calculation as performed by DarkSUSY [18]. We see that the Z exchange contribution is a reasonable estimator of the total annihilation rate in most cases, and predicts an annihilation rate which is 2 orders below the value of < σv > during freeze-out (and therefore below the value needed to explain the FGCE). Although some points have significantly larger annihilation rates than predicted by the Z exchange diagram, these models annihilate almost exclusively to τ + τ − through t-channel exchange of light staus. Bottom squarks light enough Figure 1: (Top) Present-day value of the thermally-averaged DM annihilation cross section, < σv >, in units of cm 3 s −1 , plotted against the LSP mass for models with Ω LSP = Ω DM . The black dots show < σv > assuming only s-channel Z exchange, calculated at tree level. (Bottom) σv for χχ → bb as a function of the velocity v χ , mediated by the SM Z boson for m χ = 35 GeV. Here, for purposes of demonstration, the coupling of the LSP to the Z boson, which is expressed through its Higgsino content via the combination of the mixing parameters for the neutralino mass matrix (|Z 13 | 2 − |Z 14 | 2 ), is taken to be 10 −4 ; the cross section can be appropriately rescaled for other values of interest.
give comparable annihilation rates are strongly excluded by collider limits. The present-day annihilation cross section decreases with increasing LSP mass because heavier LSPs are able to annihilate through the Z pole in the early universe, reducing the Higgsino content required to achieve the correct relic density. Since the Z pole becomes kinematically inaccessible at late times (due to low LSP velocities), the present-day annihilation cross section for heavier LSPs is suppressed by their lower Higgsino content.
In the lower panel of this Figure, we show the behavior of σv (without thermal averaging) for a 35 GeV Majorana LSP annihilating through Z exchange, plotted as a function of the LSP velocity. Since the 35 GeV χ 0 1 annihilates through the Z pole in this scenario, the cross section exhibits a strong velocity dependence. In the early universe, the LSP velocity is sufficient to push the collision center-of-mass energy towards the Z pole, and the annihilation cross section can take full advantage of the resonance to generate the observed relic density. In the present era, the LSPs are moving very slowly so that their annihilation takes place away from the Z pole with a cross section that is correspondingly far smaller. LSPs with a slightly lower (higher) mass require a somewhat larger (smaller) coupling to the Z, since they receive, on average, less (more) of a boost from the pole. At present-day velocities, the LSPs are highly non-relativistic, and the kinematic enhancement of the annihilation rate for LSPs close the Z mass (but with a mass difference larger than the Z width) essentially disappears.
We have learned two lessons from this exercise: (i) The FGCE cannot be explained in the pMSSM. Given the parameter freedom in the pMSSM, this result is likely to hold for any model in the MSSM with a Majorana neutralino LSP. (ii) Since the neutralinos are Majorana fields and are required to annihilate through the Z, the cross section has a very strong velocity dependence due to the absence of vector couplings as well as the proximity of the Z pole (M Z /2m χ ∼ 1.3). This renders it difficult, if not impossible, to generate the observed relic density while explaining the FGCE. Clearly it would be advantageous to avoid this strong dependence of < σv > on the velocity by building a model wherein the SM Z resonance is not the dominant channel for both early universe and present-day DM annihilation.

Dirac Gauginos
Another possible scenario is the MSSM with Dirac gauginos, which has recently gained in popularity [13] in light of the null LHC SUSY search results. In this scenario, a ∼ 35 GeV LSP must be almost pure bino (to the level of ∼ 10 −4 or smaller) since any significant Higgsino component will produce a vector coupling of the LSP to the Z boson. This, in turn, yields a very large spin-independent (SI) direct detection cross section [19] which is excluded by experiment. A Higgsino admixture is thus not allowed, ensuring that the annihilation of the LSP does not proceed through the Z boson, which avoids the problems encountered above in the pMSSM.
Essentially pure Dirac binos can have a significantly larger t− and u−channel induced co-annihilation cross section via sfermion exchange than in the corresponding case of a Majorana LSP, even when the sfermion mass satisfies the LEP bounds. The behavior of this cross section is easily understood once we recall that both the conventional µ and Aterms are absent in R-symmetric models, so that left and right sfermions do not mix. Since the sfermions are pure electroweak eigenstates, the cross section has two contributions -one from each of thef L,R exchanges -and the corresponding rates will be proportional to the combinations N c Y 4 L,R , with Y L,R being the relevant hypercharges of thef L,R . Since both sbottoms have a smaller hypercharge than do staus (by a factor of 3), annihilation to the bb final state will be suppressed by a factor of 27 relative to annihilation to τ + τ − when the stau and sbottom masses are equal (although we note that collider searches allow for staus to be lighter than sbottoms). Specifically, in the v 2 → 0 limit, one obtains the cross section [19] for the process χχ → ff : where β 2 f = 1 − m 2 f /m 2 χ and g 1 is the usual SM U (1) hypercharge coupling. This expression shows the expected behavior. Despite the hypercharge suppression, sbottom exchange (and therefore the bb final state) can dominate if sbottoms are much lighter than the other sfermions. The problem now becomes one of rate since collider searches [20] constrain sbottoms to be somewhat heavy, particularly for the light LSP masses required to explain the FGCE, with a limit of m(b 1 ) ≥ 300−400 GeV being extremely conservative. Taking m χ 35 GeV and both sbottoms degenerate with masses at their smallest possible value of 300 GeV, we obtain an annihilation cross section that is too small, by a factor of almost 10 3 , to explain the FGCE or to obtain the observed relic density. (Note that the annihilation cross section during freeze-out is ∼ 10 − 15% larger than its present-day value due to the O(v 2 ) terms that were neglected in the above expression.) Explicitly, in the v → 0 limit we obtain the result: σv 3.7 × 10 −29 cm 3 s −1 . The annihilation rate would decrease further if heavier, more realistic sbottom masses were employed. Clearly this scenario cannot be responsible for the galactic center excess.

SO(10)/SU (5) Singlets in E 6 GUTS
Another approach is to consider a SUSY Grand Unified Theory (GUT), extending both the gauge group and matter content of the MSSM. In this case, one of the new weakly interacting, SM-singlet neutral fermions can be the LSP and potentially a DM candidate. One attractive possibility is the SUSY version of E 6 , which can lead to an additional U (1) gauge group at the TeV scale [14]. Here, the GUT-scale breaking pattern can take the form of E 6 → SO(10) × U (1) ψ , followed by SO(10) → SU (5) × U (1) χ and then, as usual, SU (5) breaks to the SM. A linear combination of the two additional U (1) groups, U (1) θ = U (1) ψ cos θ − U (1) χ sin θ, survives to lower (i.e., TeV-scale) energies. In addition to the extended gauge group, the standard MSSM matter and Higgs content is also enlarged as both are contained in the fundamental 27 (and possibly 27) representation of E 6 . The 27 decomposes as 16+10+1 under SO(10) with the SM fermions and a RH-neutrino composing the usual 16. Within this scenario, one of the new light neutral SM singlet fields provides an appealing candidate for the LSP; a particular realization of this is given by the choice of the SO(10)-singlet field, S. Since the LSP is a SM singlet, it must annihilate through the new Z gauge boson. The mass of the Z is required to be quite heavy, M Z 2 TeV, by null results from LHC resonance searches [21]. Although the heavy Z mass avoids issues associated with direct detection, the invisible width of the Z, or the significant velocity dependence observed above in the case of the pMSSM, we will see that it also leads to a problematic suppression of the annihilation rate. Another drawback of annihilation through the Z is that bb is not necessarily the dominant final state for LSP annihilations. Depending on model details, the LSP may be either a Majorana or Dirac fermion. Once this choice is made, the model contains only two free parameters, the Z mass and the value of the mixing angle θ, and is therefore very predictive. Once these parameters are known, the freeze-out and present-day annihilation cross sections are easily determined. TeV as a function of the mixing angle θ. The red (blue) curve is the result obtained for the SO(10) singlet field S when it is Majorana (Dirac), while the green (cyan) curve represents the corresponding result for the SU (5) singlet within the SO(10) 16 representation. Figure 2 shows the thermal relic cross section to O(v 2 ) as a function of the mixing angle θ for 35 GeV Majorana or Dirac S LSPs. For our numerical results, we have assumed that x = m χ /T f = 25, with T f being the freeze-out temperature, and have taken M Z = 2 TeV, saturating the lower limits from the LHC. Note that in this case the expansion of the thermally-averaged cross section in powers of velocity is fully valid as the annihilation takes place far away from the Z pole. We see that the cross section vanishes in the limit when the Z is purely a Z χ state as S is a SO(10) singlet. Whether the LSP is Dirac or Majorana, we see that the resulting cross section is far too small (by at least several orders of magnitude) to predict the correct relic abundance or explain the FGCE. This results from the rather small gauge couplings in this class of models and from the requirement of a heavy Z paired with a light LSP, since the cross section roughly scales as ∼ m 2 χ /M 4 Z . Although the annihilation rate is much larger for a Dirac LSP than for a Majorana LSP, both results are far too small to produce the observed relic abundance or explain the FGCE. We obtain a qualitatively similar result in the case were the LSP transforms instead as the SU (5) singlet within the SO(10) 16 representation, although the detailed dependence of the cross section on θ is quite different and the cross section is now seen to vanish near θ −14.48 o as can be observed in Fig. 2. Once again, the Dirac LSP case leads to a much larger cross section but also falls far short of what is required.

The NMSSM
In the Next-to-Minimal Supersymmetric Standard Model (NMSSM), the Higgs sector is extended by an additional gauge singlet superfield,Ŝ. The scalar component of this superfield obtains a vev, s, and subsequently generates the usual MSSM µ parameter dynamically through a term λĤ uĤdŜ in the superpotential, i.e., µ ef f = λs. The presence of the fieldŜ also implies the existence of an additional neutralino (the singlino), as well as one additional pair of weak isosinglet, CP-odd and CP-even Higgs bosons that mix with their usual MSSM counterparts. This results in two physical pseudoscalar bosons, A 1 = a and A 2 = A and three physical neutral CP-even Higgs bosons, H 1 = h, H 2 , and H 3 , where a and h are the lightest CP-odd and CP-even states, respectively. Here we assume that h ∼ 125 GeV is the SM-like Higgs boson discovered at the LHC. We identify the ∼ 35 GeV LSP as being principally composed of the singlino, with small mixings with the other gauginos. The LSP can annihilate to heavy fermions (mostly bb) by exchanging a relatively light pseudoscalar a, which couples to SM fields through its mixing with the MSSM pseudoscalar A. Reasonably large values for tan β are generally required to provide a large enough abb coupling. In this discussion, we assume that s-channel a exchange is solely responsible for producing both the observed relic density and the FGCE, although this need not be the case [11]. With this premise, the LSP coupling to the SM Z boson must be extremely weak in order to avoid Zmediated LSP annihilation; this requires the Higgsino content of the LSP to be quite small. In principle, such a 'simple' scenario appears to avoid [10] all of the pitfalls we have discussed in the previous section, since the LSP pair annihilation is no longer velocity suppressed, while the SI direct detection cross section has a very strong velocity suppression [9]. The flux from the galactic center appears to require, if anything, a present-day annihilation rate which is slightly smaller than the early universe freeze-out value, implying that a must be heavier than 2m χ 70 GeV so that the average center-of-mass energy of the annihilation is closer to the a pole at the larger DM velocities present in the early universe. This assertion further implies that the process h → aa (where as noted above h is the 125 GeV SM-like Higgs boson discovered at the LHC) cannot occur on-shell. Since the difference between the early universe and present-day annihilation rates is small, m a must be large enough so that early universe annihilations (with v 2 0.10 − 0.15) don't benefit too much from the a pole enhancement after thermal averaging. On the other hand, if m a becomes too large, the annihilation cross sections will be too small to account for the observed relic density or the FGCE, indicating that a balance between these constraints is necessary. As we will see below, satisfying all of these conditions simultaneously is quite challenging even in the NMSSM.

The NMSSM with Z 3 Symmetry
The simplest and most frequently studied version of the NMSSM includes a Z 3 symmetry that eliminates all terms with mass dimension from the superpotential. With µ ef f generated as described above, the only other term [15] in the superpotential not found in the MSSM is κ 3Ŝ 3 . As a result, only a few parameters are required beyond those found in the MSSM: λ, κ, s and A λ,κ , with the latter being the corresponding soft breaking A-terms, i.e., and, similarly, Interestingly, it is easy to argue that even at tree level these few parameters alone do not allow enough flexibility to construct a successful model along the lines discussed above where only the light pseudoscalar a participates significantly in DM annihilation. This can be most easily seen by examining the tree-level neutralino mass matrix in the Z 3 -invariant case as is given in Eq.(2.32) of Ref. [15] which we repeat here for clarity: Here, v u (v d ) is the vev of the Higgs doublet giving mass to the up (down) sector, µ = 0 in the Z 3 invariant case, and g 1,2 are the usual SM gauge couplings. Recall that the effective µ parameter is given by µ ef f = λs.
There are several constraints that need to be imposed on this neutralino mass matrix. First, the mixing of the mostly singlino LSP with the Higgsino must be highly restricted, otherwise the LSP will have significant couplings to the Z boson (violating our assumption that only a is responsible for DM annihilation) and all of the problems present in the pMSSM will be encountered. Using the results in Fig. 1, we can easily determine the present-day value (v → 0) of the cross section for LSP annihilation through the SM Z boson. Specifically, we obtain the result σv 3.39 × 10 −32 cm 3 s −1 when the combination of mixing parameters for the neutralino mass matrix (|Z 13 | 2 − |Z 14 | 2 ) 2 = 10 −5 and m χ = 35 GeV are assumed, where Z ij is the matrix that diagonalizes Eq. (4). (Of course, for other possible Higgsino admixtures these numerical results are simply rescaled.) Similarly, employing x = m χ /T 20 (25) in order to calculate the thermally-averaged cross section in the early universe, we obtain the significantly larger result of σv 6.56(4.67) × 10 −30 cm 3 s −1 under the same assumptions since the Z pole is partially probed by the thermal LSP energy distribution. Here we see explicitly that the cross section in the early universe is significantly larger than it is today by a factor of over 100. We thus require that a exchange remains dominant over Z exchange at all times, and impose the constraint that (|Z 13 | 2 − |Z 14 | 2 ) 2 ≤ 10 − (3−4) . Second, the Higgsino mass, given by µ ef f = λs, must be > 100 GeV due to LEP constraints on the mass of the corresponding charged state, while the parameters κ, λ must be less than 0.7 to satisfy perturbativity requirements. Finally, the LSP mass in this scenario is roughly given by m LSP 2κs 35 GeV, constraining the value of κ.
At this point, we note that the DM annihilation rate through the a boson is essentially set by the combination of couplings √ 2κ sin θ a cos θ a Z 2 15 , where θ a measures the weak doublet content (i.e., the MSSM A content) of a and Z 2 15 1 measures the singlino content of the LSP. This implies that the value of κ cannot be too small, otherwise the calculated relic density will be too large. We can see this explicitly by examining the annihilation cross section. Before any thermal averaging, we find that σv for the process χχ → a → bb is given by the expression where with θ a being the CP-odd mixing angle as described above, and λ χ √ 2κ cos θ a Z 2 15 /g controlling the aχχ coupling [15]. To get an idea of the size of this cross section let us assume that the product of couplings Q = λ b λ χ = 1 and rescale the result later as needed; the resulting cross section is shown in Fig. 3 as a function of v χ for several different input parameter choices, as labeled.
From this Figure we learn several crucial lessons: (i) the cross section has a rather simple sensitivity to both tan β and m a , (ii) a lower (higher) value of m a increases (decreases) the sensitivity to the a resonance at large values of v χ , and (iii) the result is not sensitive to the precise value of Γ a /m a as long as this value is small and the center-of-mass energy of the annihilation is not too close to the resonance. To be specific for this scenario, we choose m a = 120 GeV, tan β = 40, and Γ a /m a = 0.01, although the latter choice has very little impact as seen above. These parameter choices yield a present-day thermallyaveraged cross section, < σv >= 7.56 × 10 −24 cm 3 s −1 , with an early universe freeze-out value of < σv >= 9.05(8.72) × 10 −24 cm 3 s −1 with x = m χ /T = 20 (25), assuming that Q = 1. This result strongly constrains the value of Q, requiring Q ∼ 1/(300 − 400). For values of m a 120 GeV we learn from [10] that sin θ a 0.1 is likely necessary to satisfy experimental constraints when m A ∼ 1 TeV; assuming that sin θ a ∼ 0.1 we see that κ cannot be too small, with values of κ ∼ 1/4 − 1/3 preferred, if Q is to be large enough to explain both the FGCE and the observed relic density.
It is instructive to examine the sensitivity of both the present-day and the early universe thermally-averaged cross sections to variations in the values of both m χ and m a in the region of interest; this exercise will yield information about the possible values for the remaining model parameters, such as κ. The results of these calculations are shown in Fig. 4, again taking the factor Q to be unity. Here we see that these cross sections are reasonably sensitive to both of these mass parameters, varying by approximately an order of magnitude over the ranges of interest. The Q-independent ratio of the present-day thermally-averaged cross section to its early universe value is somewhat less sensitive to these parameters, but displays interesting mass dependencies. This ratio is less than unity by construction since m a > 2m χ as discussed above. For fixed m a , annihilating LSPs of increasing mass start to probe the kinematic region near the a pole with their thermal energy distribution, decreasing the ratio of the present-day annihilation cross section to its early universe value. For fixed m χ , increasing m a decreases the sensitivity to the a pole region and the ratio tends towards unity. For most of the parameter values this ratio is seen to be roughly ∼ 0.8 with a spread of approximately ∼ 10%, which agrees with the FGCE observations.
Although we can avoid a troublesome suppression of the present-day annihilation rate, several model building requirements combine to restrict the parameter space in this scenario. We start by considering the neutralino mixing matrix more carefully. As mentioned above, we want the LSP to be mostly singlino in order to have a large coupling to the pseudoscalar singlet. In this case, the bino and wino components of the LSP are typically small since bino-singlino and wino-singlino mixing terms are absent from the neutralino mass matrix, and because LEP results require the charged wino to be heavier than 100 GeV. We therefore assume that the bino and wino are irrelevant and consider only the Higgsino-singlino sector, which is completely specified by 4 parameters (s, κ, λ and tan β). Recalling that the LSP must have a small Higgsino component to suppress the Zχχ coupling, we see that for large values of tan β, λ must be small enough that the off-diagonal entries are much smaller than the diagonal entries, µ ef f = λs and 2κs 35 GeV. Specifically, the mixing is roughly set by the ratio 8λ 2 v u v d /(λs + 2κs) 2 . LEP limits require µ ef f ≥ 100 GeV, which demands a large value of s in the (necessary) case of small λ. Since the LSP mass is set by 2κs m LSP 35 GeV, large values of s demand small values of |κ| to produce the correct LSP mass. This is in tension with the the requirement that |κ| should be large enough to predict the correct annihilation cross section. This tension indicates the failure of this approach. While it is likely that this complete argument can be demonstrated algebraically, we have instead performed large numerical scans (totaling ∼ 10 10 points in the space) over the 4 parameters relevant for singlino-Higgsino mixing in order to verify our conclusion. No successful solutions were obtained, and we therefore conclude that the NMSSM with Z 3 symmetry where only the a exchange is responsible for both the DM relic density and the FGCE also fails to reproduce the observations. , and the ratio of these annihilation cross sections (bottom), which is Q independent. In the top panels, from top to bottom m a ranges from 110 to 130 GeV in steps of 2 GeV. This mass ordering is reversed in the bottom panel.

The General NMSSM
To go further we drop the requirement of the Z 3 symmetry on the NMSSM superpotential, providing greater parameter freedom. Although the scale-invariant feature of the Z 3 model is very attractive, it also poses potential domain wall problems [22] that we now avoid. Additionally, some of the attractive features of the Z 3 model can be recovered or improved upon in the general NMSSM with model-building effort [23]. In order to most easily discern whether this more general scenario can provide a successful description of the FGCE, we first present a (mostly) tree-level discussion, which is important for providing the required transparency given the relative complexity of the parameter space. A more detailed numerical analysis including all corrections and a full scan over the parameter space follows in the next section.
The number of additional parameters is significantly increased in the full NMSSM: First, the superpotential is augmented by terms which are either linear or quadratic in the singlet superfieldŜ. Second, corresponding additional terms are necessarily generated in the softbreaking Lagrangian. Thus we find, instead of the compact expressions above, the augmented results [15] (ignoring a possible bare µ term) and From these expressions we see that the general NMSSM scenario introduces 5 additional new parameters (µ , ξ F , ξ S , m 3 , m S ) beyond those encountered in the Z 3 model; only a subset of these additional parameters is required (at least at tree level) to describe the FGCE. In the brief discussion that follows, we assume for simplicity that m 2 3 = ξ F = 0 in order to derive a workable model, but this need not be the case; viable models can be found when these parameters are non-zero and we include them below in our full scan over the parameter space. As was shown above, the neutralino mass matrix depends on only one of these new parameters, µ , which is sufficient to resolve the difficulties we found in the case of the Z 3invariant NMSSM. To see this, we again assume that M 1,2 are large so that the bino and wino are decoupled, and similarly assume that the value of µ ef f is large enough to suppress Higgsino-singlino mixing while the value of λ is somewhat small. This implies that large values of s are necessary, meaning that 2κs is also relatively large if κ is not too small as previously discussed. We can, however, choose µ to have a similar magnitude to 2κs with the opposite sign, so that these terms cancel to produced a small remainder, i.e., we assume the combination 2κs + µ 35 GeV. Note that a fine-tuning at the few percent level is thus required to obtain the desired LSP mass. We then perform a simple (flat) scan over parameters relevant for singlino-Higgsino mixing, employing the scan ranges 0 ≤ λ ≤ 0.6, 30 ≤ tan β ≤ 50, 1 ≤ s ≤ 10 TeV with 0 ≤ κ/g 2 ≤ 2 and imposing the requirements discussed above. We obtain multiple valid solutions, demonstrating the viability of this approach.
Given values for the parameters λ, κ, tan β, s, µ resulting from our scan that yield a neutralino mass matrix with the desired properties, we next consider the effective 2 × 2 CPodd scalar mass matrix at tree level. We use the convenient basis where the conventional neutral Goldstone boson has been removed so that the elements of this matrix are given by [15] where the combinations are defined as Note that in the absence of mixing M P, 11 gives the mass of the usual CP-odd field A in the MSSM. Recall that for the purposes of this discussion, we have assumed m 2 3 = ξ F = 0, and now, for further simplicity, we will also assume that the combination 2m 2 S + ξ S s = 0. With these assumptions, the only new parameters present in this matrix but irrelevant for Higgsino-singlino mixing are the A-terms, A λ and A κ . With these assumptions, we can set θ a = 0.1 and fix the smaller eigenvalue of this matrix to be m 2 a (120GeV) 2 ; we then obtain numerical values for the soft parameters A λ,κ and the mass of the heavier CP-odd state A. Since sin θ a is small, this state is similar in nature to the MSSM pseudoscalar boson. It is therefore constrained by searches for resonant di-tau production [24,25], and also potentially by flavor constraints [26] such as those arising from measurements of the branching fraction BR(B s → µ + µ − ), for the large values of tan β (∼ 40) considered here. We therefore avoid these constraints by requiring M A > 1 TeV. This approach works quite well and multiple solutions are again easily obtainable.
Lastly, we turn to the CP-even Higgs mass-squared matrix, where we provide the treelevel elements here for completeness [15]: Here g 2 ≡ (g 2 1 +g 2 2 )/2 with g 1,2 being the usual SM gauge couplings. In order to make contact with experiment we must incorporate the important large stop/top and sbottom/bottom 1and 2-loop radiative corrections to this matrix as given in [15]; this is particularly important for the 11, 22 and 23 elements. Given our assumptions outlined above, m 2 S is the only new free parameter entering this matrix that we need to scan. Here we employ the range 0 ≤ m 2 S ≤ s 2 , which is specifically chosen to avoid any tachyonic solutions for the eigenvalues of this matrix, for our scan. We require the lightest eigenvalue, which we identify as the Higgs boson observed at the LHC, to have a mass m h = 125.5 ± 3 GeV. We further demand that this state must have very little mixing with the singlet Higgs, i.e., |S 13 | 2 << 1 in the notation of Ref. [15]. This requirement ensures that the observed Higgs state will have properties that are very close to that of the MSSM and SM Higgs and, in particular, will have a small branching fraction for the decay h → χχ (generally ∼ 1%), which can only proceed through the various neutralino and Higgs mixings. This also reduces the observed Higgs boson's contribution to the Spin-Independent Direct Detection cross section for the almost pure singlino LSP.
Performing a scan over the parameters λ, κ, t β , s, m 2 S yields many sets of parameter values (i.e., 'models') that produce completely successful phenomenology. Certainly allowing the other remaining parameters that we have ignored in our discussion to be non-zero or to float freely will also yield viable scenarios. Thus, at least at (mostly) tree level, the general NMSSM can provide a phenomenologically viable description of the FGCE!

Full Scan of the General NMSSM
To go further in our exploration of the NMSSM we perform a more thorough and complete analysis incorporating, e.g., the higher-order, loop-generated radiative corrections to the various gaugino, CP-even and CP-odd Higgs mass matrices as well as the associated couplings in the NMSSM. In performing this analysis we make use of a modified version of NMSSMTools4.3.0 [27, 28] 1 and perform a scan over the general NMSSM parameter space. We subject each of the randomly chosen parameter space points to the various experimental and theoretical constraints mentioned above, along with additional constraints which will be detailed in the next section, in order to find phenomenologically successful models that describe both the observed relic density and the FGCE.
As was noted in the previous discussion, the general NMSSM contains 10 new parameters beyond those already present in the MSSM which we now include as part of our scan. An outline of our scan procedure is as follows. For simplicity, we choose fixed values for the usual gaugino masses (M 1,2,3 ), the slepton soft masses, and the leptonic A-terms, and scan over the ratio of the two Higgs vevs, tan β. Additionally, we set the NMSSM parameters m 2 S = m 2 3 = 0 since they play no important role in this analysis. Motivated by the FGCE measurements and our previous discussion, we are specifically interested in particular ranges for the physical masses of the LSP, the lightest CP-odd scalar and the lightest CP-even scalar (which we equate with the ∼ 125 GeV Higgs as discussed above). We also want to constrain the mass of the heavier CP-odd field from below to satisfy constraints from LHC searches and flavor physics. We therefore eschew scanning over the fundamental Lagrangian parameters µ and ξ F,S in favor of a scan over the interesting ranges of the physical masses m a , m A , m h and m LSP . The (common) squark masses and the corresponding squark Aterms are also set by choosing values for these physical masses, since suitably large radiative corrections are necessary in the NMSSM (as in the MSSM) in order to obtain a Higgs mass in the ∼ 125 GeV range. After performing this substitution, the remaining scan parameters are just λ, κ, µ ef f and the soft parameters A λ,κ . The approach that we have just outlined mimics the procedure that we followed in our tree-level discussion as closely as possible.
We now discuss our scan parameters and their ranges in a bit more detail. Some of the parameters are more relevant to our aim of demonstrating the existence of points that can explain the FGCE than others. We therefore divide the set of parameters into "fixed parameters," which are not expected to have a significant effect on the relevant physics and hence are set to reasonable values, and "scanned parameters," which are scanned using a uniform distribution (i.e., employing flat priors) within the ranges specified below. In NMSSMTools, these parameters are defined at the scale )/4 and are then run to Q 2 = mQ 3 mũ 3 (with the exception of tan β, which is input at Q = M Z ). Since we choose a common squark mass parameter, mQ = m (ũ,d) 1,2,3 , for our analysis, these scales are both equal to mQ. These scales are slightly different from the geometric mean of the stop masses, which is often employed as the SUSY scale, but we do not expect this difference to have any effect on our results. Fixed Parameters: As the gauginos do not play an important role in our scenario, provided they are sufficiently heavy to evade any experimental constraints, we set Similarly we fix the slepton masses to be and the third generation lepton trilinear coupling is set to We do not expect our results to be particularly sensitive to any of these choices. Additionally, in the general NMSSM, one could assign non-zero values to the Lagrangian parameters m 2 3 and m 2 S . However, as noted above, we have fixed these values to zero at the input scale. Scanned Parameters: The remaining Lagrangian parameters are chosen via a scanning procedure. These parameters, together with the ranges from which they were chosen using uniform distributions (flat priors), are listed in Table 1. Note that we perform a much broader scan of the parameters than for our tree-level study described above, since we are now interested in characterizing the phenomenologically interesting region in addition to demonstrating its existence.
As noted in Table 1, we have substituted several of the Lagrangian parameters with the corresponding physical parameters m h , m a , m A , and m χ 0 1 . This increases the efficiency of our scan since we are only interested in a restricted range of values for m h , m a , and m χ 0 1 .

Value Lower Bound Upper Bound
Replaced -- Table 1: Summary of how the parameters of the general NMSSM were chosen in our scan. The "Value" column gives the parameter value for fixed parameters. A "Scanned" entry in this column indicates that the parameter value was chosen randomly from a uniform distribution within the range specified by the subsequent two columns. A "Replaced" entry indicates that the parameter value is solved for numerically to obtain desired values of the physical parameters listed in Table 2.
In addition, m A is generically strongly constrained by experiment and is not allowed to be small. This efficiency boost is particularly necessary since we only obtain acceptable models at somewhat finely tuned regions of the parameter space, with a significant cancellation between µ and 2κs (= 2κµ eff λ using our scan parameters). As we are interested in parameter space points for which both of these terms are O(TeV), yet 30 GeV < |m χ 0 1 | < 40 GeV, it is obvious that scanning over all the parameters κ, λ, µ eff , and µ will be inefficient.
In order to proceed, we first set all squark mass parameters to a common value: as noted above, and the third generation squark trilinear couplings to be proportional to the squark mass parameter  Table 2: Ranges employed in our scan for the physical values of the listed parameters. Here m h is the mass of the lightest CP-even Higgs boson, m a,(A) is the lighter (heavier) CP-odd Higgs boson, and m χ 0 1 is the signed mass of the lightest neutralino, which is also the LSP. We scan over both signs of the LSP mass. Again note that these enlarged ranges allow for a full exploration of the phenomenologically interesting parameter regimes. and ξ S with the variables MP and MA, which represent the first and second diagonal elements, respectively, of the effective 2 × 2 CP-odd Higgs mass matrix (i.e., with the Goldstone boson removed) at tree level.
To begin our scan, we chose target values of m h , m a , m A , and m χ 0 1 , along with values for the remaining scanned parameters. To find the Lagrangian parameters (mQ, MP, MA, µ ) corresponding to these target values, we chose a seed point in the (mQ, MP, MA, µ ) parameter space that would be unlikely to lead to a tachyonic particle in the spectrum. We then used a simple variant of Newton's method to numerically search for a point giving the correct physical masses. In this search method, the 'next step' in parameter space was scaled by a small constant to reduce the chances of obtaining a tachyonic spectrum. If the numerical search for appropriate physical masses failed, a new "safe guess" for the initial parameter values was made and the procedure restarted. However, only a finite number of such "restarts" were allowed. If too many searches ended unsuccessfully as a result of, e.g., any of the particles being tachyonic, or if the search failed to obtain a point with the required values for the physical masses after a large number of iterations, the scan point was rejected.
In performing this general parameter scan a number of issues arose that required modifications of the default NMSSMTools4.3.0 package and that partially influenced the constraints we have applied in exploring the NMSSM parameter space. For example, as noted above, the calculation of Higgs masses in the NMSSMTools package includes radiative corrections to both the scalar and pseudoscalar neutral Higgs mass matrices, as well as to the charged Higgs mass. We find that for much of the parameter space relevant to our study, the radiative correction to the (1, 2) element of the neutral scalar Higgs mass matrix arising from chargino loops lead to spuriously large corrections, resulting in a tachyonic Higgs sector. This is most likely due to a failure to properly resum the large logarithms. In order to remove this unphysical effect in a consistent manner, we neglected the radiative corrections from chargino loops to the scalar and pseudoscalar mass matrices, as well as to the charged Higgs mass, if in any of these matrices the absolute value of the modification of any entry arising from the application of this specific radiative correction was greater than the initial absolute value of that entry. None of the other radiative corrections included in NMSSMTools were altered in any way as in these cases the corrections were under control.
Furthermore, NMSSMTools employs micrOMEGAs [29][30][31][32][33] to calculate the thermal relic density and other dark matter properties. As micrOMEGAs is based on CalcHEP [34], a description of all the vertices in the theory under consideration can be found in the appropriate model file. In this file we found that the Higgs-neutralino-neutralino couplings are written in terms of neutralino masses. This is acceptable in the Z 3 -invariant NMSSM, but in the general NMSSM we found, e.g., an additional contribution to the singlino mass from the µ term, as shown above. This term does not contribute to the Higgs coupling with the singlino. This is extremely important for our purposes as we are considering parameter space regions with a relatively light singlino LSP, but with an unsuppressed singlet-singlino-singlino coupling. Therefore we needed to replace the expressions for the hχ 0 1 χ 0 1 and aχ 0 1 χ 0 1 vertices in the relevant model file with those found in Ref. [15]. In implementing the expressions from [15], we were careful to choose signs such that the new vertex expressions agree with the previous ones in the limit of the Z 3 -invariant NMSSM. To repeat, in order to correct for this error, we have replaced the LSP couplings to the two lightest CP-even and CP-odd scalars in micrOMEGAs with the appropriate expressions from Ref [15]. We note that these couplings determine the dominant contributions to both the indirect and direct detection cross sections in the scenarios considered here.

Theoretical Constraints Applied to the Parameter Space
In performing this analysis, we generated a set of 6·10 8 points in the NMSSM parameter space according to the procedure discussed above and passed them through various experimental and theoretical constraints. Some of the details are worth a short description.
To ensure theoretical consistency of our models, we remove those that are flagged by NMSSMTools as having any of the following problems: (i) Negative mass-squared values for the heavy scalar Higgses, pseudoscalar Higgses, or sfermions, (ii) zero values for tan β, λ, or µ, (iii) numerical problems in the spectrum calculation and (iv) the presence of an unphysical minimum in the scalar potential. In addition to these constraints, we also check for cases where the radiative corrections are potentially non-perturbative. One example of this is the chargino corrections to the Higgs couplings, described above. Another case in which these couplings can become extremely large is in the couplings of the Higgses to bottom quarks. Here, it is possible to generate extremely large radiative corrections when µ ef f is large and negative and tan β is also large. This type of correction is well known, and can in fact be an O(1) effect, as described in [35]. The size of this correction for the scalar φ i can be estimated by taking the ratio Although the tree-level prediction for these r i is approximately 10, a sizable set of models have r i above 100 for one or more Higgses, suggesting that the corrections may need to be resummed. Since this resummation is beyond the scope of this work, we simply discarded models with potentially unphysical corrections. We estimated the range of r i that might result from valid corrections by examining our set of pMSSM models whose Higgs properties are studied in [35], and for which the hbb coupling corrections can be large but are under Observable Experiment SM Prediction Estimated Uncertainty ∆M B d [ps −1 ] 0.507 ± 0.005 [36] 0.545 ± 16.8% [20] [37] Parameter dependent ∆M Bs [ps −1 ] 17.719 ± 0.043 [36] 17.70 ± 6.3% [20] [37] Parameter dependent BR(B s → µ + µ − ) ×10 9 2.9 ± 0.7 [26] 3.74 ± 4.1 % [20] [37] Parameter dependent BR(B → X s γ) ×10 4 3.43 ± 6.7 [36] 3.17 ± 7.3% [20] [37] Parameter dependent BR(B → τ ν τ ) ×10 4 1.14 ± 0.22 [36] 0.779 ± 8.6% [20] [37] Parameter dependent Process dependent [20] 0 σ N P , LHC Process dependent [20] 0 - control. We define the corresponding minimal and maximal values according to the extremal values appearing in the pMSSM model set as follows: .
With these definitions, we find r max = 60.0 and r min = 2.42, which we then use as upper and lower bounds on r i in the NMSSM. We note that only very large corrections to r i modify our results significantly, and that our results are therefore insensitive to the precise values of r max and r min that we choose.

Experimental Constraints Applied to the Parameter Space
To create a selection of viable NMSSM models, we apply the set of experimental bounds listed in Table 3. These can be divided into constraints from flavor physics, Higgs measurements, dark matter observables, and collider searches. The first 5 observables in the Table are flavor observables that are computed by NMSSMTools. This calculation includes an estimate of the theoretical uncertainty, including both parametric errors and uncertainties from missing terms in the diagrammatic expansion. For a given model, the combined uncertainty is the linear sum of the (model-specific) theory uncertainty and the experimental error. We designate a model as being viable if it agrees with the experimental value within this combined uncertainty. Since the experimental ranges in NMSSMTools are supplanted by more recent results for the flavor observables that we consider, we have modified NMSSMTools to use the updated experimental results listed in Table 3. We also note that the theoretical precision of flavor observable calculations within NMSSMTools does not match that of the most precise SM calculations, with the result that the SM value predicted by NMSSMTools can differ from the SM values listed in Table 3. For the purposes of determining which models are excluded by flavor measurements, this fact is included in the estimated theoretical uncertainty. In our discussion, however, we will normalize the NMSSMTools prediction for each observable by the corresponding SM value as predicted by NMSSMTools in order to clearly show the effects which result from BSM physics, and not from e.g. different values for the CKM parameters. Although NMSSMTools computes the NMSSM contributions to the muon magnetic moment, we do not use this observable as a constraint given the unclear status of the relevant theoretical uncertainties and the experimental agreement with the SM. We note, however, that all of the models in our final selection are within the range (−17.7 · 10 −10 − 43.8 · 10 −10 ) for ∆(g − 2) µ , which is the region we allowed in our study of the pMSSM [17].
As noted above, we assume that the LSP is a thermal relic and accounts for all of the observed dark matter abundance. As a result, the relic density measurement from Planck and the null results from direct detection experiments place important constraints on our models. We employ the default treatment of the relic density within NMSSMTools, requiring it to be within 10% of the Planck observation to account for theoretical uncertainties in the relic density calculation (the experimental errors are negligible by comparison). For direct detection, we use the LUX limit on the spin-independent cross section, but rescale it by a factor of 4 to account for matrix element uncertainties in the nucleon, e.g., the strange content of the proton. Including such uncertainties has minimal impact on our results. Spin-dependent experiments are not yet sensitive to models that remain allowed by the other constraints, so we make no explicit requirement on the spin-dependent cross section.
Constraints from observations of the SM-like Higgs boson also set important limits on our model parameters. First, as noted above, we require that the lightest Higgs mass be in the range 122-128 GeV, assuming a 3 GeV theoretical uncertainty in the Higgs mass calculation [41] with the experimental value of ∼ 125 GeV. 3 Secondly, we note that the couplings of the 125 GeV Higgs to SM final states (other than bb) are very close to their SM values in the models we consider. Since the bb channel is presently poorly constrained by direct measurements, the most important effect of increasing (decreasing) the hbb coupling in our model set is to simultaneously reduce (enhance) the Higgs branching fractions to the other SM final states, thus altering the total observed signal strength. This same effect also provides the dominant constraint on the Higgs invisible branching fraction, given the SM-like nature of the light Higgs couplings. Reproducing the CMS measurement of µ = 1±0.26 (i.e., 2σ confidence), where µ denotes the relative strength of the measured Higgs couplings to their SM expectations averaged over all final states, requires that the sum of the branching fractions for bb and χχ final states lies between 0.468 and 0.686. Here, we neglect the direct observation of the bb final state, since the presently large errors on this measurement make this a reasonable approximation.
Finally, we apply constraints arising from the non-observation of SUSY particles at LEP and the LHC. These restrictions include an upper limit on the invisible width of the Z boson from LEP, bounds on superpartner masses from both LEP and LHC searches, and direct searches for scalar resonances decaying to ditau final states at the LHC. Given our scan ranges for the NMSSM parameters, only charginos, stops and sbottoms can a priori be in tension with either the LEP or LHC data. Since the masses of these sparticles are generally unimportant for the dark matter scenario we consider here, we simply avoid these collider searches by requiring charginos to be heavier than 100 GeV and stops/sbottoms to be heavier than 700 GeV. Finally, we make use of the NMSSMTools constraints on the ditau signal strength in the gluon fusion and bb-associated production channels, which are derived by comparing the calculated signal strengths to the limits given in [24].
After applying these experimental constraints, we make the final requirement that our models provide a good explanation of the galactic center excess. As noted above, we require a LSP mass in the 30-40 GeV range (set by our choice of scan range) and then further require that the present-day annihilation rate exceeds 0.7·10 −26 cm 3 s −1 in order to obtain the required FGCE described by Ref. [4]. All of our models will necessarily predict a present-day annihilation cross section which is slightly suppressed compared to that for the early universe as described earlier, and are therefore consistent with upper limits on the annihilation flux.
After all of these constraints are applied, a set of ∼ 52.8k NMSSM model points remain viable, whose properties we now discuss in some detail.

NMSSM Numerical Results and Discussion
We begin by examining the regions for the fundamental parameters of the superpotential and soft-breaking Lagrangian that reproduce the FGCE and are consistent with all other constraints we have applied. Recall that the values of both the superpotential and soft parameters were taken from initially flat distributions. The first three panels of Fig. 5 and all four panels of Fig. 6 display the density of models in a region, as projected onto planes formed by pairs of these fundamental parameters. There are several observations to note: (i) The model density distribution for λ reaches maximum near ∼ 0.15 while κ prefers to take on larger absolute values with its model density being asymmetric with respect to its sign. 4 (ii) For positive values of κ, we note that a reduction in model density appears when λ is large and κ is small; this does not occur for negative values of κ. (iii) The distribution of A κ is quite uniform over the chosen scan range while that of A λ peaks below ∼ 10 TeV and is seen to be sign-asymmetric with positive values being preferred. Note that since A λ always appears together with λ itself in the soft-breaking Lagrangian, their values are correlated so that the product |λA λ | always lies below ∼ 15 TeV. (iv) The density of tan β values is relatively constant above ∼ 10 but on rare occasions models are seen to extend down to as low as tan β ∼ 1.5. For large tan β, the set of models with large negative values of µ ef f are seen to be relatively depleted. (v) Unsurprisingly the values of 2κs and µ are very highly correlated as their sum at tree level produces the rather small LSP mass (which may take either sign). While |µ | tends to have values of several TeV or more (with better than 1% tuning in the LSP mass as described above), there is a long tail stretching down to much lower values, even below 100 GeV. There is, in fact, a single model in our set with |µ | ∼ 1 GeV wherein the LSP picks up a substantial Higgsino content and the mass of the a is not far from 2m χ . This case exhibits behavior similar to the Z 3 -symmetric scenario, but is exceedingly fine tuned, with the relic density being generated by annihilation near the Z pole while the FGCE is obtained through LSP pair annihilation through the a pole, as in the model of [11]. (vi) We note that the fundamental parameters ξ (F,S) have dimensions of mass (2,3) , respectively, so that it is most useful to present their values as a mass in TeV (2,3) . Although a sprinkling of models have values of ξ F,S that extend to both very large positive and negative mass (2,3) , both of these distributions peak in the negative region with values of order a few TeV (2,3) as might be expected. (vii) Lastly, µ ef f is seen to slightly prefer positive values of a few TeV with the mass range 100 GeV being excluded by the lower bound on the Higgsino mass from LEP. For positive values of µ ef f we see that slightly larger values of λ are preferred.
We now investigate how the restricted ranges of the fundamental parameters influence the values of the physical quantities of interest. Phenomenologically, one of the most important aspects of our general NMSSM scenario is the mixing within the CP-odd and CP-even Higgs sectors; this is described by the angle sin θ a as well as the closely related quantity sin θ h , which essentially measures the Higgs singlet content of the ∼ 125 GeV Higgs boson, h. In particular, we are interested in the correlation of sin θ a with the masses of the CP-odd scalars a, A. Fig. 7 and the last panel of Fig. 5 show these results in some detail and there are several important observations to make: (i) The light pseudoscalar mass can range up to ∼ 500 GeV, although the available parameter space becomes increasingly restricted as the pseudoscalar mass increases. We will examine this point in more detail below. (ii) Although we expect that | sin θ a | ∼ 0.1 based on the simpler analysis above, we see that values as small as ∼ 0.001 are possible, as are values as large as ∼ 0.7 in rare cases. Of course the value of the mixing is highly correlated with m a as can be seen in the top panels of Fig. 7. For fixed values of κ, as m a decreases we see that smaller values of | sin θ a | are necessary to obtain the observed relic density. Similarly, the larger the value of |κ| for fixed m a the smaller | sin θ a | is allowed to be. (iii) Furthermore, the region with large tan β and low m a allows for both |κ| and | sin θ a | to be relatively small, whereas larger values of tan β and |κ| permit a large m a . This is not too surprising as the dark matter annihilation amplitude scales roughly as ∼ κ tan β sin θ a /(m 2 a − 2m 2 χ ) away from the a pole. Since models with m a near its upper limit require nearly maximal values of | sin θ a |, |κ|, and tan β, it is not possible to increase the pseudoscalar mass significantly beyond this limit while maintaining coupling perturbativity. (iv) The last panel in Fig. 5 shows that sin θ a is also correlated with the values of λ and A λ ; this is reasonable since in the limit that the LSP mass can be neglected, the off-diagonal  entry in the effective 2 × 2 CP-odd Higgs matrix is λA λ v as can be seen above. (v) Of course, both eigenvalues of the CP-odd Higgs mass matrix are correlated with sin θ a . First, we see in the left panel of Fig. 8 that the region with the highest model density in the m a − m A plane occurs at relatively small m a ∼ 100 GeV, and very large values of m A of order a few TeV. The second panel in Fig. 8 displays the values of | sin θ a | in the m a,A plane. Obviously for fixed m A , as m a increases | sin θ a | must increase accordingly in order to satisfy the relic density measurement. On the other hand, increasing m A with m a fixed lowers the value of | sin θ a | as we would expect from the diagonalization of a 2 × 2 real symmetric matrix. The lower bound on M A from LHC searches increases with m a , since heavy pseudoscalar masses require larger values of tan β (for which the LHC searches are most effective), as we saw in Fig. 7. As a result, the heavy pseudoscalar mass becomes increasingly constrained as m a increases.
The observed value of the branching fraction BR(B s → µ + µ − ) is quite close to the SM prediction and provides a significant constraint on the NMSSM parameter space, particularly for the CP-odd Higgs sector. When this sector supplies the dominant SUSY contribution to this decay the amplitude behaves as ∼ tan 3 β (sin 2 θ a /m 2 a + cos 2 θ a /M 2 A ) and the MSSM limit is recovered when θ a → 0. The first two panels in Fig. 9 present the value of BR(B s → µ + µ − ) as a function of both pseudoscalar masses, m a,A , taken as a ratio to the SM value and colorcoded according to the density of allowed models in each bin. Here we see that the bulk of our models predict essentially the same BF as does the SM, but we do have tails stretching both below and above this value. Note that as m A increases the distribution of predicted branching fractions narrows appreciably and becomes more localized about the SM expectation, while for large m a there is a bifurcation in the set of predictions producing a gap below the SM expectation. The deviation from the SM prediction depends strongly on µ ef f , with positive values of µ ef f enhancing the decay rate and negative values suppressing it compared withe the SM expectation, probably as a result of interference with the SM amplitude. The lower panel in Fig. 9 shows the behavior of BR(B s → µ + µ − ) with respect to the value of the mixing sin θ a . As expected, we see that as m a increases larger values of | sin θ a | are allowed, and that for a fixed value of m a larger values of | sin θ a | generally increase the branching fraction.
Mixing between the doublet and singlet fields within both the CP-odd and CP-even Higgs sectors influences the decay properties of the lighter states a and h. The top two panels in Fig. 10 display histograms of the branching fractions for the decays h, a → bb, τ + τ − and LSP pairs χχ, along with the a → tt, Zh decays when kinematically allowed. In the case of the light CP-even Higgs, we see that the branching fractions for both the bb and τ + τ − final states are generally very close to their SM values since | sin θ h | is generally very small, below ≤ 10 −2 . When this mixing takes on its larger values, the 125 GeV Higgs picks up a significant singlet component which is sufficient to allow for a sizable branching fraction into LSP pairs. However, for most models the BR(h → χχ) usually lies below 1% and is thus likely to be inaccessible, even to measurements made at the ILC [42]. Models with larger branching fractions to the χχ final state have correspondingly reduced branching fractions to bb and τ + τ − , resulting in the tails in the corresponding distributions as observed in the top-  on the value of |κ| and the mixing | sin θ h |. We see that increasing the magnitude of κ, which controls the strength of the singlet-LSP coupling, influences the branching fraction of h into LSP pairs; in particular, for a fixed range of values of | sin θ h |, larger values of |κ| result in a corresponding increase in the h branching fraction to LSP pairs. For the light pseudoscalar a, the reverse situation holds -when | sin θ a | is small, a is mostly singlet and nearly always decays to LSP pairs as seen in the top-right panel of Fig. 10. However, as observed above | sin θ a | can be quite large, 0.1 − 0.5, in which case the a can have significant branching fractions to the bb and τ + τ − final states, although the branching fractions to Zh and tt remain small. For most models, however, even the bb and τ + τ − branching fractions are quite small, with, e.g., BR(a → τ + τ − ) typically below ∼ 1%. Figure 11 shows the correlations among the h and a branching fractions into the bb, τ + τ − and χχ final states; these results exhibit the behavior discussed above. We see that the branching fractions for the bb and τ + τ − modes are generally highly correlated for both h and a, although there are a few exceptions due to the presence of large radiative corrections. The lower two panels of this Figure show the anti-correlation of both the h and a decays to LSP pairs with their decays to SM final states. Here we see, e.g., that the h → bb branching fraction remains in the SM range until the singlet admixture becomes large enough to generate a sizable decay rate into LSP pairs. Note that BR(a → χχ) is nearly unity for most of the models in this set and thus the decay rates to SM final states are generally very small. Now that we understand the decays of the light pseudoscalar a, we can examine the potential constraints on our NMSSM model set arising from searches for additional Higgs fields at the LHC. We note that NMSSMTools implements the well-known constraints from   searches for the MSSM states H and A decaying to τ + τ − , in both the gluon fusion and bb-associated production channels, so that the current limits were applied during the model generation process. In Fig. 12 we examine constraints on bb-associated production of a and their dependence on |κ|. The first panel of the Figure shows the rate for associated production with a decaying to τ + τ − . We see the rate decreases dramatically with increasing |κ|. This rate decrease results from the corresponding increase in the three-point coupling aχχ, which has three important effects. First of all, a large aχχ coupling reduces the a branching fraction to τ + τ − for a given aτ + τ − coupling. Second, models with a large aχχ coupling need a smaller value of the abb coupling to obtain the correct annihilation cross section, resulting in a diminished bb-associated production cross section. Third, since both the abb and aτ + τ − couplings are proportional to tan β × sin θ a , reducing the abb coupling also reduces the aτ + τ − coupling, further suppressing the branching fraction to τ + τ − . We see that the LHC search provides an important constraint on models with small values of |κ|, especially at larger values of m a . We do not show the constraints from production of a via gluon fusion since they are generally weaker than the constraints from bb-associated production.
Of course, if a decays mostly into LSP pairs then bba production may be observable as bb+MET; this channel has been recently discussed in Ref. [43]. In the limit where the aχχ coupling is much larger than the abb coupling (so that BR(a → χχ) 1), which is valid for most models in our sample, the signal strength in this channel is completely described by the a mass and the abb coupling. With this assumption, we can compare the abb coupling in our models to the limit in the lower panel of Figure 5 in [43] with appropriate rescaling. We show the result, including its |κ| dependence, in the right panel of Fig. 12. Once again we see that models with large values of |κ| are more difficult to observe, however the dependence is much more mild since the branching fraction for a to LSP pairs is typically ∼ 1 and therefore essentially independent of |κ|, meaning that the only impact of |κ| is through its effect on the abb coupling. We see that our models are safely below the bound, although searches specifically targeted at this signal could conceivably be competitive with the ditau searches in the future.
We now study the properties of the dark matter sector for our NMSSM model points. The top-left panel in Fig. 13 shows the distribution for the present-day total LSP pair annihilation cross section for our NMSSM model set. Recall that the lower limit of this histogram was set by the requirement that we reproduce the observed FGCE based on the analysis of Ref. [4]. Note that the annihilation rate for the majority of the models lies above < σv >= 1.5 × 10 −26 cm 3 s −1 . The top-right panel presents a histogram of the annihilation fraction to the bb final state. By design, the most common final state for the LSP annihilation process is indeed bb, which typically comprises ∼ 90% of the total rate. However, as shown in the bottom-left panel of Fig. 13, other final states, such as τ + τ − and gg, can also have significant rates, with the τ + τ − channel accounting for ∼ 10% of the annihilation rate in many models. The bottom-right panel shows the correlation between the present-day dark matter annihilation cross section and that at freeze-out, which accounts for the observed relic density (recall that we employed the NMSSMTools version of micrOMEGAs). We see that the The abb coupling as a function of m a . A limit on the abb coupling, assuming a ∼ 100% BF to LSP pairs, was obtained by rescaling the limit in Figure 5 of [43], and is shown as the solid black line. Both panels are color-coded by the value of |κ| bulk of the models reside in the diamond-shaped region in the upper left corner of the panel, where the two annihilation cross sections have comparable values with the present-day rate being slightly smaller as expected. However, we also find a long tail of model points with larger freeze-out cross sections and smaller present-day annihilation rates. In these cases, we find that m a is close enough to 2m χ for the annihilation rate to benefit significantly from the a pole during freeze-out, resulting in a suppression of the present-day annihilation rate which is not resonantly enhanced. Some of these models also have a slightly elevated Higgsino content of the LSP, so that Z exchange can also play a small role during the annihilation process in the early universe.
In Fig. 14 we further examine the correlation between the present-day and freeze-out dark matter annihilation cross sections. Here, the models are colored-coded by the values of the LSP mass relative to both the a and the Z, so that we can ascertain the effects of these specific resonance exchanges. As expected, we see in the left panel that the tail of the distribution that descends to lower present-day annihilation rates and larger rates during freeze-out occurs as 2m χ approaches m a . This indicates that the a pole is indeed enhancing the early universe annihilation rate for these models as we expected. Interestingly, at the bottom of this long tail we find a set of models (colored blue and green) where the LSP is quite light relative to the a, but where the present-day cross section is still suppressed. The right panel of Fig. 14 displays the same two annihilation cross sections with the points now color-coded by the ratio of the LSP and Z masses. Here, the range of mass ratios is a priori narrow as the LSP was required to lie in the 30-40 GeV mass region, however some variation is noticeable with models having lighter LSPs tending to cluster on the top-left side of the distribution. In this panel we see that most of the models in the tail have masses not too far Histogram of the present-day annihilation rate in the bb final states, taken as a ratio to the inclusive annihilation rate. Bottom-left panel: Histograms of the present-day annihilation rate into bb, τ + τ − , and gg final states, taken as a ratio to the inclusive annihilation rate. Bottom-right panel: Correlation between the inclusive present-day annihilation rate and its value at freeze-out. from ∼ 0.4m Z so that the influence of the Z pole is likely non-negligible, accounting for the models which are far from the a pole but still have suppressed present-day cross sections. Figure 14: Correlation between the inclusive present-day annihilation rate and its value during freeze-out, color-coded by the ratio of the LSP mass to the mass of the a (left) and Z (right) bosons.
Lastly, we consider both the spin-independent (SI) and spin-dependent (SD) direct detection cross sections for our model points in Fig. 15, along with the current LUX and COUPP constraints. The SI cross section is a measure of the interaction between the LSP and the nucleus, and in this scenario is mostly due to CP-even Higgs exchange resulting from mixing via, e.g., a non-zero value for sin θ h or a small Higgsino admixture in the LSP, as well as through a possible loop-induced hχχ coupling. Note that the t-channel exchange contribution is small since the first and second generation squarks are quite heavy in the cases we consider here. The SD interaction, on the other hand, can arise from Z exchange via the Higgsino content of the LSP and/or through the exchange of the CP-odd Higgs fields, although the latter process is suppressed by the values of both the mixing angles well as v 2 . We thus expect that both the SI and SD direct detection cross sections for the singlino-like LSPs in our model set will be quite small. Indeed, this is what we observe in the top-left panel of Fig. 15, which shows the model density distribution in the SI-SD cross section plane.
Here we see that typical values for the SI (SD) cross sections are of order ∼ 10 −13 (10 −10 ) pb, both of which are several orders of magnitude away from the current experimental bounds. Unfortunately only a small fraction of our model points will be accessible to future direct detection searches. The top-right panel in this Figure shows how the magnitude of the SD cross section is strongly correlated with the Higgsino content of the LSP, denoted here as Z 2 13 + Z 2 14 and represented by the point coloration. Models where the LSP has a larger Higgsino content have a larger coupling to the Z boson, resulting in sizable SD scattering through Z exchange. On the other hand, we find that the SI cross section is apparently not sensitive to the Higgsino content of the LSP for this set of models. The bottom-left panel shows the dependence of the SI and SD cross sections on the value of | sin θ a |, which could be significant if a, A exchange are the dominant contribution. Here we see that, at most, there is a modest sensitivity implying that the SD cross section is indeed dominated by Z exchange. The bottom-right panel shows the dependence of the SI cross section on the value of | sin θ h |. The strong dependence of the SI cross section on | sin θ h | shows that the SI scattering proceeds mainly through h exchange, via the mixing between h and the scalar singlet Higgs. However, since there are, in fact, three CP-even Higgs bosons of various masses and couplings contributing simultaneously to the SI cross section, certain values of their masses and mixings can result in interference between the various contributions. The effect of this interference can be seen as the orange model points in the top left side of the Figure, for which relatively large values of | sin θ h | still result in rather small SI cross sections.

Conclusions
In this paper we have considered a set of SUSY scenarios, searching for a model that might lead to a consistent explanation of the FGCE while also satisfying other known experimental and theoretical constraints. We began by exploring a variety of SUSY models which failed to describe this excess, all for different reasons. These attempts demonstrate the significant model building challenge posed by the FGCE within the SUSY context.
In the (p)MSSM a ∼ 35 GeV neutralino can lead to the correct relic density only if it is a bino-Higgsino admixture that annihilates near the Z pole in the early universe. 5 At the present time, dark matter annihilation takes place far from the Z pole and so the prediction for the FGCE flux is too small by several orders of magnitude. If the LSP were instead a Dirac gaugino, this would help to increase the annihilation cross section by avoiding the velocity/helicity suppression that occurs in the Majorana scenario. Of course, in this case the LSP would need to be a very pure bino to avoid the direct detection constraints and so the observed relic density could only be obtained via t-channel sfermion exchange, in particular, through a right-handedτ with mass near the LEP limit. Given the hypercharge of the right-handed sbottom, as well as the lower bound on the sbottom mass, the desired bb final state would not occur at high enough rates to account for either the observed relic density or the FGCE. In the case of an extended gauge symmetry, such as the GUT theory E 6 , new fermions and new gauge bosons are present so that the Z-induced direct detection constraint can be avoided if the LSP were a Dirac (or even Majorana) SM singlet annihilating through a heavy new gauge boson Z . Unfortunately, the LHC lower bounds on a Z mass lead to a highly suppressed DM annihilation cross section for either of the SM singlet DM candidates present in E 6 .
We next turned our attention to the NMSSM, identifying the LSP as an almost pure singlino state. We first investigated the Z 3 version of the NMSSM (with only four new parameters beyond the MSSM) under the assumption that the DM annihilation process that produces the observed relic density (s-channel light CP-odd Higgs exchange) is the same as Figure 15: Scatter plots of the spin-dependent vs spin-independent scattering cross sections for our NMSSM model points, color-coded by the density of allowed models in each bin (top left), the Higgsino content of the LSP, (|Z 13 | 2 + |Z 14 | 2 ) (top right), the pseudoscalar singlet-doublet mixing | sin θ a | (bottom left), and the scalar singlet-doublet mixing | sin θ h | (bottom right). Current limits from COUPP [44] and LUX [39] are shown as horizontal and vertical lines, respectively. that which leads to the presently observed FGCE. Regrettably, we showed that this scenario also proves to be inadequate to the task. This is due to the large Higgsino admixture which is necessarily picked up by the LSP due to the structure of the neutralino mixing matrix, and the insufficient parameter freedom which prevents us from avoiding this mixing in this scenario. Recall that once the LSP has picked up a significant Higgsino component then the Z will play an important role in the dark matter annihilation process and one has to fine-tune the parameters to generate the observed relic density as well as the present-day annihilation in the galactic center. We wished to avoid such fine-tuning, motivating our assumption that CP-odd Higgs exchange is solely responsible for both. Lastly, we examined the general NMSSM which has nine additional superpotential plus soft-breaking parameters beyond those of the usual (p)MSSM and thus has far greater flexibility. To better understand the various phenomenological requirements that the FGCE imposes on this version of the NMSSM, we began with an examination of this model at tree level. Even in this rather simplified case we found that it was non-trivial to meet both the dark matter requirements as well as other experimental constraints. However, viable solutions were shown to exist, selecting a specific region of the parameter space. In order to perform a more complete and detailed analysis, we made use of a modified version of the NMSSMTools4.3.0 code which includes radiative corrections to the various mass spectra and corresponding couplings, and also provided a means to calculate the relic density as well as the dark matter direct detection and present-day pair-annihilation cross sections in this scenario.
We then generated a set of 6 · 10 8 points in the general NMSSM parameter space with rather loose range requirements, subjecting each of them to a set of global experimental constraints (collider, flavor, astrophysical, etc.) and also demanding a sufficiently large present-day dark matter annihilation cross section into bb in order to explain the FGCE. Of this set, ∼ 52.8k model points survived these basic requirements, demonstrating that the general NMSSM with a single dark matter annihilation mechanism provides a successful SUSY scenario for the FGCE. We then studied the physical properties of these surviving model points in some detail. First, we examined the regions of the NMSSM parameter space that led to viable models, studying in particular the role played by correlations among these parameters and the properties of the physical Higgs states. One of the interesting results of this investigation is the rather wide range of possible values for the light CP-odd Higgs mass, m a , (which is correlated with the values of κ and sin θ a , and also to some extent tan β) allowed by the constraints. Since the a has a large singlet component and the LSP is dominantly singlino, a principally decays to LSP pairs. This makes it difficult to observe at the LHC, where searches for b-jets and missing energy face large backgrounds. Although the best constraints come from di-tau resonance searches, future searches for invisible a decays could also provide some sensitivity. In contrast, the light CP-even, ∼ 125 GeV Higgs generally has a rather weak, mixing-induced coupling to the LSP, so that the branching fraction for h → χχ typically lies far below 1%; hence h has properties quite close to that of the SM Higgs, as is the case in the MSSM.
Lastly, we studied the dark matter and flavor properties in our NMSSM model set. Since the LSP only couples to the SM via mixing with the Higgsinos, and/or mixing amongst the Higgs fields, both the SI and SD direct detection cross sections are found to be quite small, i.e., several orders of magnitude away from the current experimental bounds, in this scenario. We find that only a small fraction of our NMSSM models would be accessible to future direct detection searches. Rare B decays, such as B s → µ + µ − , can probe this parameter space to some extent. Although most of the models yield branching fractions quite close to the SM prediction, there are a significant number of cases where the predictions differ from the SM by ∼ 50%. Once the theoretical and experimental uncertainties are under control, B s → µ + µ − could provide a signal for this subset of models.
In summary, we find that it is non-trivial to build a UV-complete model based on SUSY that describes the FGCE. However, we find that the general NMSSM easily accommodates the data and that the data selects a specific region of parameter space. Unfortunately, the signatures elsewhere, such as in Higgs physics at the LHC and in direct detection, are small and will be difficult, if not impossible, to detect. It is possible that a signature could be observed in the flavor sector.
Indirect detection thus provides the best probe of this scenario. Given that the FGCE dictates the dark matter annihilation process cross section, Fermi should soon see some gamma ray excesses elsewhere, e.g., in the dark matter searches using Dwarf Galaxies. We hope that further data will confirm the signature observed at the galactic center and that dark matter will soon be discovered!