Interpreting the galactic center gamma-ray excess in the NMSSM

In the Next-to-Minimal Supersymmetric Standard Model (NMSSM), all singlet-dominated particles including one neutralino, one CP-odd Higgs boson and one CP-even Higgs boson can be simultaneously lighter than about 100 GeV. Consequently, dark matter (DM) in the NMSSM can annihilate into multiple final states to explain the galactic center gamma-ray excess (GCE). In this work we take into account the foreground and background uncertainties for the GCE and investigate these explanations. We carry out a sophisticated scan over the NMSSM parameter space by considering various experimental constraints such as the Higgs data, $B$-physics observables, DM relic desnity, LUX experiment and the dSphs constraints. Then for each surviving parameter point we perform a fit to the GCE spectrum by using the correlation matrix that incorporates both the statistical and systematic uncertainties of the measured excess. After examining the properties of the obtained GCE solutions, we conclude that the GCE can be well explained by the pure annihilations $\tilde{\chi}_1^0 \tilde{\chi}_1^0 \to b \bar{b} $ and $\tilde{\chi}_1^0 \tilde{\chi}_1^0 \to A_1 H_i $ with $A_1$ being the lighter singlet-dominated CP-odd Higgs boson and $H_i$ denoting the singlet-dominated CP-even Higgs boson or SM-like Higgs boson, and it can also be explained by the mixed annihilation $\tilde{\chi}_1^0 \tilde{\chi}_1^0 \to W^+ W^-, A_1 H_1$. Among these annihilation channels, $\tilde{\chi}_1^0 \tilde{\chi}_1^0 \to A_1 H_i $ can provide the best interpretation with the corresponding $p$-value reaching 0.55. We also discuss to what extent the future DM direct detection experiments can explore the GCE solutions and conclude that the XENON-1T experiment is very promising in testing nearly all the solutions.


Introduction
The compelling evidences for the existence of Dark Matter (DM) from various cosmological and astrophysical observations have provided us a good portal in the search for new physics beyond the Standard Model (SM). One possible method to explore DM in the present Universe is the indirect detection, which looks for the particles produced when DM annihilates in the DM halo. These particles include photons, antiparticles and neutrinos, and among them gamma rays have often been defined as the golden channel for DM indirect detection since the signal can be traced back to the source. The Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope, due to its unprecedented angular and energy resolutions, has produced the most detailed maps of the gamma ray sky for a wide range of energies. Intriguingly, as was reported by several independent groups [1][2][3][4][5][6][7][8][9] and also by Fermi Collaboration itself [10], the Fermi-LAT data have revealed the presence of an extended excess of gamma rays over the modeled foreground and background emissions towards the Galactic Center (GC). Although several astrophysical mechanisms, such as the thousands of unresolved millisecond pulsars [11][12][13] and the interactions between comic rays (CR) and interstellar gases [14][15][16][17][18], have been proposed to interpret this Galactic Center Excess (GCE), they usually fail to generate the morphology and energy spectrum of the GCE simultaneously 1 . So in this work, we instead consider another possibility that the GCE is produced by the annihilation of DM. Although this interpretation has been constrained by the measurements of CR such as the Fermi-LAT detection of the gamma-rays from dwarf spheroidal galaxies (dSphs) [20][21][22][23], the non-observation of spectral features in the AMS-02 measurements of CR positron [24][25][26][27], and PAMELA observations of the CR anti-protons [28][29][30][31][32][33][34], it still remains a most attractive one not only because the excess emission shows spectral and morphological properties consistent with a telltale sign from DM annihilation, but also because in such an interpretation, the annihilation cross section required to explain the GCE is of the right size to account for the DM density from thermal freeze-out.
So far there have been a large number of attempts to explain the GCE by DM annihilation in various new physics models . In the early analyses of the annihilations, great efforts were focused on the channelsχχ → bb with mχ ∼ 35GeV andχχ → ττ with mχ ∼ 10GeV since they can reproduce well the GCE spectrum obtained at that time. Recently a critical reassessment of the DM interpretation was made by examining in a comprehensive way the foreground and background uncertainties [9]. It was found that taking the estimated uncertainty in the high-energy tail of the spectrum into account, a much larger number of DM annihilations are able to fit well the γ-ray data than previously noted [73,79]. Explicitly speaking, as far as the annihilationχχ → bb is concerned, now the mass of DM is extended to a broader range from 30GeV to 70GeV in explaining the GCE [73,79]. Other annihilation channels such as DM annihilation into light quark pairs and even gluon pair are also able to provide a good fit to the GCE [79]. More strikingly, this new analysis opens up a very good solution usually neglected before, namely DM annihilation into a pair of light non-standard Higgs bosons [49,50]. This important progress motivates us to renew the solutions to the GCE in supersymmetric theories, which usually predict the lightest neutralinoχ 0 1 as a natural DM candidate. As the most economical realization of supersymmetry, the Minimal Supersymmetric Standard Model (MSSM) is unsatisfactory in explaining the GCE due to the following four reasons [66,93]. First, the relic density of DM has required its mass to be larger than about 40GeV [103]. In this case, the annihilationsχ 0 1χ 0 1 → ττ , qq with q denoting a light quark can not provide a good fit any more. Second, except for excessive fine-tuning cases the LHC experiments have pushed the lower mass bounds for the CP-odd Higgs boson and the bottom squarks up to several hundred GeV. As a result, the cross section of DM annihilation into bb in present day is too small to significantly contribute to the GCE [66,93]. Third, due to the small velocity of DM in our galaxy, the annihilation rate for DM into SM-like Higgs pairs is p-wave suppressed. Consequently this channel is not large enough to generate the GCE. Finally, as for the annihilationsχ 0 1χ 0 1 → W W, ZZ, their fits to the GCE spectrum indicate that regardless of their annihilation rates the corresponding p-values are always less than 0.04 [73,79]. This means that the annihilations can not generate the proper spectrum shape for the GCE. We note that for a given parameter point of the MSSM, DM usually annihilates into multiple final states. In this case, the situation can not be improved greatly because, due to the particle spectrum of the MSSM allowed by the current experiments, either the total cross section falls short for the GCE, or the dominant annihilation channel can not reproduce the GCE spectrum well [92].
Given the problems of the MSSM, we consider to interpret the GCE in the Next-to-Minimal Supersymmetric Standard Model (NMSSM) with a Z 3 symmetry, which is the simplest gauge singlet Higgs extension of the MSSM [104]. Distinguished from the MSSM, the NMSSM predicts three singlet-dominated particles: one neutralino, one CP-even and one CP-odd Higgs bosons. These particles are rather special in that all of them can be simultaneously lighter than about 100 GeV, and that the couplings for the interactions among themselves are determined by the parameter κ, which alone is able to predict the right rates for some annihilation channels to explain the GCE (see the following discussion). These features make the NMSSM with a singlet-dominated DM well suit to account for the GCE because, as we will show below, some golden channels for the GCE need light particles to act as the DM, the mediator and/or the annihilation final state.
We note that the interpretations of the GCE in the NMSSM have been intensively discussed in [57,58,66,68,69,93]. However, in [57,58,66,68,69] the authors did not consider the systematic uncertainties mentioned above. As a result, the model parameter space they considered is much narrower than that of this work and the obtained conclusions were incomplete. While for [93], although the authors have taken the uncertainties into account, they considered the parameter space characterized by a large λ which is different from our discussion.
The aim of this work is to explore any possible solution to the GCE in the Z 3 NMSSM. For this end, we perform a sophisticated scan over the model parameters by considering various experimental constraints such as the DM relic density, the Higgs data as well as the observation of dwarf galaxies. We use the correlation matrix presented in [9] to include the systematic uncertainties on the GCE spectrum and only keep the parameter points that can reproduce well the spectrum. In our study we mainly consider a singlino-like DM which is believed to interpret the GCE without excessive fine tuning. As we will show below, the annihilationχ 0 1χ 0 1 → H i A 1 with H i A 1 denoting a scalar-pseudoscalar Higgs pair may provide the best fit to GCE, and the canonical annihilationχ 0 1χ 0 1 → bb still remains a satisfactory solution except that mχ0 1 is now allowed to vary within a broader range. Moreover, it is interesting to see that the mixed annihilation into W + W − and H i A 1 final states is also able to generate a spectrum consistent with the GCE. These conclusions are quite different from previous studies in the NMSSM. This paper is organized as follows. In Section II, we introduce some of the characteristic features of NMSSM, the basic knowledge about the GCE and our strategy for the parameter scan. In Section III, we discuss in detail the interpretations of the GCE when H 2 is the SM-like Higgs boson, and in Section IV, we carry out a similar study but for the case that H 1 acts as the SM-like Higgs boson. We draw our conclusion in Section V and provide more information of the NMSSM couplings in the Appendix.

