A global fit of the $\gamma$-ray galactic center excess within the scalar singlet Higgs portal model

We analyse the excess in the $\gamma$-ray emission from the center of our galaxy observed by Fermi-LAT in terms of dark matter annihilation within the scalar Higgs portal model. In particular, we include the astrophysical uncertainties from the dark matter distribution and allow for unspecified additional dark matter components. We demonstrate through a detailed numerical fit that the strength and shape of the $\gamma$-ray spectrum can indeed be described by the model in various regions of dark matter masses and couplings. Constraints from invisible Higgs decays, direct dark matter searches, indirect searches in dwarf galaxies and for $\gamma$-ray lines, and constraints from the dark matter relic density reduce the parameter space to dark matter masses near the Higgs resonance. We find two viable regions: one where the Higgs-dark matter coupling is of ${\cal O}(10^{-2})$, and an additional dark matter component beyond the scalar WIMP of our model is preferred, and one region where the Higgs-dark matter coupling may be significantly smaller, but where the scalar WIMP constitutes a significant fraction or even all of dark matter. Both viable regions are hard to probe in future direct detection and collider experiments.


Introduction
Weakly interacting massive particles (WIMPs) are promising candidates for dark matter (DM), and can be searched for at colliders and through direct and indirect detection experiments [1][2][3]. The scalar singlet Higgs portal model is among the simplest WIMP DM models. It comprises the Standard Model (SM) and a real singlet scalar DM field, S, which interacts with the SM Higgs field H through the operator S 2 H † H [4][5][6]. The scalar Higgs portal model can accommodate the DM relic density, would contribute to the invisible Higgs width, and it can be probed in direct and indirect DM searches.
As compared to previous Higgs portal model interpretations of the GCE [26][27][28][29][30], we provide a detailed numerical fit of the GCE signal within the scalar Higgs portal model, taking properly into account the theoretical uncertainty from the DM distribution. Furthermore, we allow for unspecified additional DM components beyond the scalar WIMP of our minimal model. We will show that, taking into account all the constraints, the scalar Higgs portal model can indeed describe the GCE signal, albeit only in a small region of parameter space near the Higgs resonance where the WIMP mass m S ∼ m h /2. The paper is organised as follows. In section 2 we introduce the scalar singlet Higgs portal model and briefly review previous collider and astroparticle analyses of this model. The DM annihilation γ-ray signatures of the scalar Higgs portal model are presented in section 3, together with a discussion of the galactic center excess signal and the astrophysical uncertainties due to the dark matter distribution. We present a detailed numerical fit of the strength and shape of the GCE γ-ray spectrum, including in particular the astrophysical uncertainties and allowing for unspecified additional DM components. Constraints on the model parameters from the Higgs invisible width, direct detection searches, independent searches for γ-rays and from the dark matter relic density are discussed in section 4. In section 5 we finally present a global fit of the GCE within the scalar Higgs portal model, taking into account the above-mentioned constraints. We conclude in section 6.

The scalar singlet Higgs portal model
The scalar singlet Higgs portal model [4][5][6] is among the simplest UV-complete WIMP DM models. The model comprises the Standard Model and a real scalar field, S, which is a singlet under all SM gauge groups. Imposing an additional Z 2 symmetry, S → −S, the scalar particle is stable and thus a WIMP DM candidate. The Lagrangian of the scalar Higgs portal model reads After electroweak symmetry breaking, the last three terms of the above Lagrangian become with H = (h + v, 0)/ √ 2 , v = 246 GeV, and where we introduced the physical mass of the singlet field, m 2 S = m 2 S,0 + λ HS v 2 /2. The scalar self coupling, λ S , is of importance for the stability of the electroweak vacuum and the perturbativity of the model, see e.g. [31], but does not affect DM phenomenology. 1 For the purpose of this paper, the model is thus fully specified by only two parameters beyond those of the SM: the mass of the scalar DM particle, m S , and the strength of the coupling between the DM and Higgs particles, λ HS .
The scalar singlet Higgs portal model defined in eq. (1) is certainly minimal, and possibly too simplistic. However, a coupling between a new gauge singlet sector and the SM through the Higgs bilinear H † H should be expected in a large class of SM extensions, as H † H is the only SM gauge singlet operator of mass dimension two. Even within the minimal scalar Higgs portal model, eq. (1), the S 2 H † H interaction term gives rise to a rich phenomenology, including invisible Higgs decays, h → SS, a DM-nucleon interaction through the exchange of a Higgs particle, and DM annihilation through s-channel Higgs, t-channel scalar exchange, and the S 2 h 2 interactions, see section 3.
The phenomenology of the singlet Higgs portal model has been extensively studied in the literature, see e.g. the recent reviews [29,33] and references therein. Other recent general analyses of the model have been presented in [34,35], while [36][37][38][39] have specifically explored the constraints from searches at the Large Hadron Collider (LHC). Astrophysical constraints, in particular from γ-lines, have been studied in [40][41][42][43]30]. Constraints on the scalar Higgs portal model from perturbativity and electroweak vacuum stability have been revisited in [44], while the possibility to drive inflation through a non-minimal coupling of the scalar to gravity has been analysed in [45] in light of current constraints. Extensions of the Higgs portal model that provide a similar phenomenology have been studied in [46][47][48][49].
3 The galactic center excess 3

