Search for Scalar Dark Matter via Pseudoscalar Portal Interactions: In Light of the Galactic Center Gamma-Ray Excess

In light of the observed Galactic center gamma-ray excess, we investigate a simplified model, for which the scalar dark matter interacts with quarks through a pseudoscalar mediator. The viable regions of the parameter space, that can also account for the relic density and evade the current searches, are identified, if the low-velocity dark matter annihilates through an $s$-channel off-shell mediator mostly into $\bar{b} b$, and/or annihilates directly into two ${\it hidden}$ on-shell mediators, which subsequently decay into the quark pairs. These two kinds of annihilations are $s$ wave. The projected monojet limit set by the high luminosity LHC sensitivity could constrain the favored parameter space, where the mediator's mass is larger than the dark matter mass by a factor of 2. We show that the projected sensitivity of 15-year Fermi-LAT observations of dwarf spheroidal galaxies can provide a stringent constraint on the most parameter space allowed in this model. If the on-shell mediator channel contributes to the dark matter annihilation cross sections over 50$\%$, this model with a lighter mediator can be probed in the projected PICO-500L experiment.


I. INTRODUCTION
The standard model (SM) so far is quite successfully tested in the current high energy physics experiments. Nevertheless, the existence of the dark matter (DM) is indicated by various astrophysical measurements and astronomical observations. The Galactic center (GC) is an excellent place to generate the DM signals, because it concentrates a large quantity of dark matter. From the data collected by Fermi Large Area Telescope (Femi-LAT), several studies have found an excess of GeV gamma-rays near the region of the Galactic center [1][2][3][4][5][6][7][8][9][10][11], where the spectrum and morphology can be interpreted as the signals generated by annihilating dark matter (DM) particles [1][2][3][4][5][6][7][8][9][10]. The interpretation is not conclusive yet. A population of millisecond pulsars (MSPs) has been proposed as a plausible origin of the GC gamma-ray excess [12][13][14][15]. However, if so, the same region should contain a much more population of low-mass X-ray binaries than that observed so far [16]. On the other hand, Hooper et al. have also argued that if these MSPs convert more than a few percent of spin-down power into very high-energy e + e − pairs, then the inverse-Compton emission would exceed the observation by HESS [17].
Although several interaction types of DM models could be responsible for the GC gammaray excess, the predictions for the DM mass and relevant parameters might be in tension with the stringent constraints from direct detection experiments and colliders. A so-called "coy dark matter" model recently stressed that, for the DM fermions interacting with SM particles via a pseudoscalar mediator , the DM annihilation cross section into b quarks can be large enough to provide a good fit to the GC gamma-ray excess, while its corresponding t-channel process (the DM-nucleus scattering), relevant to the direct detection, is suppressed by four powers of momentum transfer, and, on the other hand, only a limited portion of the allowed parameter region can be constrained at the 14 TeV LHC run [18,19].
Another idea, similar to the secluded dark matter scenario [42], is to introduce a model with hidden sector mediators, in which the DM first annihilates into on-shell mediators, and subsequently the mediators decay to the SM particles via a very small coupling [29,[43][44][45]. Because the low-velocity annihilation cross section is highly insensitive to the mediator's coupling to the SM particles, it can thus explain the GC gamma-ray excess and easily evade the stringent constraints from the direct detection and collider experiments.
Motivated by these results as mentioned above, we consider the scalar DM interactions with SM quarks via a pseudoscalar mediator. A pseudoscalar particle is interesting from a model-building point of view, because a model with an extension of the Higgs sector, e.g. a two-Higgs doublet model, can naturally contain a such state. For this model, interactions through the s-channel exchange of a pseudoscalar with SM quarks could account for the GC gamma-ray excess [46,47], and constraints from current direct detection experiments and the LHC could be obviated. We find that, if the annihilation is only given by the pseudoscalar mediated s-channel process, the parameter region m A 2m φ , where m A and m φ are the masses of the mediator and dark matter, respectively, is ruled out by a combination analysis of various experiments, especially by the relic density constraint and observations of dwarf spheroidal galaxies (dSphs) [48]. The viable regions of the parameter space will be discussed.
On the other hand, because not only the scalar DM s-channel annihilation but also the DM annihilation into two hidden on-shell pseudoscalar mediators, via t and u channels, are s-wave processes, we find that, for m A < m φ , broad parameter regions that can provide a good fit to the observed GC gamma-ray excess and evade the current searches. Unlike the previous works, where the hidden sector mediator model is characterized by very small values of the mediator-SM couplings, the strong suppression of signals for direct detections and colliders in the present hidden sector is mainly due to the coupling structure for which the mediator interacts with the SM quarks via pseudoscalar couplings. We will show that, for m A < m φ , the DM annihilation into on-shell mediators gives sizable contributions to the observed GC gamma-ray excess.
In this paper, we adopt the framework of the simplified model, where the scalar dark matter annihilates through a spin-0 mediator with pseudoscalar couplings to SM quarks, which are assumed to be proportional to the Yukawa couplings motivated by minimal flavor violation [49]. Using a minimal set of parameters, the simplified model can not only capture the feature of a specific ultraviolet (UV) complete model, but also provide a generic framework to perform the experimental data analysis.
We further consider the updated and projected bounds set by the gamma-ray observations of dSphs [48], direct detection experiments, and LHC monojet result [50]. The dSphs, containing little dust or gas, are believed to be DM dominated. So far, no gamma-ray emission has been measured from dSphs by Fermi-LAT. For the direct detections, because the DM-nucleus interaction in this model is spin dependent, the LUX [51] signals mostly arise from the unpaired neutrons inside the abundant Xe isotopes, which are the LUX detector materials, while the PICO-60 [52][53][54][55], using the CF 3 I and C 3 F 8 as targets, detects the signals mainly via unpaired protons inside the target nuclei. We find that in the direct searches the PICO results set the most stringent bound which is also insensitive to the choice of the parameter set. At the LHC, the production of the scalar DM pair via the spin-0 mediator in the monojet accompanied by the missing transverse energy ( E T ) is dominated by the gluon fusion processes. We give the monojet constraint on the relevant parameter space, by taking into account the recent CMS 13 TeV results with 12.9 fb −1 and projected sensitivity for the high luminosity LHC.
The layout of this paper is as follows. In Sec. II, we introduce the simplified scalar DM model with a spin-0 mediator that couples to SM quarks via pseudoscalar couplings. This section includes the formulas about decay widths of the mediator, that are relevant to calculations for the DM relic density, indirect detection, and monojet result at the LHC. In Sec. III, we describe the approaches in detail for model constraints due to observations, and their implementations. For this model, the constraints on the parameter space are presented in Sec. IV. The conclusions are summarized in Sec. V. The brief descriptions for the relic abundance and results of thermally averaged annihilation cross sections are given in Appendices A and B, respectively.

II. THE SCALAR DARK MATTER MODEL
We consider a simplified model, where the scalar dark matter (φ ( * ) ) annihilates through a spin-0 mediator (A) with pseudoscalar couplings to SM quarks (q). The effective Lagrangian is given by for the complex scalar DM, or for the real scalar DM, where the latter one with the factor of 1/2 gives the identical expression for the direct detection rate and annihilation cross section as the former case. Motivated by the minimal flavor violation ansatz [49], we assume the pseudoscalar-SM quark couplings are proportional to the associated Yukawa couplings, given by g q = g √ 2m q /v with m q being quark's mass and v = 246 GeV. For simplicity, in the following, we consider the DM particle to be a real scalar, while the generalization to a complex scalar one is straightforward. Because this work is motivated by the fermionic case, we will assume the mediator is a pseudoscalar, such that the A-φ ( * ) -φ coupling is CP violating; however, our conclusions are independent of the mediator's parity. For a pseudoscalar mediator, the UV completion of the simplified model could be built up to relate to the Higgs-extension portal [28,34,35,[37][38][39][40][41]56] or axion-portal DM models [57], which might contain more (mediator) particles in addition to an extra pseudoscalar compared with the SM and thus have richer phenomenologies (see also discussions in Refs. [41,56] and LHC constraints in Refs. [37][38][39][40]).
We have collected the results for the annihilation cross sections, that are relevant to the relic density and indirect detection searches, in Appendix B. The partial decay widths of the pseudoscalar particle are where The DM particles annihilating into the quark pair, via the s channel, is s wave, while if kinematically accessible (m φ > m A ), the annihilation, via t and u channels, into on-shell mediators, is also s wave. We will show that the latter cannot be neglected, when it is allowed. Motivated by the phenomenological considerations, we choose the following three scenarios to explain the GC gamma-ray excess: (i) scenario s1: completely due to the DM annihilations into quark pairs, i.e., q σv qq , (ii) scenario s2: due to the annihilations equally into quark pairs and into on-shell mediators, i.e., q σv qq = σv AA , and (iii) scenario s3: mainly due to the annihilation into on-shell mediators, assuming that σv AA = 20 q σv qq .

III. DESCRIPTION FOR EXPERIMENTAL CONSTRAINTS
A. Indirect searches

Galactic Center Gamma-Ray Excess
The differential flux of a gamma ray from a given angular region ∆Ω, originating from the annihilation of scalar DM particles, is where η ≡ "2" is for the self-conjugated DM (e.g. real scalar DM) and "4" for non-self-conjugated DM (e.g. complex scalar DM), dN f γ /dE is the energy spectrum of DM prompt gamma rays produced per annihilation into the final state f, and the scalar DM mass is denoted by m φ . The factorJ Ω is an average of the J factor over the solid angle ∆Ω, covering the region of interest (ROI), given byJ where the integral is performed along the line of sight (l.o.s.), ρ(r) is the DM halo profile with r being the distance from the Galactic center, and ψ is the angle observed away from the Galactic center and dΩ ≡ cos b d db satisfying cos ψ = cos b · cos with and b being the longitude and latitude in the Galactic polar coordinate, respectively.J c Ω is the canonical value ofJ Ω , while J γ parametrizes the deviation from the canonical profile due to variation of the DM distribution slope γ. The values ofJ c Ω and J , sensitive to astrophysical uncertainties, depend on the observational ROI in a particular analysis. Following [8,9] where the tail of the spectrum has extended to higher energy, we employ the ROI of | | < 20 • and 2 • < |b| < 20 • (i.e., 40 • × 40 • square centered on the GC with the latitude |b| > 2 • ) to study GC gamma-ray emission.
We adopt the Galactic DM density distribution described by a generalized Navarro-Frenk-White (gNFW) halo profile [58,59], In the numerical analysis of φφ → A * → qq, we use the two-body spectra dN q γ /dE from the PPPC4DMID results, which, generated using PYTHIA 8.1 [60], have included the electroweak corrections [61,62]. When the DM annihilation into two on-shell mediators occurs, the process has two final states, i.e., four quarks produced, and the photon spectrum dN q γ /dE defined in the DM center of mass frame can be written in terms of the spectrum (dN q γ /dE ) A described in the rest frame of the mediator by considering a Lorentz boost [63], where t max = min 1, We fit the DM mass and its corresponding annihilation cross section σv to the prompt gamma energy spectrum of the GC excess extracted by Calore, Cholis, and Weniger (CCW) [8]. CCW studied Fermi-LAT data covering the energy range 300 MeV−500 GeV in the inner Galaxy, where the ROI extended to a 40 • × 40 • square region around the Galactic center with |b| ≤ 2 • masked out. We perform the goodness of fit with a χ 2 test statistic for the total annihilation cross section σv and m φ , where a total of 24 bins are used in the energy range 300 MeV− 500 GeV, dΦ γ /dE i and (dΦ γ /dE i ) obs stand for the model-predicted and observed GCE flux in the ith energy bin, respectively. Here the covariance matrix Σ contains statistical error, empirical model systematics and residual systematics, where the latter two are correlated across different energy bins. Way dSphs based on 6 years of data [48,64]. The Fermi-LAT data contain 24 bins, and the bin energy range spans from 500 MeV to 500 GeV, and no significant gamma-ray excess was measured from the dSphs. The observation can offer stringent constraints on the annihilation cross section of DM particles. We will use combined results of the 15 dSphs, recently reported by Fermi-LAT [48], in the theoretical analysis. Using the joint likelihood method, we compute and then add the bin-by-bin delta-log-likelihood for the 15 dSphs. The profile likelihood ratio test statistic (TS), following a χ 2 distribution, is given by with the profile likelihood for target (each dSph) k described as where N dSph = 15 and N bin = 24 are the numbers of dSphs and bins, respectively. Here the binned likelihood, L ki , is a function of the gamma-ray's energy flux within the "ith" bin 1 , E i ∆Φ γ , and the J factor likelihood for a target k is modeled by a normal distribution [65], For the J factor of a target k, J k is its expected value, and J o,k is the measured nominal value of an error σ k . For a given m φ , σv andJ k are the respective best-fit values for the DM annihilation cross section and the J factor, corresponding to the minimum value of −2 k=N dSph k=1 ln L k , whileĴ k are the conditional maximum likelihood estimators (MLEs) of the nuisance parameters when σv is fixed to a given value. σv 0 is the upper limit of the cross section, corresponding to the null measurement, and its 95% confidence level (C.L.) limit can be obtained by increasing the value of σv from σv until TS = 2.71.
We will use the results for J o,k and σ k given in Ref. [48]. These values were obtained assuming an NFW halo profile and shown to be insensitive to the models of dark matter density profile if the central value of the slope is less than 1.2. On the other hand, we also consistently use the PPPC4DMID results to generate the relevant two-body and four-body dN f γ /dE spectra as the studies of the GC gamma-ray excess.
B. Direct detections

The effective Lagrangian at the nucleon level
To obtain the differential DM-nucleus scattering rate at direct detection experiments, we rewrite the pseudoscalar-quark interacting Lagrangian at the nucleon level with the replacement, where the pseudoscalar coupling with the nucleon (labeled by p ≡ proton and n ≡ neutron) can be expressed in terms of quark spin contents of the nucleon [66], ∆q (N ) , and is given by The ratio c p /c n is sensitive not only to m u /m d , but also to the values of ∆q (N ) 's [67]. For illustration, numerically we will adopt the following two sets of parameters [66]: Meanwhile, we will use m u /m d = 0.48 for set 1 and m u /m d = 0.59 for set 2, where the ratios are consistent with m u /m d = 0.48 ± 0.10 given in PDG [68]. In addition to that, all quark masses are also used from PDG and consistently rescaled to µ = 1 GeV in the direct search studies.

The nuclear recoil rate and DM velocity distribution function
In direct detection, the scattering rate of the DM particles off target nuclei is given by where N T is the atomic number of the target, f ⊕ ( v, t) is the DM velocity distribution in the Earth frame, and v min is the minimal DM velocity needed for a nucleus to scatter with a recoil energy E R . Throughout this paper, I will consistently use the local DM density ρ 0.4 GeV/cm 3 . Here, where v ⊕ is the relative velocity of the Earth in this (Galactic) frame, and its magnitude can be with v 232 km/s being due to the motion of the Sun relative to the Galactic frame, u E 30 km/s being the relative speed between the Earth and Sun, and γ 60 • being the angle between the Milky Way's disk and Solar System [69][70][71].
The gNFW halo profile given in Eq. (9) exhibits a double-power law density; the inner log slope of the halo density near the core is −γ, while the log slope at large radii is −3. However, the isotropic Maxwellian velocity distribution, which is usually used in the direct detection analysis, arises from the density slope of −2. To have a velocity distribution function consistent with the gNFW halo profile given in Eq. (9), we adopt the isotropic velocity distribution ansatz, which can reproduce the Eddington formula with double power-law density, given by [72], where v 0 is the dispersion, v esc is the escape velocity, and k 2, the best fit to the gNFW profile model, is controlled by the outer slope of the halo density. We use v 0 = 220 km/s and v esc = 544 km/s [73].

The differential rate
At the leading order, the relevant nucleon matrix element is related to the following nonrelativistic operator: where the momentum transfer is q = p − p, S N is the nucleon spin, and O N is given by for the complex scalar DM case, The differential rate can be expressed as where with j being the nuclear spin, "T " denoting the target nucleus, and N + and N − being nucleon's nonrelativistic fields involving only creation and annihilation operators, respectively. The nuclear form factors F N,N Σ for various nuclei, given in Refs. [74,75], are functions of The nuclear shell model calculation showed that the expectation values of the nucleon S N and the spin of the initial target nucleus J are mainly due to the unpaired nucleon [76].

Null results in direct-detection experiments: LUX and PICO
To determine the stringent exclusion bounds on physical parameters due to the null results obtained in direct detection experiments, we use a Poisson likelihood function to model the distribution of the observed events and adopt the likelihood ratio test statistic [68,[77][78][79], where the likelihood of data is given by the product of Poisson distributions, the MLE of n, such that 0 ≤ λ(n) ≤ 1 for any n th i ≥ 0. Here n i = n th i + n b i with n th i and n b i being the event numbers for the theoretical prediction and expected background, respectively. The TS of goodness of fit has an asymptotical χ 2 distribution and can be rewritten as where the last term in Eq. (32) is zero when n obs i = 0. We take Λ ≡ m A /(g φ g) as the relevant parameter. Thus, adopting TS = 2.71 yields a one-side 95% C.L. upper limit for Λ with respect to any given m φ . In general, for each bin (or each module) i, the number of events theoretically expected at a direct detection experiment can be expressed by where E i is the exposure of the experiment and i T (E R ) is the efficiency and acceptance that a nucleus T with recoil energy E R is detected. Considering the background comes with uncertainty in the form n b i =n b i ± σ b i , the likelihood function is modified as where the probability density function of n b i is modeled as a Gaussian distribution. In the following, we briefly describe the very recent LUX and PICO data that are relevant to the TS calculation.
The very recent LUX data, using a xenon target, released is based on a complete run of 3.35×10 4 kg·days exposure, called WS2014-16 [51]. The events passing the cut with distance to the wall larger than 4 cm are selected, but those with 3 cm< r < 4 cm are neglected [80], where the fiducial boundary is defined as 3 cm inwards from the observed wall in S2 space, and the radius is about 6−19 cm corresponding to the drifted electrons' drift time between 40 µs and 300 µs [80,81].
There are about 85% of events selected, if the number of events is roughly proportional to the fiducial volume. Assuming that the DM events distribute evenly below and above the red solid curve in Fig. 1 of Ref. [51], we restrict ourselves to the signal region only below this curve. We therefore multiply the efficiency by an additional factor 1/2 × 0.85. On the other hand, we quote the efficiency from Fig. 2 of Ref. [51] and take 3 phd ≤ S1 ≤ 50 phd (with phd ≡ photons detected), such that four events were observed compared with 3.3 background events predicted [80,81], where the latter due to leakage from the electron recoil band are assumed to be equally distributed in S1. The PICO Collaboration has also reported the result with a total exposure after cuts of 1167 kg·days at a thermodynamic threshold energy of 3.3 keV using the PICO-60 dark matter detector and the bubble chamber filled with C 3 F 8 [53]. The PICO-60 C 3 F 8 improves the constraints on the DM parameters, compared with PICO-2L run 2 experiment [54]. We use the best fit efficiency curves for C and F, as given by the solid lines in Fig. 4 of Ref. [55]. The background is predicted to be 0.25 ± 0.09 single bubble events from neutrons, 0.026 ± 0.007 events from electron recoils, and 0.055 ± 0.007 events from the coherent scattering of 8 B solar neutrinos. No single-scattering nuclear recoil candidates are observed.

C. Monojet searches
The studies for monojet plus missing transverse energy (MET) are one of the important searches for dark matter at the LHC. Here, we employ the very recent CMS 13 TeV results, corresponding to an integrated luminosity of 12.9 fb −1 [50]. Using a profile likelihood ratio, we employed the CL s method [82] to calculate the 95% confidence level (C.L.) upper limit on the E T signal events for pp → φφj at the reconstructed level, where the Standard Model background within a bin is modeled as a Gaussian distribution, but the correlations between uncertainties of the background yields across different E T bins are neglected.
To obtain the constraints of parameters in the simplified model from the upper limit of monojet signals involving the DM missing transverse energy at the trigger (reconstructed) level, we implement the present model into FeynRules [83,84] to get a UFO output [85], which is then used in MadGraph5 aMC@NLO [86] for the simulation of the relevant monojet events. We set the renormalization scale (µ R ) and factorization scale (µ F ) to be ξH T /2, where H T is the sum of the missing transverse energy and the transverse momentum of the jet (j), and ξ ∈ [1/2, 2] denotes the scale uncertainties.
We use MadAnalysis 5 to analyze the events of simulations [87]. The next-to-leading order (NLO) NNPDF3.0 [88] parton distribution functions (PDFs) with the corresponding α s (M Z ) = 0.118 are used to generate the signal events, which are further hadronized by using PYTHIA8 [60].
For jet clustering, consistent with the CMS study, we construct the (AK4) jets by employing the anti-k T algorithm [89] with the distance parameter R = 0.4, as implemented in FastJet [90]. The selection cuts for jets at the reconstructed level are p T > 20 GeV and |η| < 5.0, while the leading one in the event is required to have p T,j 1 > 100 GeV and |η j 1 | < 2.5. We also impose the jet veto that rejects events if the azimuthal separation between p miss T and the directions of each of the four highest p T jets with p T > 30 GeV is smaller than 0.5 radians. This criterion has been used by the CMS to suppress the QCD multijet background.
Because the triggers for events with E T > 300 GeV become full efficient at the CMS, we will take into account the missing transverse energy within the range E T = 350−590 GeV in the numerical analysis. Within these bin widths, the 95% C.L. upper limit for the total number of signal events due to the generation of the DM particles will be used in the analysis.

A. Indirect detections
We first fit the DM mass and its corresponding annihilation cross section σv at an average velocity of v ∼ 10 −3 c to the GC prompt gamma energy spectrum [8]. The results of fits are plotted in Fig. 1, and summarized in Table I, together with ±1σ errors, χ 2 min /dof , and p value. The results that we show correspond to ρ = 0.4 GeV/cm 3 and J γ = 1 which refers to the inner log slope to be −1.26(= −γ). It seems that the model results shown in the left panels of Fig. 1 do not appear to be a good fit to the prompt gamma spectrum. However, due to the strong correlations of the systematical errors among different energy bins, the best-fit values give good fits with p values > 0.35 2 .
The DM s-channel annihilation is dominated by thebb final state for m φ > m b , because the mediator interacts with quarks via the Yukawa-like coupling. If the DM annihilation into the onshell mediator pair is kinematically allowed, the fitted regions depend on the mass of the mediator, as shown in Fig. 1, where, for illustration, we have shown results for three values of m A /m φ = 0.2, 0.5, and 0.8 in scenarios s2 and s3.
If the low velocity DM annihilation cross section is dominated by σv AA (scenario s3), and the produced A particles are nonrelativistic, then two final states, i.e., four b quarks, are generated in the decays of two A particles, such that the best-fit GC excess results for the DM mass and annihilation cross section are therefore larger by a factor of ∼ 2 compared to that in the pure schannel annihilation case (scenario s1). The best-fit values of m φ and the cross section will further decrease either for a smaller m A due to the fact that the PPPC4DMID gamma-ray spectrum needs to be boosted from the A particle rest frame to the dark matter rest frame in the fit or for a larger contribution arising from the DM s-channel annihilation intobb (scenario s2).
For the scenario s2, the parameter region with m φ ∼ 50−60 GeV and m A ∼ 10−12 GeV yields a good fit, for which the gamma-ray spectrum produced from the on-shell A decaying into the final state bb can be negligible due to the smallness of the phase space, and therefore the result is dominated by the DM s-channel annihilation, via an off-shell A, into bb. On the other hand, for scenarios s2 and s3, the GC gamma-ray excess data can also be fitted by the parameter region m φ ∼ 30 GeV, where in addition to the contribution arising from the DM annihilation into two on-shell mediators, which dominantly decay into thecc pairs in the final state, s2 receives a sizable contribution from the DM s-channel annihilation intobb.
The results show that some GC gamma-ray excess allowed regions are excluded by the observations of dSphs. It should be noted that the J values of dSphs are obtained subject to the assumptions that the dSphs are spherically symmetric and have negligible binary motions [65].
Note also that the uncertainty due to the DM profile of the Galactic center is not included in the fit of the GC gamma-ray excess; the allowed region of the annihilation cross section shown in  The corresponding p value of χ 2 min is given. σv = q σv φφ→qq + σv φφ→AA is the total cross section, using ρ = 0.4 GeV/cm 3  square root of the number of observed dSphs [92]. Assuming that the Fermi-LAT Collaboration can successfully collect 15-yr gamma-ray emission data about 60 dSphs [93], the σv limits will be further improved by a factor of ∼ 3.16 in the future, as shown by the dot-dashed curves in Fig. 1, and this model is very likely to be testable.

B. Direct detections
In Fig. 2  [94] 2σ and 3σ allowed regions, and exclusion results extracted from earlier measurements by PI-CASSO [95], COUPP [96], XENON100 [97], and SuperCDMS [98], where the method for treating     these null data can be referred to Ref. [67]. The approach of the DAMA modulation analysis is similar to that given in Ref. [67]. For the DAMA signals, two regions at m φ ∼ 10 GeV and at m φ ∼ 40 GeV can be interpreted if the DM particle scatters on the sodium for the former region and iodine for the latter. However, the PICO measurements seem to strongly disfavor the parameter space fitted from the DAMA modulation data.
Our results are summarized as follows. Although the PandaX-II [99] and XENON1T [100], using also the xenon target, have recently obtained a slightly stronger bound on the (spin-independent) cross section by a factor of 2.5 compared with LUX WS2014-16, the PICO-60 measurements still give the most stringent exclusion bound, which is insensitive to the choice of the parameter set, among the current direct detection experiments 3 . In the following, we will therefore use the PICO-60 results on the model analysis.

C. Monojet and scenario s1
The monojet+MET search can provide a relatively stronger constraint on parameters when the mediator is produced on shell at the colliders. For this case, the monojet cross section can be approximated as σ(pp → jφφ) ∼ σ(pp → j + A) × BR(A → φφ). Therefore, for m A > 2m φ , the monojet search may constrain the GC excess region where only the DM s-channel annihilation into the SM quark pair is relevant; in other words, the monojet constraint on the GC excess region can be categorized to the scenario s1. For the parameter range with m A < 2m t , the total width of the mediator (Γ A ) obtained using Eqs. (3), (4), and (5), is always small; for instance we have Γ A /m A < 0.02 for g φ , g < 5, consistent with the narrow width approximation.
This monojet constraint depends on the value of g. For illustration, take m A > 2m φ with m A = 100 GeV and m φ = 46.4 GeV as an example. For g 1.3, the present monojet result cannot provide a sufficient constraint on g φ in the range of g φ ∈ [0, 4π], because BR(A → φφ) is already larger than ∼ 90% for g φ = 1. On the other hand, for a large g limit, the monojet cross is independent of g, such that, for g 3, we get that g φ 0.36 is excluded at 95% C.L. by the CMS monojet search.
Numerically, we find that g = 2, which will be used in the analysis, can provide a stronger limit on the value of gg φ .
In Fig. 3, taking an illustrative DM mass of m φ = 46.4 GeV, we show results in the (m A , g φ g) plane, where the GC gamma-ray excess allowed region at the 3σ C.L. is given for the scenario s1.
There, adopting g = 2, we show the (hatched magenta) region excluded by the very recent CMS monojet search with 12.9 fb −1 of data at 13 TeV, and the projected limit (magenta dashed curve) for the high luminosity LHC (HL-LHC) with an integrated luminosity of 3000 fb −1 at 14 TeV [103]. The current CMS monojet constraint is shown to be much less restrictive for this model.
Because we have assumed that the mediator-quark couplings are proportional to the quark's mass, the production of A, mainly via the top loop, is dominated by gluon fusion at the LHC. Therefore, we may expect that the S/B (number of signal to number of background ratio) is approximately the same as the present value, such that the projected limit on g φ g corresponds to an improvement The dashed lines with magenta, red, orange, and green colors are the 95% C.L. upper limits due to projected sensitivities for HL-LHC monojet, dSphs, PICO-500L C 3 F 8 , and PICO-500L CF 3 I. In general, the thermal relic density can be accounted for by the gray regions. In (a), (b), and (c), the monojet constraints and thermal relic density shown in the black, brown, and cyan solid curves correspond to the chosen values: g = 4π, g = 2, and g φ = 4π, respectively. by a factor ∼ 4.3. As shown in Fig. 3, the projected HL-LHC limit could constrain the favored parameter region. However, the projected limit might become less restrictive if g is much different from 2.

Scenario s1
In the scenario s1, we consider that the DM particles annihilate only via an s-channel mediator to the SM quark pair at a velocity v ∼ 10 −3 c. In Fig. 3, in addition to the parameter region favored by the GC gamma-ray excess and the 95% C.L. upper limits placed by LHC monojet searches, as discussed previously, we also show regions excluded from the PICO-60 CF 3 I measurements (hatched green) and Fermi-LAT dSphs observations (hatched red). The PICO Collaboration has proposed a ton-scale PICO-500L detector, having an active volume of about 800 L [104,105]. The predicted bound from the projected sensitivity of the PICO-500L, assuming an exposure of 500 kg·yrs and the same detection efficiency used in PICO-60, is plotted in the dashed orange or dashed green line in Fig. 3, if the bubble chamber is filled with C 3 F 8 or CF 3 I [106]. The projected PICO experiments will constrained the region with m A 25 GeV. However, this region is already excluded by the combined analyses of the relic abundance and the observed dSphs (see Fig. 3 and also discussions below).
Note that the DM annihilation cross section obtained in the GC excess could be revised by a factor ∈ [0. 33, 5.88] due to uncertainties of the DM profile near the Galactic center and local DM density; including these uncertainties, the new GC excess boundaries are denoted by the thick dashed blue lines in Fig. 3.
As for the relic density constraint, Ω DM h 2 = 0.1198 ± 0.0026 [68,107], for comparison, we show all allowed parameter regions in Fig. 3, considering the parameter range of g φ , g ∈ [0, 4π]. If m A > m φ and the temperature at the decoupling time is not high enough to produce the on-shell A, i.e., m A > m φ + K φ (with K φ the thermally kinetic energy of φ), then the DM annihilation is only relevant to the s channel φφ →qq, dominated by the bb pair in the final state. We find that only the parameter region with 2m φ < m A 2.7m φ is compatible with the combination of all data as well as the GC gamma-ray excess within errors; near the resonance region, if using the canonical astrophysical parameterJ c Ω as input, we need to adopt a much larger value of g φ ( 10) to have a large width of the mediator and then to suppress the resonant enhancement on the annihilation cross section, so that we can have a good fit to the combination of the relic density and GC excess, otherwise, we need a large revision to the adopted astrophysical parameters for accounting for a smaller coupling g φ . The projected sensitivity of dSphs can strongly constrain this allowed region, as shown in Fig. 3.
For m A < m φ +K φ , the s channel φφ →qq and (t,u) channels φφ → AA are relevant to the relic abundance (as well as the GC gamma-ray excess), because these two annihilation processes are s wave. To investigate the channel dependence of the relic abundance, we display results in Fig. 3 for the three coupling values of (i) g = 4π, (ii) g = 2, and (iii) g φ = 4π, respectively. Only the case (i), dominated by φφ →qq, is consistent with the GC excess allowed region under the requirement of the scenario s1, but, however, is completely excluded by the gamma-ray measurement from dSphs.
Although cases (ii) and (iii) seem to be reconciled with the observed dSphs and GC gamma-ray excess in the (m A , g φ g) plane, they however contain a sizable contribution from the DM annihilation into on-shell mediators. This contribution is comparable with that from the DM s-channel annihilation into the bb pair, in contrast with the assumption of the scenario s1. This may imply that for m A < m φ , the DM annihilation into on-shell mediators plays an important role in the phenomenology of the GC gamma-ray excess.
In concluding this subsection, it is interesting to note that the constraints arising fromtt + E T (with A → E T ) channel and the mediator's visible decay channels at the LHC: pp → A → τ + τ − , γγ, could be comparable with and/or complementary to the monojet result. 4 This was recently stressed by Banerjee et al. [33].

Scenarios s2 and s3
In Figs. 4 and 5, we, respectively, consider two alternative scenarios, s2 and s3, for which, when DM particles move at an average velocity v ∼ 10 −3 c, the former is described by σv AA = q σv qq for the DM annihilation around the GC, and the latter is assumed to respect σv AA = 20 q σv qq .
The constraints due to various experiments are given in the (m φ , g φ ) plane, where g φ is the only coupling relevant to the DM annihilation into the on-shell mediator pair. We also display results for the three values of m A /m φ = 0.2, 0.5, 0.8 to illustrate the mass dependence of the mediator.
In the scenario s1, the value of the low-velocity σv b b is consistent with the thermally averaged cross section required by the relic density (∼ 1.78 × 10 −26 cm 3 /s corresponding to m φ /T f 20 with T f the freeze-out temperature for the real scalar DM), whereas, when the φφ → AA channel is open and m A > 2m b , the fitted cross section becomes larger for the GC gamma-ray excess as shown in Figs. 1, 4, and 5. For the latter, at the face value, one may worry the resulting relic density is too small. However, we find that the thermally averaged annihilation cross sections for σv AA and q σv qq at the freeze-out temperature are, respectively, smaller by a factor of ∼ 0.77 and ∼ 0.75, compared to the corresponding low-velocity annihilation ones. Thus, we can obtain the parameters that produce the correct relic abundance and also provide a good fit to the GC gamma-ray excess. The relevant formulas are collected in Appendices A and B.
We observe that the allowed parameter regions, constrained by the current measurements, are, respectively, in the ranges of g φ ∈ [0. Although, the present PICO-60 results cannot provide sufficient constraints on the parameter space, the PICO group is planning to run PICO-500L using C 3 F 8 as the target material, for which CF 3 I can be a substitution [104][105][106]. The constraints from PICO-500L, assumed to have a run of 500 kg·yr exposure, can be considerably more restrictive for a light mediator with m A 0.5m φ in scenario s2 and m A 0.2m φ in the scenario s3, if, in particular, the detect chamber is filled with CF 3 I (dashed green line).

V. CONCLUSIONS
In light of observations of the GC gamma-ray excess, we have investigated a simplified model, in which the scalar dark matter interacts with SM quarks through a pseudoscalar mediator, for which the s-wave DM annihilation can occur through an s-channel pseudoscalar exchange into the quark pair, or directly into two on-shell mediators if kinematically allowed.
If the contribution of the low-velocity DM annihilation is mainly due to the interaction through an s-channel off-shell mediator, we have found that only the parameter region with 2m φ < m A 2.7m φ is allowed by the combination of all data. For the allowed parameter space near the resonance region, if using the canonical astrophysical parameterJ c Ω as input, we need to have a much larger width of the mediator, corresponding to g φ 10, to suppress the so-called resonant enhancement on the annihilation cross section, such that a consistently good fit to the GC excess and relic abundance can become likely; otherwise, the canonical astrophysical parameters need be revised largely for accounting for a smaller coupling g φ .
For the region 2m b < m A < m φ (scenarios s2 and s3), although a larger low-velocity annihila- tion cross section is obtained to fit the GC gamma-ray excess, the thermally averaged annihilation cross sections at the freeze-out temperature are smaller by a factor of about 0.75−0.77 compared to the low-velocity annihilation ones. As a result, we find that the DM annihilation into two hidden on-shell mediators, which may be accompanied by an s-channel annihilation into the bb pair via an off-shell mediator, can be capable of accounting for the GC gamma-ray excess and relic abundance, and evade the current constraints from direct detections, observations of dSphs, and monojets results at the LHC. In this model, the signal suppression of the hidden sector is mainly due to the coupling of the pseudoscalar mediator to the SM quarks.
The current constraint from the CMS monojet plus missing transverse energy search are shown to be very weak for this model. The projected sensitivity of the monojet search at the high For a temperature T 3m φ , the thermally averaged annihilation cross section is given by [108], where K 1,2 are the modified Bessel functions and s is the center-of-mass energy squared. For the DM particles satisfying the condition x(≡ m φ /T ) 1, the annihilation cross section can be further approximated as where = (s − 4m 2 φ )/(4m 2 φ ) v 2 lab /4, and the cross sections, if kinematically allowed, are given by with n c = 3 being the number of the quark's colors. The former cross section is the s-channel process, while the latter contains the u and t channels. For the A's decay width, if kinematically allowed, we will consider the main channels, Γ A = q Γ(A →qq)+Γ(A → gg)+Γ(A → φφ), where the partial decay widths are explicitly listed in Eqs. (3), (4), and (5). As for m A ∼ Λ QCD , although its decay width depends on the channels of hadronization [26], it was found that the hadron decay widths can be neglected in the numerical calculation of the annihilation cross section [67].
In analogy to the relic abundance, for the indirect search case, the dark matter particles can be assumed to be at rest as a whole in the Galactic frame, and x in Eq. (B2) equals to 2/v 2 p , with v p the most probable speed of the dark matter distribution. Note that both σv Møl φφ→qq and σv Møl φφ→AA are s-wave dominant.