Theoretical setup for the GCE in the NMSSM
We start our analysis by recapitulating the basics of the NMSSM. As one of the most economical extensions of the MSSM, the NMSSM introduces one gauge singlet Higgs superfield in its matter content, and since one purpose of the extension is to solve the µ-problem of the MSSM, a Z 3 symmetry is usually adopted in the construction of the superpotential to avoid the appearance of parameters with mass dimension. As a result, the superpotential of the NMSSM and the soft breaking terms in Higgs sector are given by [104] where W F is the superpotential of the MSSM without the µ-term,Ĥ u ,Ĥ d andŜ are Higgs superfields with H u , H d and S acting as their scalar components respectively, the dimensionless coefficients λ and κ parameterize the strengthes of the Higgs self couplings, and m u ,m d ,m S , A λ and A κ are soft-breaking parameters. In practice, after the electroweak symmetry breaking the soft-breaking squared massesm 2 u ,m 2 d andm 2 s are traded for m Z , tan β ≡ v u /v d and µ ≡ λv s as theoretical inputs.
Due to the presence of the superfieldŜ, the NMSSM contains a singlino field which is the fermion component ofŜ, and one more complex Higgs field S compared to the MSSM. As a result, the neutralino mass eigenstatesχ 0 i (with i ranging from 1 to 5) are the mixtures of bino, wino, higgsinos and singlino, and the CP-even (odd) Higgs mass eigenstates H i with i = 1, 2, 3 (A i with i = 1, 2) are mixtures of the real (imaginary) parts of H u , H d and S. Throughout this paper, we assume the mass order mχ0 1 < mχ0 2 < · · · < mχ0 5 for neutralinos, and m H 1 < m H 2 < m H 3 , m A 1 < m A 2 for Higgs bosons.
There are three distinguished features in the NMSSM. One is that DM in the NMSSM may be either singlino-dominated or bino-dominated. As expected, the properties of a singlino-dominated DM are quite different from those of a bino-dominated DM, which makes the DM physics in the NMSSM much richer than that in the MSSM [105]. Another feature is that, in the presence of a singlino-dominated DM with mass below 100GeV, the singlet-dominated CP-even and CP-odd Higgs bosons can be simultaneously lighter than about 100GeV [105,106], and the strengthes for the interactions among these particles are determined by the parameter κ which may be as large as 0.1. This feature, as we will show below, makes the NMSSM with a singlino-dominated DM well suit to explain the GCE. In the appendix, we list the properties of these particles used in our analysis. The other feature is that either H 1 or H 2 in the NMSSM can act as the SM-like Higgs boson [107] and generally speaking, H 2 as the SM-like boson is more attractive from phenomenological point of view and also from naturalness argument.
In the DM explanation of the GCE, the observed γ-ray originates mainly from the cascade decays of the annihilation final states. In the NMSSM, the possible annihilation final states include ff , V V , H i H j , A i A j and H i A j [108], where f (V ) denotes any of the fermions (vector bosons) in the SM, and H i (A j ) denotes a CP-even (CP-odd) Higgs boson. In this work, we are particularly interested in the annihilationsχ 0 1χ 0 These annihilations proceed through the s-channel mediator of a Z boson or a Higgs boson with an appropriate CP quantum number, and also proceed through the t/u-channel exchange of a sbottom, a chargino and a neutralino respectively. The complete expressions of the annihilation cross sections are rather complicated, but in non-realistic limit, i.e. the velocity of DM approaching zero, some contributions become unimportant. In this case, the velocity weighted annihilation cross section can be approximated by [108] where C XY Z denotes the coupling of the interaction involving the particles X, Y and Z, 1 . In getting Eq.(2.5), we note that a good fit to the GCE requires that the H i A 1 final state is produced close to threshold, i.e. δ 0, so we can expand σ H i A 1 v 0 in terms of δ. Then the first two terms on the right hand of Eq.(2.5) come from the left diagram of Fig.1, and the last term comes from the right diagram of Fig.1.
The flux per unit solid angle at some photon energy E γ , which is observed by Fermi- LAT, is then given by where dN γ XY /dE is the photon spectrum generated by the annihilationχ 0 1χ 0 1 → XY , ρ DM is the DM profile and the integral over ρ 2 DM is along the light-of-sight (LOS) at an angle θ towards GC. In the DM interpretation of the GCE, a generalized Navarro, Frenk & White (NFW) DM profile is usually adopted, and its expression is given by [109,110] ρ(r) = ρ r r with slope parameter γ = 1.26, scale radius R s = 20 kpc and the local DM density ρ = 0.4 GeV/cm 3 at the radial distance of the sun from the galactic centre r . Here the coordinate r is centered on the galactic centre and can be expressed as r 2 (s, θ) = r 2 + s 2 − 2r s cos θ with s and θ being the LOS distance and the aperture angle between the axis connecting the earth with the galactic centre and the LOS respectively. In our study, we use the package micrOMEGAs-3.6.9.2 [111] to calculate the DM relic density and with the help of PYTHIA [112] to generate the flux in Eq.(2.6). Note that in any explicit model, DM usually annihilates into multiple final states. In this case, the different fluxes are summed over.