.1 The Fermi-LAT observation
The presence of a GCE has been reported by several groups in the last few years [7][8][9][10][11][12][13][14][15]. The GCE seems compatible with a spherical morphology, extending up to at least 10 • away from the galactic center, and with a steep 'cuspy' radial profile [14,15]. The inferred energy spectrum is peaked at a few GeV in the usual E 2 × flux representation. Various astrophysical mechanisms and scenarios have been proposed to explain the excess [17,18,23]. On the other hand, intriguingly, it has been shown that the excess is also compatible with an interpretation in terms of DM annihilation, with a cross section close to the thermal value and with a DM mass around 50 GeV. Recently, the GCE has also been confirmed by the Fermi-LAT collaboration [16]. In the present analysis, we will use the results from [15], where a detailed spectral and morphological analysis of the excess has been performed and where the inferred energy spectrum has been made available together with an error covariance matrix. The covariance includes an estimate of systematic uncertainties related to the galactic foreground emission, inferred from a grid of different foreground models and from a scan of the typical model residuals along the galactic plane. Below m S = m h , only the s-channel Higgs diagram, figure 1 a), contributes to the annihilation cross section, and the relative strength of the different SM final states (f = t, h, Z, W, b, τ, c, g, γ) is determined by the SM Higgs branching ratios, independent of the Higgs-scalar coupling λ HS . Above the Higgs threshold, m S ≥ m h , all diagrams depicted in figure 1 contribute, and, in particular, the hh final state opens up. The strength of the annihilation into Higgs pairs, as compared to W, Z or top-quark pairs, depends on the size of the Higgs-scalar coupling λ HS .

Annihilation cross section and photon spectrum
We have implemented the scalar Higgs portal model into FeynRules [50] and used mi-crOMEGAs [51] linked to CalcHEP [52] to compute the velocity-averaged annihilation cross sections σv f for the various SM final states f . The loop-induced annihilation processes SS → gg, γγ have been included using the effective Lagrangian [53][54][55] where and Equation (4) takes into account the QCD corrections to the ggh vertex [54,55]. Here, v = 246 GeV is the vacuum expectation value of the Higgs, α s and α are the strong and electromagnetic couplings, respectively, s is the center-of-mass energy of the process and Q q is the electric charge of the quark q. We take into account the contribution from the bottom and top quarks in the sum, q = b, t. The strong coupling is evaluated at s, and we consider the one-loop running of α s . We checked the accuracy of our implementation by comparing our results for the Higgs branching ratios as a function of the Higgs mass to those of the LHC Higgs Cross Section Working Group [56]. We find agreement within a relative error below 5% for hgg and 15% for hγγ in the mass region between m h = 90 GeV (the minimal Higgs mass considered in [56]) and m h ≃ 300 GeV (above which the respective contributions to the annihilation are completely irrelevant). Note that differences are expected as the results from Ref. [56] include further higher-order QCD and electroweak effects not taken into account in our calculation. For m S > m h /2 the total Higgs width in our model is identical to the SM Higgs width and we take the theoretical prediction provided by the Higgs working group, Γ SM = 4.03 MeV [56]. For m S ≤ m h /2 we compute the invisible width Γ inv = Γ (h → SS) first and run micrOMEGAs with Γ h = Γ inv + Γ SM as an input parameter. We take the Higgs mass from the combined analysis of ATLAS and CMS [57], m h = 125.09 GeV. The relative contributions of the different SM final states to the annihilation cross section is displayed in figure 2 for two different choices of the Higgs-scalar coupling, λ HS = 1 and λ HS = 0.01, respectively. The fragmentation, hadronization and decay of the SM particles from the primary annihilation process, figure 1, produces a spectrum of γ-rays, predominantly from π 0 → γγ. The loop-induced process SS → h → γγ results in γ-ray lines at E γ ≈ m S , and is suppressed as compared to the continuum photon spectrum from pion decays. Searches for γ-ray lines will be discussed in section 4.3. In general, γ-rays can also be produced from electrons and positrons through inverse Compton scattering and synchrotron radiation. These contributions are important for DM particles annihilating predominantly into e + e − and µ + µ − final states [58]. In the scalar Higgs portal model, however, these annihilation channels are strongly suppressed, so that γ-rays from inverse Compton scattering and synchrotron radiation can be neglected in our analysis.
We have generated the γ-ray spectrum from SM particles, produced with a centre-of-mass energy E = 2m S , with the Pythia 8.209 event generator [59]. The contributions of 3-body final states from annihilation into off-shell gauge bosons, W W * and ZZ * , have been calculated using Madgraph5 aMC@NLO [60,61]. Comparing the spectra calculated with the default Pythia 8.209 generator to those obtained with Pythia 6 [62] we find differences of typically less than about 10%, c.f. [63,64]. Larger uncertainties occur for the gg final state, which is however less significant for the overall γ-ray flux. We have compared the γ-ray spectra for the  2-body final states to those presented in Ref. [63]. We find very good agreement in general, with some small deviations in the flux from annihilation into gg and cc final states. The spectra from the various final states, f = t, h, Z, W, b, τ, c, g, are combined, weighted by their relative strength as predicted within the scalar Higgs portal model as a function of the DM mass, m S , and the Higgs-scalar coupling, λ HS , see figure 2. The resulting γ-ray flux per unit solid angle at a photon energy E γ is where dN f /dE is the photon spectrum per annihilation for a given final state f , σv f is the corresponding velocity-averaged annihilation cross section, and ρ is the DM density. The integral has to be evaluated along the line-of-sight (l.o.s.) at an observational angle θ towards the galactic center. The l.o.s. integral of the DM density-squared over the solid angle dΩ is called the Jfactor, and is discussed in more detail in section 3.3.
The scalar Higgs portal model prediction for the γ-ray spectrum per annihilation, f dN f /dE × σv f / σv , is shown in figure 3 for different choices of λ HS and for dark matter masses m S that are of particular relevance in describing the GCE, as discussed below.

