The inert doublet model in the light of Fermi-LAT gamma-ray data: a global fit analysis

TL;DR

This study conducts a comprehensive global fit of the inert doublet model, integrating collider, dark matter search data, and gamma-ray observations to identify parameter regions compatible with current constraints and the Galactic center excess.

Contribution

It presents the first combined analysis considering recent gamma-ray data and relaxes the assumption that the inert doublet accounts for all dark matter, identifying viable parameter spaces.

Findings

Two parameter regions fit all constraints and explain the Galactic center excess.

One region near the Higgs resonance, another around 72 GeV with specific annihilation channels.

Future experiments could further test these viable scenarios.

Abstract

We perform a global fit within the inert doublet model taking into account experimental observables from colliders, direct and indirect dark matter searches and theoretical constraints. In particular, we consider recent results from searches for dark matter annihilation-induced gamma-rays in dwarf spheroidal galaxies and relax the assumption that the inert doublet model should account for the entire dark matter in the Universe. We moreover study in how far the model is compatible with a possible dark matter explanation of the so-called Galactic center excess. We find two distinct parameter space regions that are consistent with existing constraints and can simultaneously explain the excess: One with dark matter masses near the Higgs resonance and one around 72 GeV where dark matter annihilates predominantly into pairs of virtual electroweak gauge bosons via the four-vertex arising from the inert doublet's kinetic term. We briefly discuss future prospects to probe these scenarios.

Tables

Table 1. Table 1: J 𝐽 J -factors and their uncertainties for the dwarf spheroidal galaxies considered in this study. For the first seven dwarfs we use the measured values from [ 87 ] . For Tucana III and Tucana IV J 𝐽 J -factors are estimated, see text for details.

dwarf	${\log }_{10}\left(J_{\text{meas}}\right)\left[{\log }_{10}\left({\text{GeV}}^2{\text{cm}}^{-5}\right)\right]$
Coma Berenices	$19.0\pm 0.4$
Draco	$18.8\pm 0.1$
Sculptor	$18.5\pm 0.1$
Segue 1	$19.4\pm 0.3$
Ursa Major II	$19.4\pm 0.4$
Ursa Minor	$18.9\pm 0.2$
Reticulum II	$18.9\pm 0.6$
Tucana III	$19.3\pm 0.6$
Tucana IV	$18.7\pm 0.6$

Table 2. Table 2: Parameters of the fit and the corresponding allowed ranges. ​ ∗ Additionally we require m A 0 , m H ± > m H 0 subscript 𝑚 superscript 𝐴 0 subscript 𝑚 superscript 𝐻 plus-or-minus subscript 𝑚 superscript 𝐻 0 m_{A^{0}},m_{H^{\pm}}>m_{H^{0}} .

parameter	range
$m_{H^0}$	$[45;1000]\text{GeV}$
$m_{A^0}$	$[{100}^{\ast };1100]\text{GeV}$
$m_{H^{\pm }}$	$[{100}^{\ast };1100]\text{GeV}$
${\lambda }_L$	$[-4\pi ;\mathrm{\hspace{0.17em}4}\pi ]$
${\lambda }_2$	$[{10}^{-6};\mathrm{\hspace{0.17em}4}\pi ]$
$\log \left(J/J_{\text{nom}}\right)$	$[-4{\sigma }_{\log J};4{\sigma }_{\log J}]$
$R$	$[{10}^{-3};1]$

Table 3. Table 3: Parameters describing the GCE best-fit regions appearing in figure 2 . The regions 1a and 1b correspond to the Higgs funnel with m H 0 ≈ m h 0 / 2 subscript 𝑚 superscript 𝐻 0 subscript 𝑚 superscript ℎ 0 2 m_{H^{0}}\approx m_{h^{0}}/2 . Region 2 refers to the case where DM annihilates predominantly into virtual gauge boson pairs. The last line indicates the dominant annihilation channel and its relative contribution to annihilation today.

	region 1a	region 1b	region 2
$m_{H^0}\left[\text{GeV}\right]$	${62.522}_{-0.004}^{+0.017}$	${62.61}_{-0.06}^{+0.23}$	${71.95}_{-0.48}^{+1.18}$
${\lambda }_L$	${0.00009}_{-0.00001}^{+0.00039}$	${0.011}_{-0.002}^{+0.003}$
$R$	${0.99}_{-0.97}^{+0.01}$	${0.02}_{-0.02}^{+0.07}$	${1.00}_{-0.19}^{+0.00}$
$\log \left(J/J_{nom}\right)$	${0.25}_{-0.38}^{+0.31}$	${0.24}_{-0.37}^{+0.33}$	${0.19}_{-0.16}^{+0.25}$
${\chi }_{\text{GCE}}^2$	$26.9$	$27.3$	$33.6$
$p\left({\chi }_{\text{GCE}}^2\right)$	$0.22$	$0.20$	$0.054$
$p\left(\text{LUX}\right)$	$1$	$0.79$	$0.67$
$p\left(\text{dSph}\right)$	$0.33$	$0.33$	$0.34$
$p\left(\text{rel.D.}\right)$	$0.87$	$0.92$	$0.99$
$p\left(\text{ST}\right)$	$0.68$	$0.67$	$0.68$
dom. channel	$b\overline{b}\left(81\%\right)$	$b\overline{b}\left(67\%\right)$	$W{W}^{\ast }\left(89\%\right)$

Equations

$$V~{}=~{}\mu_{1}^{2}|H|^{2}+\mu_{2}^{2}|\Phi|^{2}+\lambda_{1}|H|^{4}+\lambda_{2}|\Phi|^{4}+\lambda_{3}|H|^{2}|\Phi|^{2}+\lambda_{4}|H^{\dagger}\Phi|^{2}+\frac{\lambda_{5}}{2}\Bigl{[}(H^{\dagger}\Phi)^{2}+\mathrm{h.c.}\Bigr{]}.$$

$$H~{}=~{}\left(\begin{array}[]{c}G^{+}\\ \frac{1}{\sqrt{2}}\left(v+h^{0}+\mathrm{i}G^{0}\right)\end{array}\right),\qquad\Phi~{}=~{}\left(\begin{array}[]{c}H^{+}\\ \frac{1}{\sqrt{2}}\left(H^{0}+\mathrm{i}A^{0}\right)\end{array}\right),$$

$$m_{h^0}^2$$

$$m_{H^0}^2$$

$$m_{A^0}^2$$

$$m_{H^{\pm }}^2$$

$${\lambda }_{L,S}$$

$$\left\{{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\mu }_2\right\},$$

$$\left\{m_{h^0},m_{H^0},m_{A^0},m_{H^{\pm }},{\lambda }_L,{\lambda }_2\right\},$$

$$R\equiv \frac{{\rho }_{\text{IDM}}}{{\rho }_{\text{total}}}$$

$$\Omega h^2{|}_{\text{DM,total}}=\frac{\Omega h^2{|}_{\text{IDM}}}{R}.$$

$${\chi }_{\Omega }^2=\frac{{\left(\Omega h^2{|}_{\text{DM,total}}-\Omega h^2{|}_{\text{Planck}}\right)}^2}{{\left({\sigma }_{\text{rel}}\times \Omega h^2{|}_{\text{DM,total}}\right)}^2},$$

$$-2\log {\mathcal{L}}_{\Omega }={\chi }_{\Omega }^2+2\log ({\sigma }_{\text{rel}}\Omega h^2{|}_{\text{DM,total}}),$$

$${\sigma }_{\text{SI}}=\frac{{\lambda }_L^2{f}_{\text{N}}^2}{4\pi }\frac{{\mu }_{\text{r}}^2{m}_{\mathrm{N}}^2}{m_{h^0}^4{m}_{H^0}^2},$$

$$E^2\frac{\mathrm{d}\varphi }{\mathrm{d}E}=\frac{1}{4\pi }\frac{⟨\sigma v⟩R^2}{2m_{\text{DM}}^2}{\int }_{E_{\min }}^{E_{\max }}\mathrm{d}{E}_{\gamma }{E}_{\gamma }\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}{E}_{\gamma }}\times {\int }_{\text{ROI}}\mathrm{d}\Omega {\int }_{\text{l.o.s}}\mathrm{d}s{\rho }_{\text{DM}}^2,$$