Parameter scan strategy for GCE solution
We simplify our scan over the NMSSM parameter space by fixing the parameters that are not closely related to the DM studies. The soft SUSY breaking parameters in the squark sector are all fixed to be 2 TeV except that we vary those for the third generation to generate a CP-even Higgs near 125 GeV. We assume A t = A b and M U 3 = M D 3 to reduce the number of free parameters. Similarly, all of the soft SUSY breaking parameters in the slepton sector are fixed to be 300 GeV to explain the discrepancy of the measured value for muon anomalous magnetic moment from its SM prediction. As for the gaugino sector we abandon the Grand Unified Theory relation and fix the wino mass and gluino mass at 1 TeV and 2 TeV respectively. Consequently, the remained free parameters include tan β, µ, λ, κ, A λ , A κ in the Higgs sector, M Q 3 , M U 3 and A t for third generation quarks and the bino mass M 1 , which are all defined at the scale of 2 TeV in the scan. We use NMSSMTools-4.3.0 [113] to scan intensively the following NMSSM parameter region: (2.8) Table 1. Favored parameter region of the NMSSM to explain the GCE, which are classified by the dominant final state in DM annihilations. These annihilations are called Solution I, II, III, IV and V respectively in the following discussion. All input parameters are defined at 2TeV and quantities with mass (annihilation cross section) dimension are in unit of GeV (10 −26 cm 3 /s). The process to retain the parameter points include the following steps: • We require the DM to be singlino-dominated and satisfy mχ0 1 ≤ 150 GeV, and impose all the experimental constraints encoded in NMSSMTools-4.3.0 [113] which include the relic abundance at 3σ level (0.107 ≤ Ωh 2 ≤ 0.131), LUX exclusion bound at 90% C.L., various B-physics measurements as well as the discrepancy of muon magnetic moment at 2σ level. We also consider various electroweak precision data calculated in [114].
• We consider the constraints on the Higgs sector with the package HiggsBounds-4.1.2 [115] which contains the data from LEP, Tevetron and LHC. For the SM-like Higgs boson, we further perform a fit to the data with the package HiggsSignal [116] and keep the 2σ samples.
• We use micrOMEGAs-3.6.9.2 [111] to calculate the DM annihilation cross section at present day, and then impose the constraints from dSphs by the data in [23] for the bb annihilation channel and with the method introduced in [93] for the H i A 1 final states.
• We also use micrOMEGAs-3.6.9.2 [111] to generate the γ-ray spectrum. Considering the astrophysical uncertainties which may come from the errors in our setting on the local DM density ρ , the scale radius R s and the inner slope parameter γ in Eq. (7), for each parameter point we allow an uncertainty factor A in the range of (0.17, 5.3) for the annihilation cross section, or equivalently for the height of the gamma-ray spectrum in Eq.(2.6) [79]. Then for the A−tuned γ-ray spectrum, we perform a fit to the residual GCE spectrum obtained in [9] by using the publicly available covariance matrix, which include both the statistical and systematic uncertainties of the measured flux. The corresponding χ 2 sp function is calculated by [9,79]: where Σ ij is the covariance matrix, dN/dE i is the measured flux in the i-th energy bin, and dN /dE i is the flux predicted by the NMSSM, which depends on the parameter point and also on the factor A.
We define the GCE χ 2 as the minimum value of χ 2 sp (A) among different choices of A, χ 2 GCE = min(χ 2 sp (A)), and keep the parameter points that satisfy χ 2 GCE ≤ 35.2. These points are assumed to have the capability to explain the GCE at 95% confidence level for 23 degree of freedom [9].
The parameter ranges of the GCE solutions are listed in Table.1, which are classified by the dominant final state in DM annihilations (see the following discussion). For the first three types of the DM annihilations H 2 acts as the SM-like Higgs boson, while for the last two types H 1 corresponds to the SM-like Higgs boson. One distinguished feature that Table  1 exhibits is that all the singlet dominated particles in the GCE solutions, including DM, the singlet-dominated CP-even and CP-odd Higgs bosons, are lighter than about 150GeV. This feature, as we will emphasized below, makes the NMSSM well suit for explaining the GCE.

