Interpretations of galactic center gamma-ray excess confronting the PandaX-II constraints on dark matter-neutron spin-dependent scatterings in the NMSSM

The Weakly Interacting Massive Particle (WIMP) has been one of the most attractive candidates for Dark Matter (DM), and the lightest neutralino (χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\chi }^0_1$$\end{document}) in the Next-to-Minimal Supersymmetric Standard Model (NMSSM) is an interesting realization of the WIMP framework. The Galactic Center Excess (GCE) indicated from the analysis of the photon data of the Fermi Large Area Telescope (Fermi-LAT) in the gamma-ray wavelength ≲1fm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim 1 \,\mathrm{fm}$$\end{document}, can be explained by WIMP DM annihilations in the sky, as shown in many existing works. In this work we consider an interesting scenario in the Z3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_3$$\end{document}-NMSSM where the singlet S and Singlino S~0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{S}^0$$\end{document} components play important roles in the Higgs and DM sector. Guided by our analytical arguments, we perform a sophisticated scan over the NMSSM parameter space by considering various observables such as the Standard Model (SM) Higgs data measured by the ATLAS and CMS experiments at the Large Hadron Collider (LHC), and the B-physics observables BR(Bs→Xsγ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$BR(B_s\rightarrow X_s\gamma )$$\end{document} and BR(Bs→μ+μ-)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$BR(B_s\rightarrow \mu ^+\mu ^-)$$\end{document}. We first collect samples which can explain the GCE well while passing all constraints we consider except for the DM direct detection (DD) bounds from XENON1T and PandaX-II experiments. We analyze the features of these samples suitable for the GCE interpretation and find that χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\chi }^0_1$$\end{document} DM are mostly Singlino-like and annihilation products are mostly the bottom quark pairs b¯b\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{b}b$$\end{document} through a light singlet-like CP-odd Higgs A1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_1$$\end{document}. Moreover, a good fit to the GCE spectrum generically requires sizable DM annihilation rates ⟨σbb¯v⟩0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \sigma _{b\bar{b}} v \rangle _{0}$$\end{document} in today’s Universe. However, the correlation between the coupling CA1bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{A_1 b\bar{b}}$$\end{document} in ⟨σbb¯v⟩0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \sigma _{b\bar{b}} v \rangle _{0}$$\end{document} and the coupling CZχ~10χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{Z \widetilde{\chi }^0_1 \widetilde{\chi }^0_1}$$\end{document} in DM-neutron Spin Dependent (SD) scattering rate σχ~10-NSD\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^{SD}_{\widetilde{\chi }^0_1-N}$$\end{document} makes all samples we obtain for GCE explanation get excluded by the PandaX-II results. Although the DM resonant annihilation scenarios may be beyond the reach of our analytical approximations and scan strategy, the aforementioned correlation can be a reasonable motivation for future experiments such as PandaX-nT to further test the NMSSM interpretation of GCE.