$$J_{40^{\circ }}=\underset{\text{ROI}}{\int }\mathrm{d}\Omega \underset{\text{l.o.s}}{\int }\mathrm{d}s{\rho }_{\text{DM}}^2,$$

$$\log \left(\frac{J_{40^{\circ }}}{{\text{GeV}}^2{\text{cm}}^{-5}}\right)=23.25\pm 0.43.$$

$${\chi }_{\text{GCE}}^2=\sum _{i,j}(d_i-t_i){\left({\Sigma }_{ij}+{\delta }_{ij}({\sigma }_{\text{rel}}{t}_i{)}^2\right)}^{-1}(d_j-t_j)+\frac{(\log J_{40^{\circ }}-\log J_{40^{\circ },\text{nom}}{)}^2}{({\sigma }_{\log J}{)}^2},$$

$$-2\log {\mathcal{L}}_{\text{GCE}}={\chi }_{\text{GCE}}^2+\log |{\Sigma }_{ij}+{\delta }_{ij}({\sigma }_{\text{rel}}{t}_i{)}^2|,$$

$$\Gamma \left(h^0\to H^0{H}^0\right)=\frac{{\lambda }_L^2{v}^2}{8\pi m_{h^0}}\sqrt{1-\frac{4m_{H^0}^2}{m_h^2}}$$

$${\chi }_{ST}^2={\mathbf{v}}^T{\mathcal{C}}^{-1}\mathbf{v}$$

$$\hat{S}=0.06,\hat{T}=0.097$$

$$\mathcal{C}=\left(\begin{array}{cc}0.0085& 0.0063\\ 0.0063& 0.0057\end{array}\right).$$

$$m_{H^{\pm }}\gtrsim m_W.$$

$$m_{A^0}\gtrsim 100\mathrm{GeV}.$$

$$\frac{{⟨\sigma v⟩}_{\mathrm{today}}}{{⟨\sigma v⟩}_{\mathrm{freeze}-\mathrm{out}}}\times R\sim \mathcal{O}(1).$$

Full text

1. 1 Introduction

2. 2 The inert doublet model

3. 3 Constraints and global fit settings

4. 3.1 Dark matter observables

5. 3.1.1 Relic density

6. 3.1.2 Direct detection

7. 3.1.3 Indirect detection constraints: dwarf spheroidal galaxies

8. 3.1.4 The Galactic center excess

9. 3.2 Non-dark matter observables

10. 3.2.1 Unitarity, perturbativity and vacuum stability

11. 3.2.2 Higgs invisible width

12. 3.2.3 Electroweak precision observables

13. 3.2.4 LEP-II bounds on the masses of the heavy ${\cal{Z}}_{2}$-odd states

14. 3.3 Scan settings

15. 4 Results and discussion

16. 4.1 Global dark matter fit

17. 4.2 Fitting a possible signal from dark matter annihilation

18. 4.3 Future prospects for direct, indirect detection and the LHC

19. 5 Conclusion

20. Acknowledgements

1 Introduction

The nature of dark matter (DM), the existence of which is corroborated by observations over a wide range of physical scales in the Universe, is one of the most important open questions in contemporary fundamental physics. An explanation in terms of a weakly interacting massive particle (WIMP) is an attractive possibility which has motivated an enormous experimental effort. Indirect detection experiments have reached sensitivity to the thermal self-annihilation cross section for DM masses around the electroweak scale and direct detection experiments have substantially improved the limits on WIMP-nucleon scattering over the past few years. Interpreting these results in terms of well-motivated theoretical models is hence an important task in order to pinpoint the nature of DM.