GCE solutions with H 2 being the SM-like Higgs boson
In this section, we exhibit the features of the GCE solutions for the case that DM is singlino-dominated and H 2 acts as the SM-like Higgs boson. All the solutions considered in this work survive the constraints listed in last section and meanwhile can explain the GCE at 95% C.L..
In Fig.1 we project the solutions on σv 0 − mχ0 1 plane (upper panel) and χ 2 GCE − mχ0 1 plane (lower panel). Solutions marked by red square, blue triangle and black asterisk correspond to the cases that DM annihilates with the largest branching ratio into bb, H 1 A 1 Table 2. Detailed information of the benchmark points used in our discussion. Quantities with mass, annihilation and scattering cross section dimension are in unit of GeV, cm 3 /s and pb respectively.   and W + W − final states respectively, which hereafter are collectively called Solution I, Solution II and Solution III correspondingly. Then the upper panel of Figure.2 indicates that, for the ranges 30GeV ≤ mχ0 1 ≤ 40GeV, 50GeV ≤ mχ0 1 ≤ 62GeV and 63GeV ≤ mχ0 1 ≤ 70GeV, Solution I is viable, while for 63GeV ≤ mχ0 1 ≤ 115GeV and 83GeV ≤ mχ0 1 ≤ 100GeV, Solution II and Solution III can account for the GCE respectively. For any of the solutions, the σv 0 is larger than 1.7 × 10 −27 cm 3 /s, and its lower bound increases monotonically asχ 0 1 becomes heavier. The reason for the latter behavior is that, for a heavier DM, its number density is smaller. So to obtain the same photon flux for the GCE, a larger cross section is needed.
The lower panel of Fig.2 indicates that the best interpretation in Solution I comes from mχ0 1 50GeV with χ 2 GCE 23 and a p-value of 0.44. This conclusion coincides with that of [79], which was obtained in a model independent way and for a pure bb annihilation channel. For Solutions II and Solutions III, the best interpretations locate at  that, in the case of mχ0 1 m Z /2 (mχ0 1 62GeV), DM annihilated in early universe mainly through a nearly on-shell Z boson (SM-like Higgs boson). Since nowadays this dominant annihilation channel is helicity (p-wave) suppressed, σv 0 can not reach the size required for the GCE.
In Table 2, we present detailed information of three benchmark points P1, P2 and P3 for Solution I, II and III respectively. This table indicates that the sole annihilation channelχ 0 1χ 0 1 → bb orχ 0 1χ 0 1 → H 1 A 1 can be responsible for the GCE; while for the channel χ 0 1χ 0 1 → W + W − , it must mix sizeably with the channelχ 0 1χ 0 1 → H 1 A 1 to account for the GCE. We will return to this issue later.
In our calculation, we found that the condition on the GCE χ 2 can reduce the number of the parameter points that survive the constraints by more than 90%. This implies that the GCE has non-trivial requirements on the parameters of the NMSSM, especially it suggests that some of the independent parameters may be correlated. Motivated by this thought, we study the correlations among the parameters λ, κ, µ, mχ0 1 and m A 1 which are important parameters in the interpretation of the GCE and show the corresponding results in Fig.3. In the following, we concentrate separately on each kind of the solutions and investigate its features. Such a study is helpful to understand the correlations in Fig.3 and also the properties of the benchmark points listed in Table.2.

Solution I -the bb annihilation channel
Among the solutions to the GCE, Solution I is the most intensively studied one. After considering the systematic uncertainties, one important improvement of Solution I over its previous version is that DM mass is now allowed in the range from 30GeV to 70GeV, which is much wider than before.
The key features of Solution I are as follows: • The lighter CP-odd Higgs boson is correlated with DM by m A 1 2mχ0 1 . This correlation is shown in the upper right panel of Fig.3 which means that the annihilation proceeds resonantly.
This feature can be understood as follows. In Solution I, the heavy CP-odd Higgs boson is doublet-dominated with its mass usually at TeV scale. Then Eq.(2.3) indicates that the main contribution to the annihilation comes from the moderately light A 1 , which is singlet-dominated. With the formula presented in Eq.(A.20) and v s ≡ µ/λ 450GeV shown in the lower left panel of Fig.3, one can get This inequation means that the couplings involved in the annihilation are highly suppressed so that the process must proceed resonantly to ensure σ bb v 0 ∼ 10 −26 cm 3 /s. Moreover, our results indicate that the width of A 1 is very small, Γ A 1 10 −2 MeV. So as m A 1 approaches 2mχ0 1 , the denominator in Eq.(2.3) tends to vanish and a small κ in Eq.(3.1) is then suffice to predict the right rate of the annihilation for the GCE. This character is illustrated in the upper left panel of Fig.3. In fact, a small κ is also favored to predict lightχ 0 1 and A 1 , which can be seen from Eq.(A.3) and Eq.(A.8).
• The parameter µ is upper bounded by about 300GeV, which is shown in the lower panels of Fig.3.
This feature is actually required by the DM relic density [69]. Generally speaking, in order to predict the measured Ωh 2 , the velocity weighted cross section σv should be around the canonical value 3 × 10 −26 cm 3 /s at freezing out (see for example points in Table 2). Since 2m χ 0 1 /m A 1 > 1 in Solution I, σv for the annihilatioñ χ 0 1χ 0 1 → A * 1 → bb at present day is usually larger than that at freezing out due to the thermal broadening [117]. Since the dwarf galaxy measurements have required σ bb v 0 2 × 10 −26 cm 3 /s (see Fig.2), new contributions such as those mediated by a Z boson or a CP-even Higgs boson must intervene for the DM annihilation in early Universe, and a moderately small µ can accelerate the annihilation [69].
• Sinceχ 0 1 60GeV for most cases in Solution I, the SM-like Higgs boson H 2 may decay intoχ 0 1 pair. We checked that Br(H 2 →χ 0 1χ 0 1 ) 18%, which is required by the Higgs data at the LHC.