Dark matter density profile and uncertainties
The DM density in the Milky Way is only directly measured in the vicinity of the solar system, and only with a quite large uncertainty, mostly systematic in nature. In the inner galaxy no direct measurements are available since the gravitational potential is dominated by the baryonic matter. Extrapolations are thus necessary together with assumptions about the shape of the DM density profile, which is typically parameterized as a cored or cuspy profile. However, in the particular case of this analysis, where we are studying the DM interpretation of the GCE, we can limit the study to the DM profiles which are compatible with the measured shape of the GCE itself, i.e., cuspy profiles.
We will thus parameterize the DM profile as a generalized Navarro-Frenk-White (NFW) profile [65]: where r is the spherical distance from the galactic center, and we will assume γ = 1.2 ± 0.08 (Gaussian error) as given in [15]. For ρ s and r s , the scale density and scale radius of the profile, respectively, we will use the recent results from [66] where the authors study a canonical NFW profile (γ=1) using up-to-date measurements of the rotation velocity of the Milky Way as a function of the galacto-centric radius. In particular they provide contours in the ρ s -r s plane. From the contours it can be seen that the single parameters are not well constrained, but they are tightly correlated so that the contours can be simply approximated as a narrow line in the plane. We found that the power law ρ s = 42.7 GeV/cm 3 · (r s /kpc) −1.59 describes well the relation among these two parameters. Formally, the above relation is valid only for γ = 1 but we will use it for each γ within the explored uncertainty of ±0.08. This is expected to be a good approximation since different values of γ change the profile mainly in the inner kpcs from the galactic center, while the analysis of [66] is anyway performed for r 2 kpc, and is thus not very sensitive to moderate changes in γ. We also note that the analysis of [66] provides a value of the DM local density ρ ⊙ = 0.471 +0.048 −0.061 GeV cm −3 which is thus implicit also in our analysis. This value has a relatively small error, which is typical of analyses based on the Milky Way rotation curve (see also [67][68][69]), while purely local analyses provide larger errors ρ ⊙ ∼ 0.2 − 0.7 GeV cm −3 [70,71].
Given the above ρ s -r s relation and the further constraint γ = 1.2 ± 0.08, we determine the error on the GCE J-factor with a Monte Carlo procedure. The quantity of interest is the Jfactor integrated over the sky region analysed in [15], i.e., a 40 • × 40 • region centred on the galactic center and with a stripe of ±2 • masked along the galactic plane, In figure 4 we show the distribution of J 40 • that we obtain from sampling γ within its Gaussian uncertainty and with ρ s and r s uniformly distributed, taking into account their correlation. It can be seen that J 40 • is well approximated by a log-normal distribution with a width of σ logJ ≃ 0.43 (see figure 4). In the following we will use this J 40 • distribution to account for its uncertainty. J 40 • , nom is the nominal value of J 40 • for γ = 1.2 ρ s = 0.74 GeV/cm 3 , r s = 19.5 kpc, i.e., J 40 • , nom = 1.79 · 10 23 GeV 2 cm −5 . We also use r ⊙ = 8.0 kpc, although we verified that J 40 • varies very little varying r ⊙ in the range 7.5-8.5 kpc, which is the typical uncertainty on r ⊙ . Note that since the authors of [15] normalise the GCE flux dividing by the angular size of the analysed region, we also need to divide J 40 • by the corresponding solid angle ∆Ω = 0.43 sr. Figure 4: Probability density function of J 40 • , eq. (10). The red histogram shows the distribution obtained from sampling the parameter γ of the generalized NFW profile, eq. (9), within its Gaussian uncertainty. The additional NFW parameters ρs and rs are sampled uniformly, taking into account their correlation. The lognormal distribution fitted to the sampled distribution is shown in blue, with parameters µ and σ as specified in the figure.

WIMP contribution to dark matter
In this study we allow for the possibility that the dark sector is more complex than containing just one DM particle species. We could, for example, imagine a second non-WIMP DM component (such as axions or primordial black holes) which does not interact weakly with the SM and which does not annihilate into SM particles today. Hence, we consider the case that the WIMP DM density is a certain fraction, R ≤ 1, of the total (gravitationally interacting) DM: Here we assume that there is no difference in the clustering properties and hence in the density profiles of the WIMP and non-WIMP DM components. The annihilation signal today thus scales as φ ∝ R 2 . We will consider R as a free parameter in the fit of the GCE signal.