The inert doublet model (IDM)444For recent accounts see, e.g. [1, 2]. is among the simplest new physics models, supplementing the standard model with an additional complex scalar field that transforms as a doublet under $SU(2)_{L}$ and is odd under a discrete ${\cal{Z}}_{2}$ symmetry, all standard model fields being taken to be even. Despite its simplicity the IDM has a rich and versatile phenomenology: it can affect electroweak symmetry breaking [3, 4, 5, 6, 7, 8], give rise to interesting, observable effects at colliders [2, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], modifiy electroweak baryogenesis [19, 20], play a role in the generation of neutrino masses [21] and, as being the focus of this work, it contains a WIMP that can account for the observed DM in the Universe with observable signatures in direct and indirect detection experiments [4, 17, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. The versatility of the IDM as a DM model introduces a fair amount of freedom to accommodate measurements and constraints from various observables, making it a non-trivial task to unfold the data and extract information about the physical parameters of the model.

In this regard, global fit techniques are of central importance. They enable the systematic study of the impact of a large number of experimental measurements while fully accounting for systematic uncertainties that affect astrophysical observables such as the DM density profile in the inner galaxy. In this paper we perform a detailed numerical fit within the IDM using MultiNest [34, 35], which allows us to comprehensively explore the model’s parameter space. Furthermore, instead of demanding that the IDM dark matter candidate should account for the entire DM abundance in the Universe, we follow a more general approach allowing for an unspecified additional DM component to contribute subdominantly or even dominantly to the total DM density. Introducing the fractional density of IDM dark matter as a free parameter in the fit enables us to extract information about the amount of DM that can be accommodated within the model.

We consider two setups. On the one hand, we fit a set of well-established observables: The DM relic density measured by Planck [36], direct detection constraints set by LUX [37], indirect detection constraints from the observation of dwarf spheroidal galaxies set by Fermi-LAT [38, 39] as well as the Higgs mass measured at the LHC [40], constraints from invisible Higgs decays [41], constraints from electroweak precision tests [42] and theoretical bounds from unitarity, perturbativity and vacuum stability.

On the other hand, over the past few years the Fermi-LAT gamma-ray data revealed an unexpected Galactic bulge emission component – the “Galactic center excess” (GCE) [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. Although various astrophysical explanations have been proposed [58, 59, 60, 61] (for statistical approaches to test the origin of the signal see further discussions in [62, 63, 64, 65, 66]), it is intriguing that the strength as well as spectral and morphological properties of the excess are compatible with a signal from DM annihilation with thermal cross section and a DM mass $m_{\text{DM}}\lesssim 100\,\text{GeV}$. Given the complexity of the Galactic center as an astrophysical environment, dwarf spheroidal galaxies (dSphs) are a much cleaner target for DM searches in gamma-rays as their dynamical and chemical properties suggest larger mass-to-light ratios. Hence they provide an important test of the DM interpretation of the GCE. Searches for a gamma-ray excess associated with dSphs have been performed and the sensitivity is competitive with that of other targets such as the Galactic center. Interestingly, slight excesses (each with a $\sim 2\sigma$ local significance) have been found in four of the recently discovered dSph targets [67, 68, 39] which are roughly compatible with a DM explanation of the GCE. Given these hints, as a second step in this paper we include the GCE in addition to the observables mentioned above in our global fit.555For other attempts to perform global fits to the GCE within UV-complete models see, e.g. [69, 70]. Note also that similarly, an excess which appears to be compatible with a signal from DM annihilations as well as with the GCE itself has recently been reported [71] in the AMS-02 antiproton data [72] (see [73] for a similar result using the boron over carbon ratio [74]). An interpretation in terms of individual DM annihilation channels but also in terms of the singlet scalar DM model was presented in [75]. Although an analysis of the antiproton flux measurements falls beyond the scope of the present paper, these hints provide additional motivation to study whether the IDM can accommodate the GCE.

The paper is organized as follows: in section 2 we briefly introduce the IDM. The various observables included in our fitting procedure, along with the method followed in order to sample the IDM parameter space, are detailed in section 3. Section 4 contains our main results and a discussion on the future sensitivities of upcoming experiments to the best-fit parameter regions. We conclude in section 5.

2 The inert doublet model

The IDM is a special case of a two-Higgs doublet model in which an exact ${\cal{Z}}_{2}$ discrete symmetry is imposed on the Lagrangian, under which all standard model fields (including the usual Higgs doublet $H$) are taken to be even, whereas the second scalar doublet $\Phi$ is taken to be odd. With respect to the standard model Lagrangian, the only modifications consist of the introduction of gauge kinetic terms for $\Phi$ and an additional piece in the scalar potential, which in the IDM reads

[TABLE]

Upon electroweak symmetry breaking, the two scalar doublets can be expanded in component fields as

[TABLE]

where $v\simeq 246$ GeV is the usual Higgs field vacuum expectation value and $G$ are the Goldstone bosons. The model contains five physical scalar states with masses given by

[TABLE]

where, following common conventions, we have defined

[TABLE]

These parameters correspond to the coupling of a pair of $H^{0}$, $A^{0}$ states, respectively, to the Higgs boson. All in all, the IDM is characterized by six free parameters:

[TABLE]

which, using equations (3), (4) and the potential minimization condition $\left.\left(\partial V/\partial h^{0}\right)\right|_{h^{0}=0}=0$ can be traded for the physically more intuitive set of parameters

[TABLE]

with $m_{h^{0}}\simeq 125.09$ GeV [76] being the measured Higgs boson mass.

The discrete symmetry imposed on the Lagrangian renders the lightest component of the inert doublet stable. If, moreover, this lightest ${\cal{Z}}_{2}$-odd particle is neutral ($m_{H^{0}/A^{0}}<m_{H^{\pm}}$), it can play the role of a DM candidate. In fact, the IDM is perhaps the simplest model in which the observed DM abundance in the Universe can be obtained through all ways that typically characterize WIMP models: adjusting couplings, approaching or taking distance from resonances/thresholds and coannihilation. Note also that the DM phenomenology of $H^{0}$ and $A^{0}$ is identical. In this respect, for simplicity in what follows we will consistently adopt the choice $m_{H^{0}}<m_{A^{0}}$. A more detailed description of the IDM phenomenology will be presented in the following sections.

3 Constraints and global fit settings

Various aspects of the IDM phenomenology have been studied in the literature. The model was first proposed as a DM model in [4], while its predicted relic abundance was analysed in more detail in [22, 23, 24, 25]. Direct detection constraints were first considered in [4] (as well as in most subsequent studies), whereas indirect detection has been studied for continuum gamma-rays [22, 32], spectral features [26, 30, 31], antimatter [29] and neutrinos [27, 28]. Other than the DM abundance, direct detection is known to impose extremely strong constraints on the IDM parameter space whereas currently the strongest indirect detection bounds stem from gamma-ray searches in dwarf spheroidal galaxies. Besides, the region of parameter space where important gamma-ray lines could be expected is severely bound by other observations.

In [17] it was shown that the invisible Higgs width imposes strong bounds on the $\lambda_{L}$ coupling if $m_{H^{0}}<m_{h^{0}}/2$, whereas the masses of the next-to-lightest states can be constrained from LEP-II searches for neutralinos [9] and charginos [10]. The new states can also induce contributions to the $S$ and $T$ electroweak parameters, as first pointed out in [4]. Constraints from searches for dileptons along with missing transverse energy at the LHC were first proposed in [12] and have been analysed in [18]. Although during Run-II they will offer the opportunity to test an interesting part of the IDM parameter space, at present their impact is subleading with respect to other searches. Lastly, a minimal set of requirements must be imposed on the parameter choices in order to ensure that the electroweak vacuum is stable (for detailed analyses cf. [6, 7, 17]) and that perturbative calculations make sense, also in the sense of perturbative unitarity of the scattering matrix.

We now proceed to discuss these constraints and the way they are incorporated in our global likelihood fit in more detail.

3.1 Dark matter observables

3.1.1 Relic density

Assuming a standard cosmological history allows us to link the relic $H^{0}$ density from thermal freeze-out, $\Omega h^{2}|_{\text{IDM}}$, to the DM density measured by Planck, $\Omega h^{2}|_{\text{Planck}}=0.1198\pm 0.0015$ [36]. In this study we allow for the possibility that the dark sector might be comprised of more than one DM component by introducing the fraction of the DM density predicted from the IDM over the total DM density in the Universe

[TABLE]

as a free (astrophysical) parameter. We assume that the clustering properties and, hence, the density profiles of the IDM and non-IDM DM components behave sufficiently similarly so that they constitute the same fraction $R$ of DM on the different scales which are relevant for the various DM observables considered here. Then, the total DM density is given by

[TABLE]

The $H^{0}$ relic density $\Omega h^{2}|_{\text{IDM}}$ has been calculated using micrOMEGAs [77], thanks to an implementation of the model in the FeynRules package [78]. Our computations take into account 3-body final state contributions to the total DM annihilation cross section, which can be extremely important in some regions of parameter space.

We compute the $\chi^{2}$ contribution for the relic density via

[TABLE]

where we assume that the dominant uncertainty originates from the theoretical prediction of the relic density, in particular from the uncertainty on the annihilation cross section. We estimate $\sigma_{\text{rel}}=10\%$ (cf. e.g. [79] for a discussion of the one-loop corrections to the Higgs portal type annihilation and [80] for a relevant discussion within the IDM). The corresponding log-likelihood finally reads

[TABLE]

up to an irrelevant constant.

3.1.2 Direct detection

In the IDM at tree-level spin-independent WIMP-nucleon scattering arises from $Z$- and $h$-exchange [4]. The former is, however, only significant for extremely small mass splittings between $H^{0}$ and $A^{0}$, which are not in the focus of this study. The cross section for WIMP-nucleon scattering via $h$-exchange reads

[TABLE]

where $f_{\text{N}}\sim 0.30$ [81] denotes the strength of the effective Higgs-nucleon interaction, and $\mu_{\rm r}=m_{\rm N}m_{H^{0}}/(m_{\rm N}+m_{H^{0}})$ is the DM-nucleon reduced mass. In this study we compute the total spin-independent scattering cross section at tree-level (including an effective vertex for the Higgs interaction with gluons) with micrOMEGAs [77].

We take into account the most recent constraints from LUX [37]. In order to estimate the respective log-likelihood we utilize the program LUXCalc [82]. However, the current version of LUXCalc is based on the results from LUX 2013 [83]. In order to account for the LUX 2016 sensitivity we proceed as follows. We first determine the (mass-dependent) gain factor in the sensitivity between the LUX 2013 and LUX 2016. Assuming that the improvement in the sensitivity is well described by a simple gain in the exposure we then rescale the signal by this factor and compute the log-likelihood with LUXCalc [82]. Computing the $p$-value in this way allows us to reproduce the limits from LUX 2016 with a relative uncertainty at a per-cent level in the mass region of interest.

It should be noted that as (11) is proportional to $\lambda_{L}^{2}$, for very small $\lambda_{L}$ electroweak corrections induced at 1-loop can be important and can eventually dominate the scattering cross section [84, 85]. However, the magnitude of the electroweak corrections is independent of $\lambda_{L}$ and below $\sim 2\times 10^{-47}\text{cm}^{2}$ for the mass regions considered in this study [84]. Hence, they are well below the current sensitivity of LUX [37] and can be neglected for the computation of the respective likelihood. We will discuss their importance for direct detection projections in section 4.3.

3.1.3 Indirect detection constraints: dwarf spheroidal galaxies

In the low mass region of the IDM where annihilation occurs predominantly into $b\bar{b}$ and/or $W^{+}W^{-}$ pairs, the most stringent limits on the velocity-averaged annihilation cross section, $\langle\sigma v\rangle$, arise from gamma-ray observations of dSphs. We use the results of the recent analysis of the Fermi-LAT data [39].

The predicted $E^{2}\times$ flux in an energy bin between $E_{\min}$ and $E_{\max}$ is

[TABLE]

where $\mathrm{d}N_{\gamma}/\mathrm{d}E_{\gamma}$ is the differential photon spectrum per annihilation and $m_{\text{DM}}$ is the mass of the DM particle. The integral of the DM density, $\rho_{\text{DM}}$, over the region of interest (ROI) and the line of sight (l.o.s) is the $J$-factor, $J_{i}$, of the considered dwarf.

In order to take into account constraints from dSphs in our fit we use tabulated likelihoods for individual dwarfs as a function of the energy flux in the considered 24 energy bins provided in [39]. The total log-likelihood is obtained by summing over the log-likelihood contributions of the individual dwarfs [39, 86]. We take into account the seven dwarfs with the largest confirmed $J$-factors [87] per default, see first seven dwarfs listed in table 1.

We compute the prediction for the binned energy flux by using the tabulated gamma-ray spectra for individual annihilation channels obtained in [88], which we combine according to their relative strength as calculated with micrOMEGAs [77]. These channels include the 3-body final states $WW^{*}$, $ZZ^{*}$, where the virtual vector boson creates a pair of fermions. The uncertainty of the $J$-factors is taken into account by profiling over the $J_{i}$ (for each dwarf) according to its uncertainty provided in table 1, i.e., for each sampled point we tabulate on-the-fly the likelihood of the considered dwarfs as a function of $J_{i}$ and take the corresponding value that provides the maximum likelihood.

In four of the recently discovered dSph targets, slight excesses (each $\sim 2\sigma$ local) have been found: specifically in Reticulum II, Tucana III, Tucana IV and Indus II [67, 68, 39] (but see also [89]). For the latter three targets the dynamical masses have not yet been spectroscopically measured and hence these targets are at present not confirmed as DM-dominated dSphs. However, in order to illustrate the impact of the respective observation on the fit we consider the case of additionally including the log-likelihood of Tucana III and Tucana IV in the fit. For these targets we use the distance-based predictions for the $J$-factors with an estimated error of 0.6 dex [39].

It should be noted that the uncertainties in the $J$-factors used above might be underestimated, comparing e.g. [87, 90] and, in particular, following the discussion in [91]. However, these uncertainties have only a minor impact on our results. Omitting the likelihood contribution from the faint dwarfs Coma Berenices, Ursa Major II and Segue 1 (that exhibit the largest uncertainties [91] among the considered ones) does not qualitatively change our results.

Finally, we note that additional constraints could stem from the Fermi-LAT searches for gamma-ray lines at the Galactic center [92]. Within the IDM line signatures have been studied in [26, 31]. Although the loop-suppression of the production cross section for two monochromatic photons is typically compensated by the higher sensitivity in searches for spectral lines, we do not expect these searches to provide constraints significantly stronger than the ones for a continuous photon spectrum in dwarf spheroidal galaxies considered above.666For instance, for the benchmark point IV in [26] ($m_{h}=120\,\text{GeV}$, $m_{H}^{0}=70\,\text{GeV}$ and $\lambda_{L}\simeq-0.07$) the predicted line signal, $\langle\sigma v\rangle_{\gamma\gamma}=7.6\times 10^{-30}\,\text{cm}^{3}\text{s}^{-1}$, falls below the upper limit, $\langle\sigma v\rangle_{\gamma\gamma}^{\text{UL}}=5.2\times 10^{29}\,\text{cm}^{3}\text{s}^{-1}$ [92], by almost an order of magnitude while the cross section into $b\bar{b}$ alone, $\langle\sigma v\rangle_{b\bar{b}}=1.6\times 10^{26}\,\text{cm}^{3}\text{s}^{-1}$, is already relatively close to the respective limit from dwarfs, $\langle\sigma v\rangle_{b\bar{b}}^{\text{UL}}=2.6\times 10^{26}\,\text{cm}^{3}\text{s}^{-1}$ [39]. Furthermore, the above quoted limit from line searches is derived for the most aggressive choice regarding the dark matter density profile considered in [92], i.e. a generalized Navarro-Frenk-White profile with an inner slope of $\gamma=1.3$, which is, however, compatible with the GCE. For choices of the dark matter density profile that are more cored the limit becomes even weaker. Note also that in [26] the 3-body final state contribution to the continuous gamma spectrum have not been taken into account. However, a full assessment of the importance of line searches in the considered parameter space of the IDM falls beyond the scope of this work.

3.1.4 The Galactic center excess

The Fermi-LAT observation

Over the last few years several groups have reported a Galactic bulge emission component in the Fermi-LAT gamma-ray data [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57] – the GCE. Extending beyond 10∘ away from the Galactic plane, the GCE appears compatible with a spherical morphology and a steep cuspy radial profile [50, 52]. In the $E^{2}\times$ flux representation the inferred energy spectrum peaks at a few GeV. Intriguingly, the excess is compatible with a signal from DM annihilation. In particular it favors an annihilation cross section close to the thermal one, $\langle\sigma v\rangle\sim 10^{-26}\,\text{cm}^{3}\text{s}^{-1}$, which could point towards explanations in terms of thermal WIMPs (but see also [93]). Besides interpretations in terms of DM, various astrophysical explanations of the excess have been proposed [58, 59, 60, 61]. Using statistical methods to indicate whether the photon-count distributions of the excess are compatible with a smooth component, exhibiting Poissonian clustering properties, evidence for an extended gamma-ray point source population have been found [62, 63, 64, 65] disfavoring the DM interpretation of the GCE, see also [94]. However, it remains difficult to control the systematic uncertainties in point source analyses that could arise due to mismodeling of the data [66]. Hence, it is premature to draw a definite conclusion about the origin of the GCE.

In this study we consider the possibility (results presented in section 4.2) that the excess could be entirely due to WIMP annihilation. We use the results of the analysis performed in [52], which provided the inferred energy spectrum along with an error covariance matrix that includes an estimate of systematic uncertainties related to the Galactic foreground emission.

Dark matter density profile and uncertainties

The DM density in the inner part of the Milky Way is subject to large uncertainties affecting the observed gamma-ray flux which is reflected by an uncertainty in the involved $J$-factor

[TABLE]

where the ROI is a $40^{\circ}\times 40^{\circ}$ region centered on the Galactic center with a stripe of $\pm 2^{\circ}$ masked along the Galactic plane [52]. The DM spatial distribution in the Milky Way has been evaluated using dynamical data, for example, in [95, 96] and the ensuing uncertainties on the $J$-factor stem dominantly from the poor knowledge of the inner slope of the DM halo profile (for a recent study of the impact of these uncertainties on the DM-induced gamma-ray flux cf. e.g. [33]). However, the GCE cannot be reproduced unless specific assumptions are made concerning the DM spatial distribution. Concretely, using a generalized Navarro-Frenk-White (NFW) profile [97] with parameters compatible with the measured shape of the GCE, i.e. an inner slope of $\gamma=1.2\pm 0.08$ [52],777For a further discussion of the morphology of the GCE see e.g. [54] and references therein. as well as recent measurements of the scale radius and the scale density [98], based on a Monte Carlo procedure it was shown in [88] that the $J$-factor is

[TABLE]

Note that since the authors of [52] normalize the GCE flux dividing by the angular size of the analyzed region, we divide $J_{40^{\circ}}$ by the corresponding solid angle, $\Delta\Omega=0.43$ sr.

Likelihood for the GCE signal

We compute the $\chi^{2}$ contribution for the GCE – including the contribution from $J_{40^{\circ}}$ – by

[TABLE]

where $d_{i}$ is the GCE measured flux in energy bin $i$ from [52] and $t_{i}$ is the respective model prediction, which depends on the model parameters, $R$ and $J_{40^{\circ}}$. $\Sigma_{ij}$ is the covariance matrix given in [52] and $\log J_{40^{\circ}\!,\,\text{nom}}$ and $\sigma_{\log\!J}$ are the nominal values of the (logarithmic) $J$-factor and its uncertainty, respectively, as given in (14). We compute the predicted flux $t_{i}$ by combining the tabulated gamma-ray spectra for individual annihilation channels(including the 3-body final states $WW^{*}$, $ZZ^{*}$) obtained in [88] weighted by their relative contributions as computed by micrOMEGAs [77]. In addition to the covariance matrix which includes statistical and systematic errors of the observed signal we include $\delta_{ij}(\sigma_{\text{rel}}\ t_{i})^{2}$ representing a diagonal error equal to a fraction $\sigma_{\text{rel}}$ of the model prediction. We choose $\sigma_{\text{rel}}$ = 10% as discussed in [99, 100, 88]. Up to an irrelevant constant factor, the resulting log-likelihood is

[TABLE]

where $|\Sigma_{ij}+\delta_{ij}(\sigma_{\text{rel}}\ t_{i})^{2}|$ is the determinant of the covariance matrix.

3.2 Non-dark matter observables

3.2.1 Unitarity, perturbativity and vacuum stability

Besides experimental constraints, it is important to impose a minimal set of theoretical requirements which ensure that the results we will obtain are reliable and physically meaningful. To this goal we use the 2HDMC code [101]. First, we demand that all couplings be perturbative, which amounts to a condition $\left|\lambda_{i}\right|<4\pi$ for all couplings appearing in (1). Secondly, in any perturbative calculation the scattering matrix should be unitary order-by-order in perturbation theory. Failure of such a condition is typically associated with the development of strong dynamics which, again, renders a perturbative treatment unreliable. The tree-level scalar and vector scattering amplitudes are required to be smaller than $16\pi$, i.e. we allow that the unitarity limit be saturated already at tree-level. Lastly, the electroweak vacuum should be sufficiently long-lived. Here we impose the condition already implemented in 2HDMC which simply requires that the vacuum be completely stable. We note that metastable vacua in the IDM have been studied in [6, 7].

In our scanning procedure, parameter space points failing at least one of these requirements are immediately discarded, i.e. these theoretical requirements are imposed as “hard” constraints. In practice, we assign them a large enough negative log-likelihood pushing them well outside the $4\sigma$ region around the best-fit points.

3.2.2 Higgs invisible width

When $m_{H}^{0}<m_{h^{0}}/2$, the decay $h\to H^{0}H^{0}$ is allowed and contributes to the invisible decay width of the Higgs boson as

[TABLE]

In this region of parameter space, the coupling of a pair of $H^{0}$ to $h^{0}$ is constrained by LHC measurements which set an upper bound, ${\rm BR}_{\rm inv}\lesssim 0.23$ [41]. For the numerical analysis we use the log-likelihood function for ${\rm BR}_{\rm inv}$ provided in [41].

3.2.3 Electroweak precision observables

The new states which are present in the IDM can contribute to the $S$, $T$ and $U$ oblique parameters [102, 103]. Deviations from the Standard Model expectations in $U$ are negligible [17]. Hence we assume the latter to be zero and consider only $S$ and $T$. We compute their $\chi^{2}$ contribution through

[TABLE]

where $\boldsymbol{v}^{T}\equiv(S-\hat{S},T-\hat{T})$ and ${\cal C}$ is the covariance matrix. For $U=0$, the electroweak fit performed in [42] gives the values

[TABLE]

for a reference Higgs mass of $125$ GeV, while the covariance matrix is given by

[TABLE]

The contributions of the new scalar states to $S$ and $T$ are computed with 2HDMC [101].

3.2.4 LEP-II bounds on the masses of the heavy ${\cal{Z}}_{2}$-odd states

The masses of the heavier ${\cal{Z}}_{2}$-odd states can be constrained by translating the corresponding mass bounds from searches for charginos and neutralinos at LEP-II. The former were recast in [10], yielding a rough bound on the mass of the charged states:

[TABLE]

Limits on the mass of the heavier neutral state (which, we remind, in our case we take to be $A^{0}$) were extracted in [9] and amount to

[TABLE]

In the subsequent analysis, we will simply restrict our scanning regions to $m_{A^{0}},m_{H^{\pm}}\gtrsim 100~{}\mathrm{GeV}$ in order for these limits to be (conservatively) satisfied, without including them in our global fit analysis.

3.3 Scan settings

In order to perform the global fit we use MultiNest [34, 35], which allows an efficient scan of the parameter space under investigation. The considered parameters and respective scan ranges are summarized in table 2.

Although MultiNest is particularly suited for Bayesian analyses, in this work we will solely adopt a frequentist interpretation. This is possible provided that the posterior, and hence the resulting likelihood, has been explored in sufficient detail. This approach has two advantages. On the one hand, the derived constraints are not dependent on the priors chosen to explore the parameters. On the other hand, we can combine the output of different MultiNest scans, which allows an efficient use of the generated chains and provides the possibility to specifically improve on the coverage of the parameter space. However, in this approach the density of points, which in the Bayesian interpretation has a precise meaning (namely, it traces the posterior distribution), does not have any physical relevance anymore.

To ensure that the likelihood is sampled in enough detail, we run multiple MultiNest scans with high-accuracy settings, using between 600 and 3000 live points, a tolerance between $\texttt{tol}=$ 0.1–0.001, and an enlargement factor between $\texttt{efr}=0.3$–$0.6$ in order to achieve a good exploration of the tails of the distribution.

For the resulting fits we perform marginalization over parameters with the profile likelihood method [104] and draw contours at a certain confidence level following the expectation of a (two-dimensional) $\chi^{2}$ distribution.

4 Results and discussion

We now proceed to discuss our main results. As a first step, we update on the status of the IDM by performing a global fit including all constraints and observables described in section 3 except the GCE spectrum and the two unconfirmed dwarfs described in section 3.1.4 (Tucana III and IV). Subsequently, we include the GCE as well as the new dwarfs and check whether the IDM can provide an explanation for the GCE whilst satisfying all other constraints.

4.1 Global dark matter fit

The results of our global fit are presented in figure 1, where we show projections of the parameter space defined in table 2 onto 2-dimensional planes of all combinations of the involved parameters. As $\log(J/J_{\text{nom}})$ only concerns the GCE it is not included in the plot. Furthermore, we do not show $\lambda_{2}$ as it is virtually featureless, since the considered observables do not have any dependence on $\lambda_{2}$ at leading order in perturbation theory. We highlight 1-, 2-, 3- and 4$\sigma$ regions of the log-likelihood (brown, red, orange and yellow points respectively) around the best-fit point (represented by a white dot). The Higgs boson mass is fixed at $m_{h^{0}}=125.09$ GeV.

A first observation that can be made by inspecting figure 1 is that the masses of the ${\cal{Z}}_{2}$-odd scalars cannot vary arbitrarily with respect to each other. Concretely, the mass splitting between the lightest ($H^{0}$ in our case) and the next-to-lightest ${\cal{Z}}_{2}$-odd particle cannot exceed a few hundred GeV: it can reach maximally $\sim 500$ GeV for $m_{H^{0}}=100$ GeV and decreases to $\sim 300$ GeV for $m_{H^{0}}=500$ GeV. This behavior is a consequence of the perturbativity requirement: from the relations (3) we observe that $m_{H^{0}}^{2}-m_{A^{0}}^{2}=\lambda_{5}v^{2}$, implying that indeed large mass splittings between the two scalars can drive the $\lambda_{5}$ coupling to non-perturbative values. At the same time, the mass difference between $A^{0}$ and $H^{\pm}$ has to be small due to the constraints from the $S$, $T$ parameters described in section 3.2.3: generically, contributions to $S$ and $T$ are due to the breaking of the custodial symmetry of the scalar potential (1) which is induced by the $\lambda_{4,5}$ terms. Assuming all parameters to be real, the symmetry can be restored if and only if $\lambda_{4}=\lambda_{5}$ [105], whereas deviations from this condition amount to contributions to the oblique parameters. At the same time, from the relations (3) we see that $m_{A^{0}}^{2}-m_{H^{\pm}}^{2}=(\lambda_{4}-\lambda_{5})v^{2}/2$. This, in turn, implies that the contributions to $S$ and $T$ increase as the mass splitting between the two states becomes large.

Secondly, from figure 1 we see that the IDM can account for the entire DM abundance in the Universe ($R\sim 1$) in three distinct $m_{H^{0}}$ regions: one centered around half the Higgs mass (the so-called “funnel region”), one around 72 GeV and, finally, for relatively large $m_{H^{0}}\gtrsim 500$ GeV. The general reasons for this behavior have been analyzed in the literature [22, 23, 24, 25, 17]: ignoring, for the moment, direct detection constraints, the IDM can reproduce the observed DM abundance in the Universe in three $m_{H^{0}}$ regions. The first corresponds to $m_{H^{0}}<m_{W}$, where annihilation proceeds through the Higgs portal - type process $H^{0}H^{0}\rightarrow h^{0}\rightarrow f\bar{f}/VV^{*}$ (with $V=W^{\pm}/Z$) as well as through direct annihilation via the point-like $H^{0}$-$H^{0}$-$V$-$V$ vertex. Annihilation into virtual gauge bosons increases in importance as the corresponding kinematic threshold is approached from below. Besides, the LEP constraints on the heavier ${\cal{Z}}_{2}$-odd scalar masses described in section 3.2.4 exclude the possibility of coannihilation in this mass range. Once $m_{H^{0}}$ becomes larger than $m_{W}$, annihilation into gauge bosons becomes dominant. If fact, it is too efficient, so destructive interference must occur between the Higgs portal - like diagram and the one involving the four-vertex, which for a Higgs mass of $\sim 125$ GeV can happen for negative (positive) values of $\lambda_{L}$ if $m_{H^{0}}>m_{h^{0}}/2$ ($m_{H^{0}}<m_{h^{0}}/2$). Above roughly 120 GeV, this interference cannot be efficient enough so that DM becomes necessarily underabundant (this upper value depends on the Higgs boson mass [25]).

The Imposition of direct detection constraints wipes out most of this parameter space, since the DM-nucleon scattering cross section is proportional to $\lambda_{L}^{2}$. The only regions surviving are those characterized by very small values of $\lambda_{L}$, which correspond either to the Higgs funnel region or to the regime where annihilation into pairs of virtual gauge bosons becomes efficient enough (but not too efficient so that cancellations are needed) without requiring large values of $\lambda_{L}$.

In the high-mass regime, additional effects come into play. As described in detail in [23], once $m_{H^{0}}\gtrsim 500$ GeV the destructive interference between the annihilation diagram involving the quartic $H^{0}$-$H^{0}$-$V$-$V$ coupling and $t$-channel diagrams involving the heavier ${\cal{Z}}_{2}$-odd particles can become efficient enough so as to bring the $H^{0}$ self-annihilation cross section down to acceptable levels. The cancellation becomes more exact the smaller the mass splitting between the inert states, and in practice the condition $\left|m_{H^{0}}-m_{A^{0}/H^{\pm}}\right|\lesssim 10$ GeV must be satisfied otherwise $\left\langle\sigma v\right\rangle$ is too large (this is actually also dictated by unitarity arguments). The predicted DM abundance then depends on the interplay of this cancellation mechanism with the contributions from $s$-channel exchange of a Higgs boson and coannihilation. If the cancellation is exact, large enough values of $\lambda_{L}$ (of ${\cal{O}}(0.3)$ [17]) are needed in order to saturate the Planck bound. Besides, since in this regime DM is relatively heavy, it is only poorly constrained by direct/indirect detection and/or collider searches.

Dropping the $R\sim 1$ requirement, i.e. going to scenarios in which the IDM only accounts for a fraction of the total DM content of the Universe, allows for more freedom in both $m_{H^{0}}$ and $\lambda_{L}$. This effect is particularly pronounced in the $H^{0}$ mass range between roughly 75 and 500 GeV, where DM naturally tends to be underabundant due to the efficiency of the direct $H^{0}$-$H^{0}$-$V$-$V$ coupling. In Higgs-portal models, it is typically the same coupling that controls the annihilation and the WIMP-nucleon scattering cross section. This implies that in such scenarios underabundant WIMP dark matter, $R<1$, does not amount to weaker direct detection bounds, since to first order both the relic density and direct detection constraints induce the same relation between the Higgs portal coupling and the WIMP fraction, $\lambda\propto R^{-1/2}$. In the IDM this is true only for regions where Higgs-mediated annihilation is clearly dominant. Above $\sim 72\,$GeV the direct $H^{0}$-$H^{0}$-$V$-$V$ interaction dominates the annihilation while it does not enter the direct detection cross section at leading order. This introduces more freedom in $m_{H^{0}}$ and $\lambda_{L}$ for $R<1$. Note that this is true only for the current direct detection constraints that have not yet reached sensitivity to the loop-induced electroweak corrections [84, 85]. For future direct detection experiments this situation can change, cf. the discussion in section 4.3.

With direct detection constraints being largely inefficient for DM masses above $\sim 75$ GeV, the behaviour of the allowed values of the relic abundance as a function of $m_{H^{0}}$ is largely determined by the interference pattern between the direct $H^{0}$-$H^{0}$-$V$-$V$ coupling and $t/u$-channel diagrams involving the heavier ${\cal{Z}}_{2}$-odd particles [23]. In particular, the largest allowed values of $R$ correspond to the smallest mass splittings between the Inert Doublet scalars. This might appear counter-intuitive, since from a Boltzmann suppression standpoint this is the regime where coannihilation effects should become the most relevant. However, this is also the regime in which the quartic couplings in the scalar potential (1) vanish. In this limit, the cancellation between four-vertex interactions and $t/u$-channel diagrams involving the heavier scalars becomes maximal, with a similar remark also applying to coannihilation processes. Turning on the quartic couplings, i.e. increasing the mass splitting between the Inert Doublet scalars, can only increase the total $H^{0}$ (co-)annihilation cross section and, hence, decrease the predicted DM abundance. We thus see that in the underabundant IDM region the upper limit on $R$ corresponds to quasi-degenerate Inert Scalars, whereas the lower one to large mass splittings within the Inert Doublet (as well as to large values of $\lambda_{L}$).

4.2 Fitting a possible signal from dark matter annihilation

We now examine whether the IDM can accommodate the GCE as a signal from DM annihilation. In addition to the observables considered in section 4.1 we include the GCE likelihood in our global fit as described in section 3.1.4. The result of the fit is shown in figure 2 (again, we omit $\lambda_{2}$ since it does not enter any of the observables we consider). We find two distinct regions in parameter space in which the IDM can explain the GCE. The respective best-fit points are summarized in table 3 and the corresponding spectra are shown in figure 3.

Higgs funnel region

The first region lies close to the Higgs funnel, where DM annihilation proceeds predominantly via $s$-channel Higgs exchange near the resonance, $m_{H^{0}}\simeq m_{h^{0}}/2$. In figure 2b we zoom into this part of the parameter space, restricting $\lambda_{L}$ to positive values in order to allow for a logarithmic scaling (we will comment on the asymmetry regarding positive and negative $\lambda_{L}$ below). Upon closer inspection we can see that this region splits up into two subregions. In the first one (hereafter referred to as region 1a), $m_{H^{0}}$ is restricted to an extremely narrow range between 62.5 GeV and 62.55 GeV (i.e. just below $m_{h^{0}}/2$) and $\lambda_{L}$ to small values $10^{-4}\lesssim\lambda_{L}\lesssim 10^{-2}$. In this case the IDM can also account for the entire DM abundance in the Universe, $R\simeq 1$. In fact, large $R$ is slightly preferred as we will explain further below. In the second subregion (region 1b), $m_{H^{0}}$ is slightly larger, $\lambda_{L}\sim 10^{-2}$ and the fit favors $R\lesssim 0.2$.

This structure is mainly driven by the interplay of two observables, namely the DM relic density and the GCE itself. In the vicinity of the resonance the annihilation cross section exhibits a large velocity dependence. Hence, a small variation in $m_{H^{0}}$ can alter the ratio $\langle\sigma v\rangle_{\text{today}}/\langle\sigma v\rangle_{\text{freeze-out}}$ by orders of magnitude. In particular, for $m_{H^{0}}$ just above $m_{h^{0}}/2$, this ratio can become larger than $10^{2}$. Moreover, the gamma-ray flux and the predicted relic DM density scale differently with the faction $R$: the relic density scales as $\Omega h^{2}|_{\text{DM,\,total}}\propto 1/(R\left\langle\sigma v\right\rangle_{\rm freeze-out})$ whereas the gamma-ray flux is proportional to $\left\langle\sigma v\right\rangle_{\rm today}\times R^{2}$. Now, the likelihood function is minimised for $\Omega h^{2}|_{\text{DM,\,total}}\sim\Omega h^{2}|_{\text{Planck}}$, which in turn implies $\left\langle\sigma v\right\rangle_{\rm freeze-out}\times R\sim\left\langle\sigma v\right\rangle_{\rm thermal}$. At the same time, in order to reproduce the GCE (which can be explained by an annihilation cross section of the order of the thermal one) a flux corresponding to $\langle\sigma v\rangle_{\rm thermal}\sim\left\langle\sigma v\right\rangle_{\rm today}\times R^{2}$ is required. Hence, our fit globally prefers parameter space regions in which

[TABLE]

As $R\leq 1$, this condition can only be satisfied if $\left\langle\sigma v\right\rangle_{\rm today}\geq\left\langle\sigma v\right\rangle_{\rm freeze-out}$, a situation occurring for $m_{H^{0}}\gtrsim m_{h^{0}}/2$.

In the best-fit region 1a $m_{H^{0}}$ is finely tuned to a value very close to (but slightly smaller than) $m_{h^{0}}/2$. The minimal mass in this region corresponds roughly to $\langle\sigma v\rangle_{\text{today}}/\langle\sigma v\rangle_{\text{freeze-out}}\simeq 1$ and, hence, $R=1$. As $m_{H^{0}}$ approaches $m_{h^{0}}/2$ from below, $\langle\sigma v\rangle_{\text{today}}/\langle\sigma v\rangle_{\text{freeze-out}}$ increases extremely rapidly for very small variations of $m_{H^{0}}$ and, hence, the gamma-ray flux remains large enough despite the reduction in the overall DM density as explained before. Region 1b sets in for slightly larger masses just above $m_{h^{0}}/2$, where $\langle\sigma v\rangle_{\text{today}}/\langle\sigma v\rangle_{\text{freeze-out}}$ is again large but gradually starts to decrease. However, this region does not extent up to masses for which $\langle\sigma v\rangle_{\text{today}}/\langle\sigma v\rangle_{\text{freeze-out}}\simeq 1$ due to an interplay of two aspects. On the one hand, for increasing mass the shape of the gamma-ray spectrum (mainly due to annihilation into $b\bar{b}$) tends to yield a worse fit of the GCE. On the other hand, the larger $\lambda_{L}$ values required start to be in tension with limits from direct detection.

This behavior is similar to the one found for the singlet scalar Higgs portal model in [88]. However, there are some differences worth commenting upon. In the singlet scalar model, around the Higgs funnel region DM annihilation only occurs via $s$-channel Higgs exchange. Hence, the relative contribution to annihilation does not depend on the coupling between the DM particle and the Higgs. In the IDM the additional four-vertex interaction $H^{0}$-$H^{0}$-$V$-$V$ interferes with $H^{0}H^{0}\to h\to VV^{(*)}$, which introduces some additional features.

In region 1a $\lambda_{L}$ becomes as small as $10^{-4}$ for $R\simeq 1$. For such small values of $\lambda_{L}$ the contributions of the four-vertex and Higgs exchange diagrams are of the same order and, hence, there is a strong destructive (constructive) interference for positive (negative) $\lambda_{L}$ in the annihilation cross section into vector bosons.888We remind that the interference is destructive (constructive) for positive (negative) $\lambda_{L}$ when $s<(m_{h^{0}}/2)^{2}$ which, in the zero-velocity limit, amounts to $m_{H^{0}}<m_{h^{0}}/2$. The situation is inversed for $s>(m_{h^{0}}/2)^{2}$ [25]. This leads to a large variation of the $WW^{*}$ contribution to the total gamma-ray spectrum, ranging from a negligible fraction for small positive $\lambda_{L}$ to around 50% for small negative $\lambda_{L}$. Since the $b\bar{b}$ annihilation spectrum provides a better fit of the GCE than the $WW^{*}$ one, the fit tends to favor a small positive $\lambda_{L}\simeq 10^{-4}$ coupling, where the $WW^{*}$ contribution is below the percent level. This is the origin of the asymmetry in $\lambda_{L}$ and furthermore explains the preference for large $R$ values in region 1a – due to the strong correlation between $\lambda_{L}$ and $R$. For larger values of $\lambda_{L}$ the Higgs-exchange diagram dominates and the $WW^{*}$ contribution becomes comparable to the one in the singlet scalar model.

In region 1b the situation is reversed and we obtain a smaller $WW^{*}$ contribution for negative values of $\lambda_{L}$. Although, following the same line of reasoning, the GCE should now favor negative $\lambda_{L}$ values, the overall likelihood is actually maximized for positive ones. This is caused by the DM relic density contribution to the total likelihood. Due to the strong velocity-dependence the various annihilation channels behave very different during freeze-out and at present times and in this case interference effects become important for larger couplings, thus giving rise to the required annihilation cross section for positive $\lambda_{L}$ only. This causes a strong preference for positive $\lambda_{L}$ in region 1b, despite the slightly worse GCE fit (given the larger $WW^{*}$ contribution). These remarks also explain the overall preference for region 1a in the fit, contrary to the case of the singlet scalar Higgs portal model analysed in [88].

Region around 72 GeV

The second mass region where the IDM can explain the GCE lies around $m_{H^{0}}=72$ GeV. In this case the fit to the GCE is considerably worse than in regions 1, $\chi^{2}/\text{dof}=33.6/22$. However, it globally remains at acceptable levels (within 2$\sigma$ from the best-fit point of region 1a), which is partly driven by the other observables considered, cf. table 3. Moreover, this region does not require a large degree of fine-tuning of $\lambda_{L}$ and $m_{H^{0}}$ and strongly favors that the IDM accommodates the full observed DM abundance ($R=1$) through sufficiently efficient annihilation of $H^{0}$ pairs into virtual gauge bosons via the quartic $H^{0}$-$H^{0}$-$V$-$V$ interaction. The corresponding coupling is a pure gauge coupling, completely independent of $\lambda_{L}$, and the latter can, hence, be tuned to sufficiently small values in order to evade current direct detection constraints. As the DM mass increases, and $m_{H^{0}}$ approaches the $WW$ threshold, this process gradually becomes too efficient which would rather imply $R<1$, cf. the discussion in section 4.1. However, the spectrum for annihilation into $WW^{*}$ provides a considerably worse fit to the GCE for masses above $m_{H^{0}}=72\,$GeV. (In fact, for pure $WW^{*}$ annihilation a DM mass around $55\,$GeV would fit best.) Note that $\langle\sigma v\rangle_{\text{today}}$ is somewhat smaller than $\langle\sigma v\rangle_{\text{freeze-out}}$ due to the kinematic suppression away from the $WW$ threshold for small velocities. As a consequence the IDM tends to undershoot the required flux for fitting the GCE which is reflected in the preference for positive $\log(J/J_{\text{nom}})$. The gamma-ray spectra resulting from the best-fit points of regions 1a, 1b and 2 are shown for comparison in figure 3, along with the corresponding Fermi-LAT measurements.

In our fit so far we only considered confirmed dSph targets, for which the $J$-factor has been measured. In four of the recently discovered dSph targets, slight excesses (each $\sim 2\sigma$ local) have been found: Reticulum II, Tucana III, Tucana IV and Indus II [68, 39]. It is therefor interesting to see in how far these excesses are compatible with the DM explanation of the GCE. To this end, we include the likelihood contribution of Tucana III and IV in our fit. The result is shown in figure 4 where we omit the GCE likelihood contribution in order to allow for a comparison. We indeed find the same two regions as before. Interestingly, the overall best-fit point now lies in region 2. However, given that the excess is mild the preference for the best-fit point with respect to, e.g., the high mass region, $m_{H^{0}}\gtrsim 600\,\text{GeV}$ (that does not fit the excess), is only at the level of 1–2$\sigma$.

Let us also note that recently a hint for a possible DM annihilation signal has also been found in the AMS-02 antiproton data, pointing to a DM mass around 60–80 GeV and a cross section around $3\times 10^{-26}\text{cm}^{2}\text{s}^{-1}$ [71]. Although for a given annihilation channel this excess favors a somewhat larger cross section than the GCE, these two observations are still compatible with each other taking into account the uncertainties on the local DM density [75], in particular for the annihilation channels $b\bar{b}$ and $WW$ which are important in the IDM.

4.3 Future prospects for direct, indirect detection and the LHC

Finally, let us briefly discuss the future prospects for direct detection, indirect detection and collider searches in the IDM, focusing in particular on the GCE best-fit regions found in section 4.2.

Direct detection

In the near future the sensitivity of direct detection experiments is expected to improve even further with the advent of new experiments. The expected 90% CL exclusion limit on the spin-independent DM-nucleon scattering cross section for a DM mass between 62 and 72 GeV is estimated at $(1.7\!-\!1.8)\times 10^{-47}\,\text{cm}^{2}$ for Xenon1T [106], $(2.4\!-\!2.5)\times 10^{-48}\,\text{cm}^{2}$ for LZ [107] and $(2.8\!-\!3.0)\times 10^{-49}\,\text{cm}^{2}$ for Darwin [108]. The latter improves by up to three orders of magnitude upon the sensitivity of LUX 2016 [37] which we took into account in our fitting procedure and which provides a limit in the vicinity of $(2.2\!-\!2.3)\times 10^{-46}\,\text{cm}^{2}$ for the same mass range. The three GCE best-fit regions found in section 4.2 are characterized by small $\lambda_{L}$ and, hence, a small DM-nucleon scattering cross section at tree-level. Especially for regions 1a and 2, our findings show that the tree-level scattering cross section is too small to be challenged by upcoming experiments. However, as mentioned in section 3.1.2, electroweak radiative corrections provide a contribution to the DM-nucleon scattering that is independent of $\lambda_{L}$, and which can actually dominate over the tree-level contribution. For the best-fit points in region 1a and 2 the radiative corrections are of the order of $1.5\times 10^{-48}\,\text{cm}^{2}$ [84] which brings them slighly outside the projected sensitivity of LZ but well within the reach of Darwin.

Note that for points with $R<1$ within the same mass range, the electroweak corrections lose importance with respect to the tree-level contribution. Such points correspond to larger values of $\lambda_{L}$, for which the decrease in the signal due to the small value of $R$ is roughly compensated by the larger $\lambda_{L}$-induced scattering cross section, which is not true for the electroweak contribution.

Besides, upcoming direct detection experiments can also provide handles for the high-mass region of the IDM. As discussed in section 4.1, in the range between 100 and 600 GeV where the IDM tends to provide too little DM, the largest $R$ values are obtained for small mass splittings between $H^{0}$ and the heavier inert doublet partners, as well as for small values of $\lambda_{L}$. This is true in particular for the points with the smallest masses that allow for $R=1$ around 600 GeV. For these points, the contribution from electroweak corrections is likely to dominate the tree-level WIMP-nucleon scattering cross section since $\lambda_{L}$ can be made sufficiently small so as to render the Higgs-mediated DM-nucleon scattering negligible. The corresponding cross section (in the limit $\lambda_{L}=0$) is about $1.4\times 10^{-47}\,\text{cm}^{2}$ [84] for a DM mass of 600 GeV. This value is just within the projected reach of LZ [107], which will, hence, start to push the lower mass limit of the IDM high-mass region (assuming $R=1$).

Indirect detection

On the side of indirect detection, Fermi-LAT searches for gamma-rays in dwarf spheroidal galaxies will ultimately probe a possible explanation of the GCE within the IDM. To illustrate this point, in figure 5 we project our scan points onto the $m_{H^{0}}$-$\langle\sigma v\rangle R^{2}$ plane and superimpose them with the projected sensitivity of Fermi-LAT, assuming 15 years of data acquisition [109]. One subtlety concerns the fact that these projections are derived assuming annihilation into a pure $b\bar{b}$ final state. In the two IDM GCE best-fit regions that we have found the dominant DM annihilation channels are either $b\bar{b}$ or $WW^{*}$, while for the high mass region the $WW$ and $ZZ$ channels dominate. Since, however, the current bounds on $\langle\sigma v\rangle$ from dwarf spheroidal galaxies are similar for the $b\bar{b}$ and $WW$ channels (the differences in the derived limits being of the order of 25%), these projections can, indeed, provide an adequate estimate of the indirect detection prospects of the IDM.

The result indicates that the entire 1–2$\sigma$ region that is compatible with a DM interpretation of the GCE can be probed with 15 years of LAT data acquisition.

For the high mass region CTA [110] is expected to provide better sensitivity. In the left panel of figure 5 we show the projected limits for 100hr of CTA observation of the Galactic Center, taken from [111], assuming annihilation into $WW$ and an Einasto DM density profile. Note that for a more cuspy generalized NFW profile (as it is found to be compatible with the GCE) those limits are stronger by up to an order of magnitude [111].

Collider searches

On the side of collider searches, given the very small $\lambda_{L}$ values that characterise the GCE best-fit regions, the standard LHC mono-X searches for $H^{0}$ pairs become inefficient. Of more interest are searches for the heavier ${\cal{Z}}_{2}$-odd states subsequently decaying into visible products and missing energy. As pointed out in [18], the production of $A^{0}H^{0}$ and $H^{\pm}H^{\pm}$ pairs does not depend on $\lambda_{L}$, as it only involves gauge couplings and could, at least for sufficiently light $A^{0}$ and $H^{\pm}$ masses, give rise to visible signals in the dileptons + MET channel. Larger mass splittings, and in particular once $m_{A^{0}}-m_{H^{0}}>m_{Z}$, are harder to probe due to the kinematic cuts imposed by the corresponding searches in order to elliminate the dominant $Z$-induced background.

Another interesting channel that was suggested in [2] concerns the high-mass region of the IDM and in particular cases where the mass splitting between $H^{0}$ and $H^{\pm}$ is small. In this case, CMS searches for disappearing charged tracks at 8 TeV [112] already exclude a range of $H^{0}$ masses between 490 and 550 GeV for $m_{H^{\pm}}-m_{H^{0}}\sim 0.2$ GeV regardless of the value of $\lambda_{L}$. It would be interesting to follow the evolution of these constraints, since they could provide a rather unique handle on the IDM high-mass regime, especially in cases of small $\lambda_{L}$ values.

5 Conclusion

In this paper we presented a global fit of the IDM taking into account state-of-the-art constraints from collider observables, direct and indirect DM searches as well as theoretical considerations. We performed a detailed exploration of the IDM parameter space and updated upon existing studies that have shown that the low-mass regime of the IDM is by now very efficiently constrained, in particular once direct detection bounds are combined with the LHC Higgs mass measurement and the condition that the IDM reproduces the total DM abundance in the Universe. Going a step further, we relaxed the latter requirement and instead demanded simply that the Universe does not get overclosed assuming a standard thermal history, which allowed us to examine the substantial – and much less frequently studied – regions of parameter space that open up.

We then examined whether the IDM can accommodate the excessive Galactic bulge emission that has been reported by numerous groups in the Fermi-LAT data (“Galactic center excess”). We found that this is indeed the case, in two distinct regions of parameter space: the first lies around the so-called “Higgs funnel”, $m_{H^{0}}\sim m_{h^{0}}/2$, in which $H^{0}$ particles annihilate mostly into $b\bar{b}$ pairs through a quasi-on-shell Higgs boson. Interestingly, the strong dependence of the total thermally averaged self-annihilation cross section on the DM velocity in this region of parameter space makes it possible to explain the GCE even if the IDM only accounts for a small fraction of the DM abundance in the Universe, as the cross section computed at velocities relevant for indirect detection can supersede the corresponding ones computed at freeze-out velocities by several orders of magnitude. The second $H^{0}$ mass range in which the IDM can explain the GCE lies around 72 GeV, close (but not too close) to the $WW$ threshold. In this case, DM annihilates predominantly into pairs of virtual $W$’s via the quartic $H^{0}$-$H^{0}$-$W^{+}$-$W^{-}$ coupling that appears in the Lagrangian gauge kinetic terms for the inert doublet. In this case the dependence of $\left\langle\sigma v\right\rangle$ on the DM velocity is milder, and the GCE can mostly be explained for $R\equiv\Omega h^{2}|_{\text{IDM}}/\Omega h^{2}|_{\text{DM,\,Planck}}\sim 1$, which is attainable without conflicting current direct DM searches since the corresponding coupling only enters the WIMP-nucleon scattering cross section at next-to-leading order. To the best of our knowledge, the existence of this region constitutes a novelty of the IDM with respect to simpler “portal” models and has not been pointed out before.

Both of these regions of parameter space involve small (${\cal{O}}(10^{-2})$ or less) values of the coupling between $H^{0}$ pairs and the standard model Higgs boson and masses $m_{H^{0}}\gtrsim m_{h^{0}}/2$. As a consequence, collider probes of the IDM GCE explanation do not appear to be particularly promising. On the other hand, our findings show that Fermi-LAT constraints from searches for DM in dwarf spheroidal galaxies should confirm or exclude this scenario with 15 years of data acquisition. At the same time, the next generation of direct detection experiments will also provide complementary information in this direction, since the LZ and especially the Darwin experiments will reach a level of sensitivity that will enable them to probe the electroweak radiative corrections to the WIMP-nucleon scattering cross section, which are fixed by the electroweak coupling strength and do not depend, hence, on the additional parameters of the model. Besides, both direct and indirect detection experiments are expected to probe a substantial fraction of the high-mass regime of the IDM, at least for scenarios that saturate the Planck bound on the DM abundance in the Universe.

Acknowledgements

We would like to thank Alessandro Cuoco and Michael Krämer for very helpful discussions. B.E. and J.H. acknowledge support by the German Research Foundation DFG through the research unit “New physics at the LHC”.