Solution II -the H 1 A 1 annihilation channel
Solution II is quite similar to the interpretations presented in [49,50,83,93,99,100], which utilize the processχχ → φ 1 φ 2 → f 1f1 f 2f2 (φ 1 and φ 2 denote scalar or pseudoscalar particles, and f 1 and f 2 are SM fermions) for the GCE. These interpretations, as were emphasized by the proposers, can easily escape the constraints from DM detection experiments and have been paid more and more attention recently.
The features of Solution II are as follows: • The singlet-dominated particles satisfy 60GeV mχ0 1 115GeV, 10GeV m A 1 110GeV, 60GeV m H 1 120GeV and δ < 0.2, and for most samples there exist following relations m H 1 mχ0 1 m A 1 . Given κ ∼ 0.1 which is required to predict the right size of the annihilationχ 0 1χ 0 1 → H 1 A 1 for the GCE (see below), the particle spectrum limits parameters such as λ, µ and A κ in certain regions (see the expressions of the tree level masses in Appendix), which are given in Table 1, and also shown in Fig.3.
Note that µ is below about 200GeV. In this case, the higgsino-dominated neutralinos χ 0 i may decay dominantly intoχ 0 1 A 1 instead of intoχ 0 1 Z since the kinematics is forbidden. In this case, the LHC search for electroweakinos by trilepton +E miss T signal is less efficient in ruling out the light higgsinos 2 . Also note that the parameters λ and µ are related by µ/GeV ≈ 60 + 260 λ for λ varying from 0.2 to 0.6 (see lower left panel of Fig.2), which means that v s ≡ µ/λ > 360GeV. This ensures that the expansions for the masses and couplings in Appendix by the power of λv/µ are good approximations.
• The s-channel contributions to the annihilation rate σ H 1 A 1 v 0 in Eq.(2.5) are usually much smaller than those from the t/u channel, and among the t/u channel contributions, the one induced by the exchange ofχ 0 1 is far dominant. As for the contributions induced by the two higgsino-like neutralinos, each of them may be sizable, but since they cancel each other, the net higgsino contribution is not important. These characters can be understood by the following approximations (see Eq. and by the fact that κ ∼ 0.1 is enough to predict theχ 0 1 contributed σ H i A 1 v 0 at the order of 10 −26 cm 3 /s (see equation (3.20) in [93]).
• Since m A 1 60GeV for most cases in Solution II (see upper right panel of Fig.2), the SM-like Higgs boson H 2 may decay into A 1 A 1 with a sizeable fraction. Given that A 1 decays dominantly into bb, this will result in 4b signal for the SM-like Higgs boson. We checked that Br(H 2 → A 1 A 1 ) 24%, where the upper bound comes from the constraints of the LHC Higgs data.
• Since a good fit to the GCE requires that H 1 A 1 is produced close to threshold, the annihilationχ 0 1χ 0 1 → H 1 A 1 will produce spectral line or box-shaped spectrum in γ-ray [83,93]. We checked that Br(H 1 → γγ) ≤ 1 × 10 −3 for most samples and 2 In doing [105], we once confronted with the situation quite similar to what we are facing now. Our detailed simulation at that time indicated that the trilepton constraint on SUSY is very weak. Moreover, in comparison with the case discussed in [118], we find that our case is more difficult to detect since the signal is smaller. Br(A 1 → γγ) < 4 × 10 −4 for all samples. So current results of the Fermi-LAT search for spectral lines [119] can not impose tight limit on Solution II (see [83] for a detailed discussion).

Solution III -the W + W − annihilation channel
In general, the pure annihilationχ 0 1χ 0 1 → W + W − is unable to explain the GCE quite well [73,79], but if it mixes sizably with other annihilation channels, the generated spectrum may be improved significantly to account for the GCE. Solution III in the NMSSM belongs to this case.
Solution III has the following features: • The W pair must be produced close to threshold to account for the GCE, which means 85GeV mχ0 1 100GeV (see right panels of Fig.3).
• From the expression of σ W W v 0 in Eq.(2.4), one can learn that, if the wino is decoupled, the annihilation rate is determined by the higgsino-dominated chargino. In this case, we have In getting these expressions, we note v s ≡ µ/λ 400GeV (see lower left panel of Fig.2), and expand N 13 and N 14 in terms of λv/µ (see Appendix). We also use the approximation mχ0 1 2κµ/λ. Then σ W W v 0 ∼ 10 −26 cm 3 /s and mχ0 1 ∼ 90GeV limit tightly the ranges of the parameters λ, κ and µ, which are shown in Table 1 and Fig.3.
Note in Solution III, the parameter µ, or equivalently the masses for the higgsinodominated chargino and neutralinos, is less than about 150GeV. Since the splitting between µ and mχ0 1 is less than about 50GeV, such a low value of µ is still allowed by the LHC search for SUSY (see footnote 2 in our discussion on Solution II).
• The upper left panel of Fig.3 indicates that the parameters λ and κ are correlated by As a result, we have mχ0 1 2µ/3.
• As we emphasized before, the annihilationχ 0 1χ 0 1 → W + W − must mix sizably with the annihilationχ 0 1χ 0 1 → H 1 A 1 to explain the GCE. This, in return, requires appropriate masses for H 1 and A 1 to improve the γ-ray spectrum generated by the W W state. In Fig.4, we plot the GCE χ 2 as a function of DM mass in Solution III with different colors denoting the branching ratio of the DM annihilation into H 1 A 1 . This figure indicates that, with the increase of the branching ratio, the GCE χ 2 tends to decrease.

GCE solutions with H 1 being the SM-like Higgs boson
In this section, we investigate the GCE solutions for the case that DM is singlino-dominated, and meanwhile H 1 acts as the SM-like Higgs boson. We carry out our study in a way similar to what we did in Section III.
In Fig.5 we project the solutions on σv 0 − mχ0 1 plane (upper panel) and χ 2 GCE − mχ0 1 plane (lower panel). For solutions marked by green lozenge, DM annihilates with the largest branching ratio into H 1 A 1 , while for those marked by red pentastar, DM annihilates mainly into H 2 A 1 . In the following, we call these two kinds of solutions Solution IV and Solution V respectively. Fig.5 then indicates that, for 70GeV ≤ mχ0 1 ≤ 87GeV, Solution IV can explain the GCE quite well with the best explanation coming from mχ0 1 81GeV with χ 2 GCE 21.4 (corresponding to a p-value of 0.55), and for 80GeV ≤ mχ0 1 ≤ 130GeV, Solution V is good in accounting for the GCE with the best explanation locating at mχ0 1 85GeV with χ 2 GCE 21.6 and a p-value of 0.54. Compared with the case that H 2 acts as the SM-like Higgs boson, we find that it is more difficult to get the GCE solutions if H 1 corresponds to the SM-like Higgs boson. One important reason is that the spectrum of the singlet-dominated particles for Solution IV and V has non-trivial requirements on the NMSSM parameters, which can not be easily satisfied due to the structure of the NMSSM itself. A good example about this argument is that we do not find any solutions where DM mainly annihilates into bb. This is due to the fact that, given a singlino-dominated DM with 30GeV ≤ mχ0 1 ≤ 70GeV and meanwhile a singlet-dominated A 1 satisfying m A 1 2mχ0 1 , the singlet-dominated CP-even Higgs boson is usually lighter than the SM-like Higgs boson [69].
In Table 2 Fig.6, we show the correlations among the parameters λ, κ, µ, mχ0 1 and m A 1 . This figure is supplement to Table 1, and as we will show below, it is helpful for our

Solution IV -the H 1 A 1 annihilation channel
Solution IV has the following features: • The H 1 A 1 state must be produced close to threshold to explain the GCE, which is reflected by δ < 0.1 from our results.
• The favored spectrum for the singlet-dominated particles is 71GeV mχ0 1 87GeV, 10GeV m A 1 40GeV and 126GeV m H 2 142GeV. Given κ ∼ 0.12 which is required to predict the right size of the annihilationχ 0 1χ 0 1 → H 1 A 1 for the GCE (see below), this spectrum limits parameters such as λ, κ, µ and A κ within rather narrow ranges, which are given in Table 1 and also shown in Fig.6.
Compared with Solution II, we find in Solution IV that, in order to predict a heavier singlet-dominated CP-even Higgs boson, the parameter µ usually takes a larger value, 220GeV µ 270GeV. As a result, λ must exceed about 0.6, which can be inferred from the relation mχ0 1 2κµ/λ 2 × 0.12 × µ/λ 80GeV. This relation also suggests that v s ≡ µ/λ 300GeV or λv/µ < 0.6, which makes the expansions listed in Appendix feasible.
• Similar to Solution II, the s-channel contributions to the annihilation rate σ H 1 A 1 v 0 in Eq.(2.5) are usually much smaller than those from the t/u channel, and among the t/u channel contributions, the one induced by the exchange ofχ 0 1 is dominant. However, since H 1 now is the SM-like Higgs boson (instead of a singlet-dominated particle in Solution II), there still exists a slight difference between the two solutions. Explicitly speaking, we find that each higgsino contribution to the annihilation is comparable in magnitude with theχ 0 1 contribution, but since the two higgsino contributions cancel each other, the total higgsino contribution is small. These features can be explained by the following formula (see Eq.(A.20)) and also by comparing Eq.(4.1) with equation (3.20) in [93] to conclude that κ ∼ 0.12 is enough to predict theχ 0 1 contributed σ H 1 A 1 v 0 at the order of 10 −26 cm 3 /s.
• Since m A 1 40GeV for all cases in Solution IV (see upper right panel of Fig.6), the SM-like Higgs boson H 1 will decay into A 1 A 1 . We checked that Br(H 1 → A 1 A 1 ) 24% as required by the Higgs data.
• We also checked that A 1 → bb is the dominant decay mode of A 1 , and Br(A 1 → γγ) < 5 × 10 −5 for all samples.

Solution V -the H 2 A 1 annihilation channel
Since H 2 in Solution V is singlet dominated, the features of Solution V should be similar to those of Solution II. The differences mainly come from the following aspects: • The spectrum of the singlet dominated particles. In Solution V, the favored spectrum is 80GeV mχ0 1 130GeV, 18GeV m A 1 100GeV and 125GeV m H 2 146GeV with m A 1 < mχ0 1 < m H 2 and δ < 0.1. Corresponding to such a spectrum, the parameter space of Solution V differs greatly from that of Solution II, which can be seen from Table 1 and also from Fig.6.
• The phenomenology of some relevant particles. For example, in both Solution IV and Solution V, the favored value of µ is uplifted in comparison with that in Solution II. As a result, the higgsino-dominated neutralinos may decay into Zχ 0 1 , which makes them to be potentially detected at 14-TeV LHC by trilepton +E miss T signals [69].

Explore the GCE solutions in future DM experiments
In this section we investigate to what extent the GCE solutions will be explored in future DM direct detection experiments such as XENON-1T and LUX experiments [120], which will improve current experimental sensitivities to DM-nucleon scattering cross sections by up to three orders. In Fig.7, we project our solutions on mχ0 1 − σ SI p and mχ0 1 − σ SD p planes with σ SI p and σ SD p denoting the spin-independent (SI) and spin-dependent (SD) cross sections respectively. The left panels in the figure are the results for the case that H 2 acts as the SM-like Higgs boson, and the right panels are those for the case that H 1 corresponds to the SM-like Higgs boson. The dotted lines, solid lines, dashed lines and dash dotted lines are the sensitivities to the cross sections set by the XENON-100, LUX, XENON-1T and LZ experiments respectively. Note that so far the XENON-100 experiment has imposed constraints on both SI and SD cross sections, while the LUX experiment only obtained limits on the SI cross section.
For σ SI p in the H 2 case, we can see from Fig.7 that the future XENON-1T experiment is able to probe a large portion of the GCE solutions, and the LZ experiment can test even more solutions. Anyhow, there still exist some solutions remaining untouched by these future experiments. This conclusion can be understood as follows. In the NMSSM after considering the current experimental constraints on sfermion masses, the main contribution to σ SI p comes from the t-channel process mediated by the CP-even Higgses H 1,2 . In this case, the Wilson coefficient f q i for the operatorχ 0 1χ 0 1q i q i is given by [121] f  1) indicates that, if κ is small or if there exists a strong cancelation between the two terms, f q i or equivalently the SI cross section will be suppressed. We numerically checked that the untouched solutions has either of the two characteristics.
On the other hand, the story for σ SD p in the H 2 case is quite different. From the lower left panel of Fig.7 we can see that the future XENON-1T experiment can test almost all of the GCE solutions, let alone the more sensitive LZ experiment. The underlying reason is that in the NMSSM with heavy sfermions, the SD cross section gets contribution mainly from the t-channel Z-mediated diagram. As a result, the size of the cross section is determined by the Zχ 0 1χ 0 1 coupling, which is given by In getting this expression, we have used the approximations for N 13 , N 14 and mχ0 1 . Then from the results presented in Fig.3, one can infer that except for some rare cases of Solution I, the SD cross section is not suppressed too much.
In a similar way, one can analyze the results for the H 1 case, which are shown on the right panels of Fig.7. For example, the upper right panel indicates that the SI cross sections in Solution IV and Solution V are usually larger than 10 −10 pb. This may be understood by a weak cancelation between the two terms in Eq.(5.1). Compared with the H 2 case, both the SI cross section and the SD cross section in the H 1 case are large and consequently, all the solutions will be tested by XENON-1T experiment.
In principle, the GCE solutions in the NMSSM may also be tested by electroweakino production processes at the LHC [69]. We will discuss such an issue in our forthcoming work.