Fit to the GCE signal
In order to perform a fit to the GCE we use MultiNest [72,73], which allows to scan the parameter space under study much more efficiently than a simple random search. The respective parameters and scan ranges are summarized in table 1. The annihilation cross sections are computed using micrOMEGAs. The χ 2 for the GCE (including the contribution from J 40 • ) is computed as: where d i is the GCE measured flux in energy bin i from [15], t i is our model prediction, which depends on the parameters m S , λ HS , R and J 40 • , Σ ij is the covariance matrix given in [15], which includes statistical and systematic errors, and J 40 • , nom and σ logJ are as defined in section 3.3. The term δ ij (σ rel t i ) 2 represents a diagonal error equal to a fraction σ rel of the model prediction itself, which we add to the original Σ ij in order to take into account the model uncertainties in the annihilation spectrum. We choose σ rel = 10%, as discussed in section 3.2. Since we include in the error a dependence from the model itself, the likelihood which we use in the fit reads slightly differently from the χ 2 expression in eq. (12). Specifically, up to a constant factor which has no influence on the fit, the log-likelihood is where MultiNest is particularly suited for Bayesian analyses, since it naturally provides a sample of the posterior distribution, i.e., the product of the likelihood times the priors for the parameters. Nonetheless, the results of the scan of the parameter space from MultiNest can also be used in the frequentist framework, provided that the posterior, and in turn the likelihood, has been explored in enough detail. The advantage of the frequentist interpretation is that the derived constraints are not dependent on the prior chosen to explore the various parameters. We will thus adopt the frequentist interpretation in the following. Within this framework marginalisation over parameters is performed with the profile likelihood method [74] and contours at a certain confidence level are drawn following the expectation of a χ 2 distribution. For example, for contour plots in two dimensions we first derive the profile likelihood in the two given parameters profiling over the remaining ones. We then draw contours around the best-fit at 1, 2, 3 and 4σ confidence level according to a two-dimensional χ 2 distribution. A further advantage of using the frequentist formalism is that the output of different MultiNest scans can be easily combined, as we indeed do in figure 5 (and subsequent figures) where different densites of points arise from different separate scans. A disadvantage of this procedure is, of course, that the density of points, which in the Bayesian interpretation has a precise meaning (i.e., it traces the posterior distribution) now loses any meaning (i.e., in figure 5 and subsequent figures only the color of the points is important, not the density). To ensure that the likelihood is well sampled, we run MultiNest with high-accuracy settings, using between 1000 and 3000 live points, depending on the scan, a typical tolerance tol = 0.001, and an enlargement factor between efr = 0.3 − 0.5 in order to ensure that also the tails of the distribution are well explored.
Previous analyses of the GCE within the scalar Higgs portal model have considered a simple dark sector with a scalar WIMP particle that constitutes all of DM, and no uncertainty in the DM density profile. In our analysis, such a simplified scenario corresponds to the special case where R = ρ WIMP /ρ total = 1 and J 40 • = J 40 • , nom . Figure 5 shows the results of the corresponding fit of the model parameters m S and λ HS , with R = 1 and J 40 • = J 40 • , nom fixed. We also show the derived parameter σv . A good fit is only achieved in a narrow Y-shaped region for which σv is of the order of 10 −26 cm 3 /s. Figure 5 also shows that the annihilation cross section σv has to increase with increasing m S (see lower left panel) as to provide an overall γ-ray flux consistent with the GCE, compensating the reduced DM density which decreases as 1/m 2 S . On the other hand, as discussed in the previous section, the dark sector may well be more complex than containing just one WIMP DM species. We thus allow R = ρ WIMP /ρ total ≤ 1 and consider R an input parameter for our fit. Moreover, we include the astrophysical uncertainties in the J-factor as explained above. Our full fit of the GCE within the Higgs scalar portal model thus contains the parameters of the model, m S and λ HS , and two additional parameters related to the astrophysical scenario, R and log(J 40 • /J 40 • , nom ). Figure 6 shows the result of the complete fit of the GCE for the four input parameters m S , λ HS , R and log(J 40 • /J 40 • , nom ). Instead of the narrow Y-shaped stripe with λ HS 0.1 which we observe for the case R = 1 (figure 5), the region of larger λ HS is now also allowed as R < 1 can compensate for the larger annihilation cross section (see upper left panel). In addition, as R = 1 relaxes the tight connection between the normalization (governed by σv R 2 ) and the shape of the photon spectrum (governed by λ HS for a given mass m S ), a new region appears for DM masses above the Higgs threshold, m S ≥ m h . At the hh-threshold the spectral shape for annihilation into hh fits better than for W W and so large λ HS are preferred. However, for large λ HS the overall γ-ray flux is too large to accommodate the GCE signal unless R is small. Hence, for m S > m h the fit prefers small R, see middle left panel in figure 6.

Constraints
The scalar singlet Higgs portal model is constrained by the Higgs invisible width, by direct detection searches, by searches for γ-rays from dwarf spheroidal galaxies, by searches for γ spectral lines, and by the DM relic density. We will discuss these constraints in turn below and quantify their impact by including the corresponding likelihoods into the global fit of the GCE signal. We have checked that constraints for CMB anisotropies [75] on the annihilation cross section are less constraining than limits from dwarf spheroidal galaxies in the relevant region of parameter space (10 GeV m S 1 TeV). Limits on the scalar Higgs portal model from CMB anisotropies have been considered in [29,33].