Introduction
An excess of gamma-rays in the direction of the Galactic Center (GCE) has been reported by several groups analyzing the data from the Large Area Telescope on board the Fermi Gamma-ray Space Telescope (Fermi-LAT) [1][2][3][4][5][6][7][8][9][10][11][12]. Various interpretations have been proposed to provide additional gamma-ray sources which can be roughly classified into two categories. One is of astrophysical origin such as a large population of unresolved millisecond pulsars (MSPs) in the Galactic bulge [13][14][15][16][17][18][19][20][21] or a series of recent leptonic cosmic-ray outbursts [22][23][24]. Another category is of particle Dark Matter (DM) origin which annihilate into the Standard Model (SM) particles subsequently producing gamma-rays in the excess energy range (see e.g.  and references therein). For the comparison between the two categories of GCE interpretation, it has been argued that [47] a sufficiently large population of MSPs would imply a large population of observable low-mass X-ray binaries, limiting the contribution of MSPs to the GCE to ∼ 4-23% and thus leaving annihilating DM as a contender. Another argument comes from [18] by comparing the GCE spectrum to that measured from 37 MSPs by Fermi experiment, which indicates that a population of unresolved MSPs exhibit a spectral shape that is too soft at sub-GeV energies to accommodate the GCE. Furthermore, [20] argued that after pulsars are expelled from a globular cluster, they continue to lose rotational kinetic energy and become less luminous. This makes the luminosity function of those MSPs depart from the steady-state distribution and thus MSPs born in globular clusters can account for only a few percent or less of the GCE. Consequently, pulsars located in the Galactic bulge or a Galactic halo of DM are both viable interpretations for the GCE.
The existence of DM has been verified by various cosmological and astrophysical observations. However, no direct evidence of DM particle interactions has been measured by any experiment. DM particle search is a promising strategy to find new particle physics principles Beyond the Standard Model (BSM). Weakly Interacting Massive Particle (WIMP) has been a popular DM candidate given that its scatterings off the SM nuclei mediated by the weak interaction may be detected [48,49]. One of the most attractive WIMPs is the Lightest Supersymmetric Particle (LSP) in supersymmetric models (SUSY) such as the lightest neutralino ( χ 0 1 ). Among the various realizations of the SUSY framework, the Minimal Supersymmetric Standard Model (MSSM) is the simplest one and provides a promising WIMP DM. However, after a series of experiments upgraded their sensitivities and set stronger bounds, MSSM is no longer a favorable framework for WIMP. Firstly, the SM-like Higgs mass measured to be around 125 GeV by the ATLAS and CMS experiment at the Large Hadron Collider (LHC) [50,51] requires a sizable loop radiative corrections through heavy top squark ( 1 TeV) which leads to a large fine tuning. Secondly, to generate the observed DM relic density while surviving the DM direct detection limits simultaneously, the Bino ( B 0 ) must be the dominant ingredient in the lightest neutralino. Consequently, the μ parameter in MSSM have to be relatively large which further aggravates the fine tuning problem.
As a minimal extension of the MSSM, the Next-to-MSSM (NMSSM) contains an additional gauge singlet superfieldŜ which can dynamically generate the μ parameter in MSSM with a small value, in which case the fine tuning problem can be alleviated significantly [52]. Recently, several experiments including LUX [53], PandaX [54] and XENON [55] updated their results of DM direct detections (DD). At this moment, the strongest constraints on Spin-Independent (SI) and Spin-Dependent (SD) WIMP-nucleon scattering cross sections come from XENON1T (σ S I χ 0 1 -P ) and PandaX-II (σ S D χ 0 1 -N ), respectively. These stringent limits have already pushed the neutralino WIMP DM scenario in NMSSM into a challenging situation: one either imposes large cancelations in the DM-nucleon interactions among different contributions, the so-called blind spot scenario [56][57][58][59][60][61], or chooses a large value of μ parameter as did in the MSSM [60]. Either choice makes the NMSSM also suffer from the difficulty of fine tuning.
In [46] we observed that a light neutralino χ 0 1 in the Z 3 -NMSSM was able to explain the GCE meanwhile satisfying various other constraints including the observed DM relic density and the Higgs data measured by the ATLAS and CMS groups at the LHC, by adjusting the χ 0 1 pair annihilation branch ratios into bb, W + W − and A 1 H i (i = 1, 2). These DM annihilations can generate additional gamma-ray sources strong enough to interpret the observed GCE and present a close connection between the astrophysical gammaray observations and DM phenomenology. More interestingly, the LSP χ 0 1 usually manifests itself as large missing transverse momentum at high energy colliders such as the LHC, which is one of the most important signals for searching SUSY.
In this work we present our dedicated analysis on a correlation pattern between the DM annihilation cross sections to explain the GCE and DM-nucleon spin-dependent scattering rates in the direct detection experiments, in the Z 3 -NMSSM where the singlet S and Singlino S 0 components play important roles in the Higgs and DM sector. Guided by the analytical arguments, we perform a sophisticated scan by employing the Markov Chain Monte Carlo (MCMC) strategy over the NMSSM parameter space, to verify our expectations and explore whether numerical calculations can reveal NMSSM parameter space that can still explain the GCE while passing the limits from XENON1T and PandaX-II.
We organize the paper as follows. In Sect. 2 we recapitulate the main features of neutralino DM in the NMSSM and the DM annihilation mechanisms to explain the GCE. Using the analytical approximations of spin-dependent DMnucleon scattering rate, we identify the correlation between the GCE interpretation and the scattering strength. In Sect. 3 we explain our scan strategies in detail. We present and analyze the results of our scan in Sect. 4 and summarize in Sect. 5.