Summary
In this work, we took into account the recently reported foreground and background uncertainties for the GCE and investigated its explanation by DM annihilation in the framework of the NMSSM. We carried out a sophisticated scan over the NMSSM parameter space by considering various experimental constraints such as the Higgs data, B−physics observables, DM relic density, LUX experiment and the dSphs constraints. Then for each surviving parameter point we performed a fit to the GCE spectrum by using the correlation matrix that incorporated both the statistical and systematic uncertainties of the measured excess. Our results indicate that due to the introduction of the gauge singlet Higgs superfield, the NMSSM with a singlino-dominated DM has multiple DM annihilation channels that are able to explain the GCE quite well, and all of these explanations require the singlet-dominated particles (including one neutralino, one CP-even and one CP-odd Higgs bosons) to be moderately light. We also discussed to what extent the future DM direct detection experiments can explore the GCE solutions, and we conclude that the XENON-1T experiment is very promising in testing nearly all the solutions.
When choosing the scenario of particle spectrum, we focused on a singlino-dominated DM and considered the cases that either H 2 or H 1 acts as the SM-like Higgs boson. For the popular situation that H 2 corresponds to the SM-like Higgs, we have the following observations on the GCE solutions: • For the annihilationχ 0 1χ 0 1 → bb, DM mass is now allowed in the range from 30GeV to 70GeV which is much wider than before. With the help of an appropriate schannel resonance, the singlet trilinear self-coupling parameter κ can be as low as 0.02 to explain the GCE. Moreover, the higgsino mass parameter µ is upper bounded by about 300GeV to ensure a correct DM relic density. Since there exist strong correlations between independent parameters, such an explanation suffers from a fine tuning problem, which is usually less than 1%.
• The annihilationχ 0 1χ 0 1 → H 1 A 1 may provide a better explanation than the channel χ 0 1χ 0 1 → bb when H 1 A 1 is produced close to threshold, and the best interpretation corresponds to a p-value of 0.55. In this kind of explanation, the singlet-dominated particles must satisfy 60GeV mχ0 1 115GeV, 10GeV m A 1 110GeV, 60GeV m H 1 120GeV and δ < 0.2. This imposes non-trivial constraints on the NMSSM parameters, especially that µ must be less than about 200 GeV. Among various contributions to the annihilation, the dominant one comes from theχ 0 1 -contributed t/u channel diagrams, in which the parameter κ plays an important role in deciding the annihilation rate.
• Apart from the necessary mixing with the H 1 A 1 final states, W + W − pair in the annihilationχ 0 1χ 0 1 → W + W − must be produced close to threshold to account for the GCE. A small µ less than about 150 GeV is necessary to increase the annihilation rate through the t/u-channel contributions induced by a higgsino-dominated chargino. The LHC search for trilepton +E miss T signal can not exclude such a possibility since the electroweakino production rates at the LHC are relatively low, and meanwhile since the splitting between µ and mχ0 1 is compressed.
• The detection of spin-independent scattering in the future XENON-1T and LUX experiments are able to cover a large portion of the GCE-favored parameter space, while the spin-dependent detection have a stronger potential to test nearly all of the relevant parameter region.
As for the case that H 1 acts as the SM-like Higgs boson, the features of the GCE solutions are quite different, which are as follows: • In comparison with the H 2 case, it is difficult to find GCE solutions when H 1 corresponds to the SM-like Higgs boson, and especially we did not find any solution that DM annihilates mainly into bb. The reason is, assuming H 1 to be the SM-like Higgs boson, there must exist sizeable mass splittings among the light singlet-dominated particles to explain the GCE, which is difficult to realize in the NMSSM due to the theoretical structure itself.
• For 80GeV mχ0 1 86GeV, DM may annihilate into H 1 A 1 and H 2 A 1 states with comparable rates to explain the GCE, while for mχ0 1 80GeV (mχ0 1 100GeV), the sole annihilationχ 0 1χ 0 can be responsible for the GCE. For all these solutions, the singlet-dominated particle H 2 and the parameter µ must satisfy 125GeV m H 2 145GeV and 210GeV µ 270GeV.
• Both the spin-independent and spin-dependent detection in the future XENON-1T experiment have a great potential to test the relevant parameter space.
Before we end our discussion, we would like to comment briefly on the interpretation of the GCE with a bino-like DM. Like the singlino-dominated DM case, a light A 1 with mass below about 140GeV is necessary for such a work, and this A 1 prefers to be singletdominated 3 . The difference is that, for the bino-like DM case, the interaction of the DM with A 1 is relatively small and consequently the annihilationχ 0 1χ 0 1 → H i A 1 can not explain the GCE any more due to its rather low annihilation rate. Also due to the suppressed interaction, m A 1 must be closer to 2mχ0 1 for the annihilationχ 0 1χ 0 1 → bb to account for the GCE, and thus the theory has to be tuned in a more elaborated way. Our sophisticated scan over the relevant NMSSM parameter space verified these conclusions.  Table 2 of [123]). In either case, mA 1 should be significantly smaller than m H ± to escape experimental constraints. Previous studies have suggested that a light doublet-dominated A1 might also explain the galactic center excess. However, due to the requirements on the elements this scenario occurs only in specific portions of the parameter space and is significantly more experimentally constrained than those we considered. In fact, in our scans for the GCE we did not find any parameter points with the doublet component of the light A1 exceeding 0.1. In summary, a light doublet-dominated A1 may exist, as suggested by e.g. Ref. [123], but it is fair to say that without a very delicate parameter tuning, it is difficult to obtain in explaining the GCE, especially when one considers more constraints than previous literatures. About this conclusion, we thank the authors of [123] for helpful discussion.