Higgs invisible width
For light scalar DM particles below the Higgs threshold, m S < m h /2, the decay h → SS results in an invisible Higgs width, Γ inv . This region of parameter space is thus constrained by the LHC limits on the Higgs invisible branching ratio, BR inv 0.23 [76].
In the scalar Higgs portal model, the invisible Higgs width is [77] Assuming that the visible Higgs decay width is given by the Standard Model width, Γ SM , so that BR inv = Γ inv /(Γ SM + Γ inv ), the upper limit on BR inv implies an upper limit on Γ inv and thus on the Higgs-scalar coupling λ HS as a function of the DM mass. For the numerical analysis we have used the value Γ SM = 4.03 MeV [56] for the Standard Model Higgs width and the log-likelihood function for BR inv provided by the ATLAS analysis [76].

Direct detection
A spin-independent DM-nucleon scattering cross section is predicted within the scalar Higgs portal model through the exchange of the SM Higgs boson. The model is therefore severely constrained by direct detection experiments. As the mass of the SM Higgs is large compared to the momentum transfer in the elastic DM-nucleon scattering, the cross section can be described by an effective interaction and is given by [29] σ SI = λ 2 Here, f N = 0.30 [29] denotes the strength of the effective Higgs-nucleon interaction, and µ r = m N m S /(m N + m S ) is the DM-nucleon reduced mass. The current best limits on σ SI come from the LUX experiment [78]. To obtain the likelihood and p-value of the LUX direct detection limits, we have used the tool LUXCalc [79], where the likelihood is constructed from a Poisson distribution. For more details we refer to Ref. [79]. In section 5 we shall also comment on the projected sensitivity of future direct detection experiments like XENON1T [80] or DARWIN [81]. Note that in contrast to indirect detection, a direct detection signal scales linearly with the DM density and hence linearly with R.

Indirect detection: dwarf spheroidal galaxies and spectral γ lines
The GCE can be tested using other independent γ-ray observations and analyses. At present, observations of dwarf satellite galaxies of the Milky Way provide the most stringent limits on σv . There is indeed a mild tension between the DM interpretation of the GCE and dwarf limits [82]. However, taking into account the respective uncertainties of the J-factors of the galactic center and of the dwarf galaxies, the GCE signal can be accommodated, as we shall quantify below.
To implement the dwarf constraints we use the tabulated likelihood as function of flux for each dwarf provided in [82]. We write the likelihood as a product of likelihoods over each single dwarf as described in [82] and [83]. In particular we consider the seven most constraining dwarfs: Willman 1, Ursa Minor, Ursa Major II, Segue 1, Draco, Coma Berenices and Bootes I. The likelihood of each dwarf contains a factor which depends on σv and the provided tabulated flux likelihood, and a further log-normal factor describing the uncertainty in the J-factor, J i , of the dwarf. For the latter, we use the nominal J i and uncertainty provided in [82,83]. The seven J i of the dwarfs are profiled during our global fit, i.e., for each point sampled in parameter space we tabulate on-the-fly the dwarf's likelihood as function of J i and take the corresponding maximum likelihood value.
We also use the results from the latest search for γ-ray lines in the inner galaxy [84], which set constraints on the annihilation cross section into mono-chromatic photons, σv γγ . In [84] different regions of interest (ROI) are analysed, each maximising the sensitivity to a possible DM annihilation signal for different DM density profiles. We use the results for both the R3 and R16 ROI, corresponding to a circular region of 3 • and 16 • of radius, respectively, see [84]. These ROIs maximise the sensitivity for a generalised NFW profile with an inner slope γ = 1.3 and for an Einasto profile which is slightly more cored, respectively. To derive the likelihood from searches for γ-ray lines as function of the integrated flux in the given ROI, we take from [84] for each energy the fluxes corresponding to the 95% upper limit (∆ log L = 2.71/2), to ∆ log L = 1.0/2 and to ∆ log L = 0.0. 2 We then approximate the log-likelihood for each energy as a parabola passing through the above points when the minimum is at a positive flux, or as a line when the minimum is at zero flux.
To translate the above likelihood for the flux into a likelihood for the annihilation cross section σv γγ , the J-factors corresponding to the ROIs J 3 • and J 16 • are needed, see eq. (10). Using the Monte Carlo procedure described in section 3.3 we find, not surprisingly, a strong correlation between J 40 • and J 3 • , J 16 • , which is shown in figure 7. This correlation is well approximated by a forth order polynomial in log-log space, also shown in figure 7. Thus, instead of considering J 3 • and J 16 • as further independent nuisance parameters of the fit, we only use J 40 • as free parameter, with the log-normal distribution described in section 3.3 accounting for the uncertainty related to the DM profile. J 3 • and J 16 • are considered functions of J 40 • with no further intrinsic uncertainties.
We calculate σv γγ using the Higgs effective Lagrangian as described in section 3.2.