Neutralino DM in the NMSSM
In this section, we briefly recapitulate the features of neutralino DM in NMSSM and the mechanisms to explain the GCE in terms of DM annihilations. Then we discuss the analytical approximations of neutralino DM-nucleon spindependent scattering rate, based on which we identify the correlation between the GCE interpretation and the scattering strength in direct detection.

Lightest neutralino in NMSSM as WIMP DM
The NMSSM introduces a gauge singlet Higgs superfieldŜ in additional to the two doublet Higgs superfieldĤ u ,Ĥ d in MSSM. A Z 3 symmetry is adopted to the superpotential of NMSSM to forbid the appearance of parameters with mass dimension. In this case, the superpotential of NMSSM is given by [62]: where W F is the superpotential of MSSM without the μterm and λ, κ are dimensionless parameters that describe the interactions among the NMSSM Higgs superfields. The soft SUSY-breaking Higgs potential is: As a result, the Higgs sector of NMSSM has six free input parameters: NMSSM neutralino sector is obtained from the mixings among five gauginos: Bino B 0 , Wino W 0 , Higgsinos H 0 d,u and Singlino S 0 . The symmetric neutralino mass matrix in the basis of (ψ where M 1 , M 2 are soft SUSY-breaking masses of Bino and Wino, respectively. After diagonalization one obtains five mass eigenstates χ 0 i = N i j ψ 0 j with ascending order of mass values, where N i j is the rotation matrix. The lightest neutralino χ 0 1 is identified as the WIMP DM in our discussion, which can be written as

Explaining GCE with NMSSM neutralino DM
The main goal of identifying the lightest neutralino in NMSSM as WIMP DM is to use their annihilations in today's Universe to provide additional gamma-ray source in the Galactic Center (GC) region to explain the GCE. Therefore an adequate cross section of present neutralino annihilations is necessary. One can express the relevant annihilation processes as In [46] we observed that the relevant final states XY can be bottom quark pair bb, W boson pair W + W − , and two scalars A 1 H 1,2 where H 1,2 and A 1 are CP-even and CP-odd Higgs bosons, respectively. These SM particles as DM annihilation products will further generate photon spectrum including the gamma-ray energy range through decay, parton showering and hadronization effects [63]. Taking the XY = bb in [46] as an example, the DM annihilation cross section σ bb v 0 in today's Universe with v → 0 can be approximated as [64] where and C A i bb are couplings. A i with i = 1, 2 are the CP-odd Higgs and Γ A i are the widths of A i . We learned from [46] that the s-channel A i -mediations dominate the χ 0 1 pair annihilation into bb.
The predicted gamma-ray spectrum from DM annihilations in our Galaxy that are observed near Earth can be expressed as the following differential flux for a given angular direction: where d N XY γ /d E γ is the prompt gamma-ray spectrum from a single DM annihilation event into XY . J is an astrophysical factor reflecting the Galactic distribution of DM obtained by integrating over ρ(r ) 2 along the line of sight (l.o.s.) with ρ(r ) being the Navarro-Frenk-White (NFW) DM profile [65]. We choose the following parameters in ρ(r ) as the same ones chosen in the Fermi-LAT analyses [12,66], which are consistent with various astrophysical observations ( e.g. the measured gravitational potential of the inner Galaxy [9]) while generally providing a good fit to the observed GCE [3,44].
where s is the l.o.s distance and θ is the angle between the l.o.s and the axis connecting Sun and GC.
2.3 Correlation of GCE explanation to neutralino DM-neutron SD scattering in NMSSM XY = bb in Eq. (7) has been indicated in many existing works as a characteristic channel of DM annihilation to interpret GCE. Equation (8) shows that the corresponding DM annihilation rate is approximately proportional to the coupling C A 1 bb in NMSSM with a generally heavy A 2 [46,62], where y b is the SM bottom quark Yukawa coupling and |P 1,A | 2 is the ingredient in the NMSSM CP-odd Higgs A 1 from the would-be CP-odd Higgs A in MSSM with mass given in Eq. (3). For a light singlet-like A 1 , one usu- In the case of relatively large tan β O(5) which is common in the NMSSM [62] and suitable for GCE explanation [46], we can have the following approximation as long as P 1,A is small but stable for a light singlet-like A 1 , which is the case we consider.
In [46] we observed the potential of future DM direct detection experiments to test the GCE explanation scenario in NMSSM, especially in terms of the DM-neutron spin dependent scattering rate σ S D χ 0 1 −N (see Fig. 7 therein). The dominant contribution comes from the t-channel mediation by Z boson where the coupling C χ 0 1 q can be approximated in the limit of non-relativistic scattering as follows where v 174 GeV is the SM Higgs vev and N i j is rotation matrix for neutralino mass diagnolization in Eq. (6). N i j can have simple approximated expressions when DM χ 0 1 is dominated by one or two components (see Appendix of [46] where |N 15 | 2 O(0.5) applies to the Singlino-like χ 0 1 , as the case in most of our samples. One can see that tan β = 1 breaks the symmetric structure with respect to the Higgsinos H 0 d , H 0 u in the neutralino mass matrix Eq. (5) and generates non-zero coupling between Z boson and χ 0 1 . More importantly, tan β > 1 in NMSSM can increase |C Z χ 0 1 χ 0 1 | and the resulting DM-neutron SD scattering rate [67] In most samples we study in this work, the Singlino component dominates in χ 0 1 . However, we note that if additional components such as Bino B 0 and Wino W 0 also contribute sizably to χ 0 1 , the above approximation may not apply very well.
To summarize this subsection, we obtained an approximately quantitative connections in terms of sizable tan β O(5) in NMSSM with a singlet-like A 1 and a Singlino-like χ 0 1 , between the GCE explanation and DM-neutron spindependent scattering rate as follows

Scan strategies
The scan over the NMSSM parameter space was performed by employing the MCMC method. We utilize NMSSMTools-5.0 [68,69] to calculate the NMSSM mass spectrums, particle decaying ratios and B-physics observables, which are further interfaced to MicrOMEGAs [67,70] to calculate the DM annihilation cross section, relic density, DM-nucleon scattering rates, and the gamma-ray spectrum to compared with GCE. We build an overall χ 2 function with which smaller values reflect better compatibility of a model sample in the NMSSM parameter space with the observables contained in χ 2 . The MCMC scan strategy is powerful in finding the NMSSM parameter regions meeting our physical goals.

NMSSM parameter simplifications and ranges
In our scan, NMSSM parameters not closely relevant to the GCE discussion are chosen with simplifications. The first and second generation squarks are required to decouple with masses fixed at 2 TeV. For up-and down-type squarks of the third generation, we simplify the soft-SUSY breaking parameters to be A t = A b and M U 3 = M D 3 . All parameters in the slepton sector are chosen to be equal, including their masses and soft-SUSY breaking couplings. The gluino mass parameter M 3 is also fixed at 2 TeV. To summarize, the ranges of free NMSSM parameters in the scan are listed below: 2 < tan β < 60, 0.001 < λ < 1, 0.001 < κ < 1, 100 GeV < μ < 1.5 TeV, 20 GeV < M 1 < 1 TeV, 100 GeV < M 2 < 2 TeV, |A top | < 6 TeV, 100 GeV < M Q 3 < 2 TeV, 100 GeV < M U 3 < 2 TeV, 100 GeV < M L = M e = A e < 2 TeV.

χ 2 method in MCMC strategy
Now we briefly recapitulate the χ 2 method based on which we embed our codes of MCMC chain into NMSSMTools-5.0. Smaller χ 2 of a NMSSM sample can be regarded as its better compatibility confronting the experimental observations built into the χ 2 . Smaller χ 2 also result in higher acceptance probability of this sample which will act as the new starting point to search nearby parameter space better fitting the experiments. Note that a model point may still be accepted as a node in the scanning chain even if it predicts one or several observables outside the experimental bounds at some confidence level, as long as its χ 2 is not very large. This strategy can help the codes avoid to get stuck in some corners of the total parameter space. When χ 2 reaches a stable minimum after a long scan chain, we will apply a set of selection rules with respect to the relevant observables and name the kept points as GCE-samples.
In this work we construct an overall χ 2 total in Eq. (22) which consists of the SM-like Higgs mass, B-physics observables, DM relic abundance, DM-nucleon scattering rates, and the shape of GCE spectrum.
In the above, F i , m h SM , B R, Ω χ 0 1 h 2 are weight factors, the SM-like Higgs mass, branching ratio and DM relic abundance, respectively. During the scan we dynamically adjust F i to explore as wide regions as possible of the NMSSM parameter space, otherwise χ 2 of one observable may overwhelm others' effects. σ S I h 2 , correspond to the case in which an observable has a theoretical prediction μ i , an observed center value μ 0 and uncertainty σ i (including both experimental and theoretical parts). In this case the χ 2 i is defined as in which the values of μ i , μ 0 , σ i will be provided when we apply the final-step selection cuts. according to [71]: with δσ S I χ 0 1 −P = σ S I,0 where σ S I,0 -Lastly, χ 2 GC E in Eq. (22) is constructed according to Eq. (5.1) of [3]: where dN /d E i = Ad N/d E i and d N 0 /d E i are the gamma-ray spectrum predicted by NMSSM samples with a scaling factor A and GCE spectrum extracted from Fermi-LAT data after background modeling, respectively, in the i-th gamma-ray energy bin. A is a tuning factor accommodating the Galactic DM halo profile uncertainties which scales the theoretical predictions of d N/d E i . We define χ 2 GC E as the minimum value of χ 2 (A) with varying A [46], where the average runs over the 22 test region of interests (ROIs) in [3]. With χ 2 GC E constructed in Eq. (28), all terms in Eq. (22) are consistent in the sense of statistic definition. We refer interested readers to [3] for more details about Σ i j .
When χ 2 total in the scan reaches a relatively stable minimum, we upgrade the requirements related to the observables in χ 2 total with hard cutting ranges from the corresponding experiments. Other constraints we further apply include the DM indirect search results from dwarf spheroidal galaxies (dSph) and the bounds from SUSY searches at LEP and LHC. More specifically, we consider: 1. The SM Higgs data [50,51,72]. We require NMSSM to accommodate a CP-even Higgs boson whose mass is near 125 GeV while decays and couplings satisfy the observed results by ATLAS and CMS. We implement this requirement by utilizing the packages HiggsBounds-5.0.0 [73] and HiggsSignal-2.0.0 [74]. The Higgs mass we use is m h SM = (125.36 ± 0.41 ± 2.0) GeV [72], indicating in order the experimental central value μ i in Eq. (23), the experimental uncertainty σ i,the at 1σ and the theoretical uncertainty σ i,ex p at 1σ . Thus the total uncertainty in Eq.
(23) is σ i = σ 2 i,the + σ 2 i,ex p , which also applies to the following observables with similar format of data. 2. B-physics constraints, i.e. B R(B s → X s γ ) = (3.43 ± 0.22 ± 0.24) * 10 −4 [75] and B R(B s → μ + μ − ) = (2.9 ± 0.7 ± 0.29) × 10 −9 [76], with the same form as m h SM . 3. Constraints from SUSY searches at the LEP and LHC, which have been encoded in the package NMSSMTools-5.0, such as the mass limits on SUSY particles. 4. Bounds on DM annihilations from Fermi-LAT gammaray observations of Milky Way dSphs [77] in which likelihood analysis codes are provided for constraining theoretical models. 1 We input the predicted gamma-ray spectrum and DM annihilation cross sections of our NMSSM samples and follow the code analysis. We apply the selection cut χ 2 d Sph < 2.71/2 suggested by the codes to pick out the NMSSM samples passing the dSph constraints. 5. Observed DM relic density Ω χ 0 1 h 2 = 0.1197×(1±10%) [78], where the 10% is included to accommodate the uncertainty of numerical calculations in MicrOMEGAs [70]. 6. For the GCE spectrum comparison between the theoretical and experimental results, we impose the final selection criterion χ 2 GC E < 35.2. This value corresponds to the explanation of GCE by our NMSSM samples at 95% confidence level with N ex p. − N the. = 24 − 1 degree of freedom (d.o.f) [46], where N ex p. = 24 is the number of GCE energy bins in the analysis of Fermi-LAT data [3] and N the. = 1 is the theoretical scaling factor A in Eq. (28).
Note that the data from DM direct detection experiments XENON1T and PandaX-II are included in Eq. (22) for the χ 2 total scan to explore NMSSM parameter space with generally small scattering rates σ S I to be under the exclusion bounds. As we mentioned in the beginning of this subsection, samples surviving all the above selection rules which can explain the GCE may not satisfy the bounds of XENON1T and PandaX-II. We name the samples obtained at this stage the GCE-samples. In the next section, we will first analyze the characteristics of the GCEsamples. Then we confront them with the latest constraints from XENON1T and PandaX-II.

Results from the scan
In this section we first present the main features of the GCEsamples obtained in Sect. 3. Then we discuss their com- .. are not accessible for DM so light in today's Universe, we checked that the dominant DM annihilation channel is χ 0 1 χ 0 1 → bb. However, we checked that some samples with χ 0 1 χ 0 1 → W + W − also passed various requirements in Sect. 3 but failed to satisfy GCE fitting criterion χ 2 GC E < 35.2, which has been noticed and pointed out in our previous work [46] with a different construction of χ 2 total . The difficulty of W + W − channel to fit the GCE comes from the over-boosted gamma-ray spectrum given the relatively large m W . -Generally speaking, for a given DM mass m χ 0 1 , larger today's DM annihilation cross sections σ v 0 correspond to smaller χ 2 GC E , i.e. better compatibility of the predicted gamma-ray spectrum with the observed GCE, except for the strong resonant region m χ 0 1 ≈ m h SM /2. Understandably, a sufficiently strong production rate of the gammaray flux from DM annihilations is helpful to act as an additional gamma-ray source to explain the GCE.
-Resonance enhancement near m χ 0 1 ≈ m h SM /2 alleviates the pressure on the NMSSM parameters of produc- and C A i bb in Eq. (8). Therefore this mass region can provide the NMSSM parameters more flexibility to fine tune the shape of the predicted gamma-ray spectrum to better fit to the GCE and result in smaller χ 2 GC E . -Generically, smaller DM masses m χ 0 1 tend to produce smaller χ 2 GC E , since the larger DM number density with smaller DM mass can provide higher probability for DM to find each other and annihilate, yielding stronger gamma-ray flux.
Note that the DM annihilation processes in Eq. (7) can also contribute to the DM annihilations in the early Universe which make DM freeze out from the SM thermal plasma. In Fig. 2 we project the GCE-samples on the plane of m χ 0 1 versus m A 1 with χ 2 GC E represented by colors, where a solid line indicating m A 1 /m χ 0 1 = 2 is provided. We can see that 2m χ 0 1 /m A 1 is near the resonant region for most samples, implying that the light CP-odd Higgs A 1 also play an important role in the s-channel annihilation in the early Universe, except in some other strong resonant scenarios such as m χ 0 We note that s-channel DM annihilations can also proceed through Z boson and the CP-even Higgs bosons H i . However, we checked that their contributions are small because the coupling C Z χ 0 1 χ 0 1 is week for a Singlino-like χ 0 1 in the DM annihilations, and the mass relation 2m χ 0 is mostly off-resonant for the H 1 -mediation.

GCE confronting DM direct detection in NMSSM
In Fig. 3 we plot all GCE-samples on the plain of R C A 1 bb versus Abs |N 13 | 2 − |N 14 | 2 with color denoting the value We note that other refined scan strategy may explore more difficult corners of NMSSM parameter space providing finely tuned resonant mechanisms, which may reconcile the GCE interpretation and the DM direct detection results. However, we believe that the approximated correlation in Eq. (20) and the GCE-samples we obtain in Fig. 3 and in Fig. 4 can serve as an illustration of the tension between these two DM phenomenologies in the Z 3 -NMSSM when the singlet superfield S play an important role in the DM and Higgs sector. The future update of PandaX-nT and XENONnT experiments will continue to challenge the above scenario.

Summary
In this work we consider an interesting scenario in the Z 3 -NMSSM where the singlet S and Singlino S 0 components play important roles in the Higgs and DM sector. Guided by the analytical argument, we perform a sophisticated scan by considering various observables such as the SM Higgs data and the B-physics observables . We first collect samples without applying the strict DM direct detection (DD) bounds and analyze their features about the GCE interpretation. We find that χ 0 1 DM are mostly Singlino-like and annihilation products are mostly the bottom quark pairsbb through a light singlet-like CP-odd Higgs A 1 . Moreover, a good fit to the GCE spectrum generically requires sizable DM annihilation rates σ bb v 0 in today's Universe. However, the correlation between the coupling C A 1 bb in σ bb v 0 and the coupling C Z χ 0 1 χ 0 1 in DM-neutron Spin Dependent (SD) scattering rate σ S D χ 0 1 −N makes all samples we obtain for GCE explanation excluded by the PandaX-II results. Although the DM resonant annihilation scenarios may be beyond the reach of our analytical approximations and scan strategy, the aforementioned correlation can be a reasonable motivation for future experiments such as PandaX-nT to further test the NMSSM interpretation of GCE.