A Properties of the singlet-dominated particles
In this appendix, we present some analytic expressions for the masses and couplings of the singlet dominated particles, such asχ 0 1 and A 1 in the NMSSM. These expressions are obtained by diagonalizing the mass matrices of the particles (like done in [57]), and are good approximations in certain cases. They are helpful in understanding the results presented in this work. In the following, we will follow notations and conventions consistent with [104] for the Z 3 NMSSM.

A.1 Neutralino masses and mixings
In the basis ψ 0 = (−iλ 1 , −iλ 3 2 , ψ 0 d , ψ 0 u , ψ S ), the neutralino mass matrix is: If the bino and wino fields are decoupled, the mass eigenstates of the neutralinos can be approximated byχ whereχ 0 1 denotes the lightest neutralino with ψ S field as its dominant component in this work, andχ 0 i represents a higgsino-like neutralino. In the limit of |µ| λv, 1 κ/λ and tan β 1, one can expand the neutralino masses and N ij by the power of λv/µ ≡ v/v s to get the following approximations: In above expressions, the Sgn and θ functions are defined by Likewise, one may consider the case that the wino and the singlino fields decouple. In this case, the mass eigenstates of the neutralinos is approximated bỹ In the limit of tan β 1, |µ| g 2 v u and |µ| M 1 , we have the following approximations:

A.2 CP-odd Higgs mass matrix
In the (A, S I ) "interaction" basis, the mass matrix for CP-odd Higgs bosons is given by In the case of m A max(v, |A κ |, |µ|), κ/λ 1 and tan β 1, the lighter CP-odd scalar A 1 is singlet dominated with its squared mass given by This approximation indicates that, without considering the radiative corrections, the singletdominated CP-odd scalar mass is determined by the parameters λ, κ, µ as well as A κ . The components of A 1 can be written as where P A 1 ,A is the active component and P A 1 ,S I is the singlet component of the A 1 .

A.4 Some properties of the singlet-dominated particles
With the assumptions that M A max(|µ|, |A κ |), |µ| λv, tan β 1 and κ/λ 1, one can approximate the masses and couplings of the singlet dominated particles, such asχ 0 1 and A 1 , by simple analytic expressions [57]. In the following, we list some of the coupling expressions used in our discussion, which are denoted by C XY Z hereafter. These expressions are actually expand the corresponding exact ones by the power of λv/µ.