Relic density
Assuming a standard cosmological history, we can link the relic WIMP density from the thermal freeze-out to the DM density as measured by Planck, Ωh 2 | Planck = 0.1198 ± 0.0015 [85]. The total DM density predicted by our model is We compute Ωh 2 | WIMP with micrOMEGAs. For details regarding the model implementation we refer to section 3.2. We include the effective Higgs-gluon coupling g eff hgg , eq. (4), which depends on the center-of-mass energy of the annihilation process, s. In contrast to the case of annihilation today, for the computation of the annihilation cross section during freeze-out √ s = 2m S is not always a good approximation: for 2m S just below m h , the dominant contribution to the thermally averaged cross section σv comes from the resonance, √ s = m h . Hence, for the computation of g eff hgg , eq. (4), we choose √ s = m h for m S < m h /2.
We compute the χ 2 for the relic density constraint from where we assume that the dominant uncertainty comes from the theoretical prediction of the relic density, σ rel = 10%. The respective log-likelihood reads − 2 log L Ω = χ 2 Ω + 2 log(σ rel Ωh 2 | DM, total ), (18) again, up to an irrelevant constant.

Results and discussion
In this section we shall present the results of a global fit to the GCE, taking into account the constraints from the invisible Higgs branching ratio, direct detection limits, independent searches for γ-rays from dwarf satellite galaxies, searches for spectral γ lines, and from the dark matter relic density, as discussed in section 4. In section 3.5 we have shown that the GCE signal can be well described by the scalar Higgs portal model. We found that a small DM mass, 30 GeV m S 100 GeV, provides the best fit, and that small values of the Higgs-scalar coupling, λ HS 10 −4 , are viable near the Higgs resonance m S ≈ m h /2. Parameter regions above the Higgs threshold, m S m h are viable, too, but only if we allow for a significant non-WIMP contribution to DM, corresponding to small values of R = ρ WIMP /ρ total . For m S m h and large λ HS of O(1), the hh final state is dominant and provides a good description of the shape of the γ-ray spectrum. Small values of R are required in this parameter region to reconcile the corresponding large annihilation cross section σv with the GCE flux ∝ σv × R 2 .
We now consider the effect of the various constraints discussed in section 4 by including the corresponding likelihoods in the GCE fit. The results of these global fits are presented in figures 8 and 9; figure 8 provides information on the viable parameter space in m S , λ HS , R = ρ WIMP /ρ total and log(J 40 • /J 40 • , nom ), when adding successively the constraints from the invisible Higgs branching ratio, the LUX direct detection limit, the limits from dwarf satellites and the limits from spectral γ lines to the GCE fit. Figure 9 shows detailed information on the results of a global fit to the GCE signal with all constraints added, including in particular the DM relic density. We shall now discuss the impact of the various constraints in turn.
• DM masses below the Higgs resonance, m S < m h /2, lead to invisible Higgs decays, h → SS, and are thus constrained by the LHC limit on the invisible Higgs branching ratio, see section 4.1. The limit on invisible Higgs decays cuts out the region λ HS 0.02 for m S m h /2, and the best fit point for the GCE signal moves to the resonant region, m S ≈ m h /2, see figure 8 upper left panel. The parameter region above the Higgs threshold, m S m h , is still viable and lies within 1σ of the best fit point.
• Severe constraints on the model parameter space are imposed by the current direct detection limits from LUX [78], see section 4.2. These limits exclude the parameter space corresponding to λ HS 0.02 (for R = 1) and λ HS 0.5 (for R ≃ 10 −3 ), and in particular remove the Higgs threshold region. As shown in the upper right panel of figure 8, only the region near the resonance, m S ≈ m h /2, and a small and barely viable parameter region at the W -threshold survives. The region at the W -threshold is, however, already between 3 and 4σ away from the best-fit point.
• Adding the constraints from dwarf galaxy γ-ray searches as described in section 4. 1 are now disfavoured, as a larger J 40 • implies an even larger value for the J-factors of the γ-lines, J 3 • and J 16 • , see figure 7. The fit presented in figure 8 takes into account the log-likelihood from R3 which is optimised for a 'cuspy' dark matter profile as considered here. Indeed we found that R3 provides a stronger constraint than R16 in the considered region of parameter space.
Finally, we discuss the constraint from the DM relic density, see section 4.4. Requiring Ωh 2 | DM, total = Ωh 2 | WIMP /R = Ωh 2 | Planck further constrains the parameter space and leads to interesting parameter correlations. Note, however, that the connection between the DM annihilation cross section σv and the relic density Ωh 2 is based on the assumption of a standard cosmological history. Deviations from the standard cosmological scenario could lead to both a reduction or an enhancement of the predicted relic density. In order to understand this structure we first consider the more simple case R = 1, and then the general case R ≤ 1.
• If we require that the scalar Higgs portal WIMP constitutes all of DM, i.e., R = 1, only a small viable region remains with λ HS ≈ 2 × 10 −4 near the very tip of the resonance at m S ≈ m h /2. At this point of the parameter space, the annihilation cross sections as of today, σv today , and at the time of freeze out, σv freeze-out , are of the same order, and the flux required to describe the GCE and the cosmological dark matter relic density can be both accommodated at the same time. For dark matter masses slightly above the resonance tip, however, the ratio σv today / σv freeze-out increases rapidly, since the annihilation cross section σv exhibits a strong velocity dependence in the vicinity of the s-channel Higgs resonance. Thus, for m S m h /2, the flux required to describe the GCE implies a reduced σv freeze-out and in turn dark matter relic densities which exceed the observed cosmological abundance, so that R = 1 is not possible. Further away from the resonance σv today / σv freeze-out approaches values of order one again, and a second region of parameter space opens up where it is possible to reconcile the GCE flux and the DM relic density. This region, however, is in tension with direct detection constraints.
• The possibility of a non-WIMP dark matter component, corresponding to R < 1, opens up more parameter space. The increase in the ratio σv today / σv freeze-out , when moving from the very tip of the resonance at m S ≈ m h /2 to slightly larger dark matter masses, can be compensated by a decrease in R.
Since Ω DM, total ∝ Ω WIMP /R ∝ 1/(R σv freeze-out ), and σv freeze-out ∝ λ 2 HS , requiring Ω Planck = Ω DM, total implies that R ∝ 1/λ 2 HS for a given mass. This anti-correlation between R and λ 2 HS is clearly visible as the narrow stripe in the R-λ HS panel of figure 9.
As mentioned above, slightly further away from the resonance, σv today / σv freeze-out approaches values of order one again, and thus R ≈ 1 would be viable. However, due to the direct detection constraint which disfavours large R × λ HS , the region with R = 1 is already 3σ away from the best fit point and R 0.5 is preferred in this region. It is interesting to note that the viable regions of parameter space extend down to R ≈ 10 −3 , i.e., the GCE signal can be explained by a thermal relic WIMP that only constitutes one per mille of the dark matter.
To describe the GCE γ-ray flux from the annihilation of the scalar WIMP, the cross section has to be σv × R 2 ≈ 1 × 10 −26 cm 3 /s . The correlation between σv × R 2 and the four fit parameters m S , λ HS , R and log(J 40 • /J 40 • , nom ) is shown in figure 10 for the case of the GCEonly fit (upper panels) and the fit including all the constraints (lower panels). Adding the constraints reduces the viable region of model parameter space to scalar masses m S ≈ m H /2 and the possible values of Higgs-scalar coupling to λ HS 0.1. Furthermore, the allowed range in log(J 40 • /J 40 • , nom ) is partly reduced. As expected there is a strong anti-correlation between σv × R 2 and log(J 40 • /J 40 • , nom ). A smaller σv × R 2 is allowed for a larger log(J 40 • /J 40 • , nom ), and vice-versa.
Besides the overall flux, the model also has to accommodate the spectral shape of the GCE. In figure 11 we show the energy spectrum of the γ-ray flux for the GCE-only fit and the fit including all the constraints, compared to the Fermi-LAT data as analysed in [15]. The correlated systematic errors from astrophysical uncertainties are shown as the red shaded areas, while the error bars denote the uncorrelated statistical errors, see [15] for details. Both fits predict γ-ray spectra which are systematically lower than the mean values of the GCE spectrum. However, as the systematic astrophysical uncertainties (red shaded areas) are strongly correlated, both the GCE only fit and the fit including all constraints provide a good description of the GCE energy spectrum, as discussed in more detail below. Note that the systematic errors depend Figure 10: The correlation between the annihilation cross section as of today, σv , and the input parameters of the fit, mS, λHS, R and log(J 40 • /J 40 • , nom ) for the case of the GCE-only fit (upper row) and taking into account all the constraint, including the relic density constraint (lower row). The white dot denotes the best-fit point. The dark-red, red, orange and yellow points lie within the 1, 2, 3 and 4σ region around the best-fit point, respectively.
on the theory prediction as explained in section 3.5. In Fig. 11 we show the errors for the fit including all constraints although the difference to the GCE-only fit would be barely noticeable in the figure. In table 2 we have collected the best fit values for the scalar Higgs portal model parameters m S , λ HS , and for the astrophysical parameters R = ρ WIMP /ρ total and J 40 • . We show results for the various global fits, including the GCE signal only, and for the GCE signal with the different constraints added successively. As discussed above, the scalar Higgs portal model can describe the GCE signal for dark matter masses near m S = m h /2, and for perturbative values of the Higgs-scalar coupling, λ HS 10 −2 . After taking into account all constraints two viable regions of parameter space emerge: one region where λ HS ≈ 2 × 10 −2 and where an additional dark matter component beyond the scalar WIMP of our model is preferred (R 0.5), and one region where the scalar WIMP constitutes all of dark matter, R 1, and where λ ≈ 2 × 10 −4 λ HS 10 −2 . For both these regions, the best fit point has χ 2 ≈ 27 with respect to the GCE signal. For the 25 data points and the 4 parameters, this corresponds to a p-value of p = 0.18, and thus indicates a reasonably good fit to the GCE.
In table 2 we also list the p-values corresponding to all the constraints we use, i.e., more precisely, the confidence level at which the best fit is compatible with the constraints coming from each single additional observable (BR inv , LUX, dwarfs, lines and DM relic density) we include in the fit. 3 We find that for both viable regions of parameter space the GCE global fit is well compatible with all the constraints. A light tension is present only with respect to the limits from the dwarf galaxies; the corresponding p-value of p = 0.23, is however well within an acceptable range. Note, nonetheless, that we have assumed nominal values from [82,83] for the uncertainties in the J-factors of the dwarfs. On the other hand, it has been pointed out that this uncertainty has been possibly underestimated [86]. This would contribute to relax the dwarf constraints and alleviate the mild tension.
Finally, we have studied the impact of future direct detection experiments like XENON1T [80] or DARWIN [81]. We simply assume a sensitivity to σv which is 10 or 50 times larger than that of the current LUX limits, and include those potential future limits in our global fit. The result is displayed in figure 12 where we show the viable regions in the R-λ HS plane given the current LUX bounds (left panel), and a 10 times (middle panel) and 50 times (right panel) larger potential future sensitivity. The future direct detection experiments will probe part of the currently allowed region of parameter space, and stronger limits would exclude regions with both λ HS and R large. On the other hand, even limits corresponding to a 50 times larger sensitivity than that of current direct detection experiments would leave viable regions of model parameter space, corresponding to either sizeable Higgs-scalar couplings λ HS ≈ 2 × 10 −2 and small R 0.1, or to smaller λ HS but values of R close to one. Note that future limits from CTA [87] will be relevant for DM masses above about 100 GeV, see e.g. [88,29,33], and would thus not constrain the GCE interpretation within the scalar Higgs portal model further. Also LHC searches for scalar Higgs portal models with dark matter masses m S m h /2 are not sensitive to the relevant model parameter space, see [36][37][38][39]. The best test of these two regions will probably come from further γ-ray searches. Limits from both dwarf galaxies and lines are expected to improve with time while larger statistics is being collected. In addition, limits from dwarfs will further improve since more dwarf galaxies, potentially with large J-factors, should be discovered in the next years by current surveys like DES [89], and future ones like LSST [90]. This will clarify if the present tension with the GCE will become more severe at the point of excluding the remaining parameter space, or, if, eventually, a signal in dwarfs and lines will emerge, confirming the GCE interpretation.   Table 2: Best fit points and corresponding 1σ error for fits to the GCE only, and including successively constraints from the invisible Higgs branching ratio, direct detection limits, independent searches for γ-rays from dwarf satellite galaxies, searches for spectral γ lines, and from the dark matter relic density. We also display the best fit in a second, viable region of parameter space (last column). Also shown are the χ 2 GCE and the p-values of the respective best-fit points taking into account the log-likelihood contributions of the observables given in the first line. p(χ 2 GCE ) represents the goodness of the GCE fit for 25 data points and 4 fitted parameters. The remaining p-values represent the confidence level at which the best fit is compatible with the constraints coming from each extra-observable we include in the fit (see text for more details). λHS R λHS R λHS R Figure 12: Results of a fit to the GCE with free parameters mS, λHS, R and log(J 40 • /J 40 • , nom ). We take into account all constraints (including the relic density constraint) and focus on the correlation between R and λHS. In the middle and right panel we include potential future direct detection limits with a sensitivity of 10 and 50 times the current LUX sensitivity, respectively. The white dot denotes the best-fit point. The dark-red, red, orange and yellow points lie within the 1, 2, 3 and 4σ region around the best-fit point, respectively.

Conclusion
We have presented a global fit of the γ-ray galactic center excess (GCE) within a minimal Higgs portal model with a scalar WIMP dark matter particle. We find that the shape and strength of the GCE is well described by dark matter annihilation in various regions of parameter space, including in particular the resonance and threshold regions where the dark matter mass, m S , is close to (half the) Higgs and W boson masses, respectively.
The parameter space of the scalar Higgs portal model is constrained by the search for invisible Higgs decays, direct dark matter searches, searches for dark matter annihilation from dwarf spheroidal galaxies, searches for mono-energetic spectral γ-lines from the Milky Way halo, and by the cosmological dark matter relic abundance. We have included these constraints into the global fit of the GCE signal and studied the implications for the model parameter space, taking properly into account the theoretical uncertainty from the dark matter distribution. Furthermore, we consider the possibility that the dark sector is more complex than that of the minimal Higgs portal model and allow for scalar WIMP dark matter densities smaller than the total, gravitationally interacting dark matter density. With this freedom, we can easily accommodate the GCE signal in regions of model parameter space which would otherwise be excluded because the model prediction would exceed the cosmologically observed dark mater relic density.
Taking into account all constraints, the scalar Higgs portal model can describe the GCE signal if the dark matter mass is near m S = m h /2, so that dark matter annihilation proceeds through resonant Higgs exchange. Two regions of parameter space are viable: one region where the Higgs-dark-matter coupling, λ HS , is of order O(10 −2 ) and where an additional dark matter component beyond a scalar WIMP is preferred, and one region where λ HS may be significantly smaller, but where the scalar WIMP constitutes a significant fraction or even all of dark matter. These regions emerge from an interplay between the different scaling of the GCE signal and the relic density with the fraction of WIMP dark matter, R = ρ WIMP /ρ total , and the strong velocity dependence of the annihilation cross section near the resonance. Note that this effect can potentially be of relevance for every model with resonant annihilation. The favoured regions of scalar dark matter masses and couplings are hard to probe in future direct detection and collider experiments. Future searches for γ-ray emission from dwarf spheroidal galaxies, however, will be able to confirm or exclude the dark matter interpretation of the galactic center excess.