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

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.


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 selfannihilation 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 wellmotivated theoretical models is hence an important task in order to pinpoint the nature of DM. The inert doublet model (IDM) 1 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 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 as regards 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 as regards 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 testing 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 DM 100 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 ∼ 2σ local significance) have been found in four of the recently discovered dSph targets [39,67,68] 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. 2 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 Sect. 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 Sect. 3. Section 4 contains our main results and a discussion of the future sensitivities of upcoming experiments to the best-fit parameter regions. We conclude in Sect. 5.

The inert doublet model
The IDM is a special case of a two-Higgs doublet model in which an exact 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 is taken to be odd. With respect to the standard model Lagrangian, the only modifications consist of the introduction of gauge kinetic terms for and an additional piece in the scalar potential, which in the IDM reads Upon electroweak symmetry breaking, the two scalar doublets can be expanded in component fields as where v 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 where, following common conventions, we have defined 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: which, using Eqs. (3), (4) and the potential minimization condition ∂ V /∂h 0 h 0 =0 = 0 can be traded for the physically more intuitive set of parameters with m h 0 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 Z 2 -odd particle is neutral (m H 0 /A 0 < m H ± ), 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.

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 analyzed 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 λ 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 analyzed 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.

Relic density
Assuming a standard cosmological history allows us to link the relic H 0 density from thermal freeze-out, h 2 | IDM , to the DM density measured by Planck, h 2 | Planck = 0.1198 ± 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 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 dif-ferent scales which are relevant for the various DM observables considered here. Then the total DM density is given by The H 0 relic density h 2 | 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 χ 2 contribution for the relic density via 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 σ 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 up to an irrelevant constant.

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 where f N ∼ 0.30 [81] denotes the strength of the effective Higgs-nucleon interaction, and 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 (massdependent) 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 loglikelihood 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 percent level in the mass region of interest.
It should be noted that as (11) is proportional to λ 2 L , for very small λ 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 λ L and below ∼ 2 × 10 −47 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 Sect. 4.3.

Indirect detection constraints: dwarf spheroidal galaxies
In the low-mass region of the IDM where annihilation occurs predominantly into bb and/or W + W − pairs, the most stringent limits on the velocity-averaged annihilation cross section, σ v , 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 × flux in an energy bin between E min and E max is where dN γ /dE γ is the differential photon spectrum per annihilation and m DM is the mass of the DM particle. The integral of the DM density, ρ 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 the 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 W W * , Z Z * , 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 ∼ 2σ local) have been found: specifically in Reticulum II, Tucana III, Tucana IV and Indus II [39,67,68] (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 DMdominated 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 distancebased 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. 3 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.
In the E 2 × 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, σ v ∼ 10 −26 cm 3 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 Sect. 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 3 For instance, for the benchmark point IV in [26] (m h = 120 GeV, m 0 H = 70 GeV and λ L −0.07) the predicted line signal, σ v γ γ = 7.6 × 10 −30 cm 3 s −1 , falls below the upper limit, σ v UL γ γ = 5.2 × 10 29 cm 3 s −1 [92], by almost an order of magnitude while the cross section into bb alone, σ v bb = 1.6 × 10 26 cm 3 s −1 , is already relatively close to the respective limit from dwarfs, σ v UL bb = 2.6 × 10 26 cm 3 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 γ = 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. 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 where the ROI is a 40 • × 40 • region centered on the Galactic center with a stripe of ±2 • 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 γ = 1.2 ± 0.08 [52], 4 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 log Note that since the authors of [52] normalize the GCE flux dividing by the angular size of the analyzed region, we divide J 40 • by the corresponding solid angle, = 0.43 sr.

Likelihood for the GCE signal
We compute the χ 2 contribution for the GCE-including the contribution from J 40 • -by 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 • . i j is the covariance matrix given in [52] and log J 40 • , nom and σ logJ 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 spec-tra for individual annihilation channels(including the 3-body final states W W * , Z Z * ) 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 δ i j (σ rel t i ) 2 representing a diagonal error equal to a fraction σ rel of the model prediction. We choose σ rel = 10% as discussed in [88,99,100]. Up to an irrelevant constant factor, the resulting log-likelihood is where 3.2 Non-dark matter observables

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 |λ i | < 4π 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π , i.e. we allow that the unitarity limit be saturated already at tree level.
Lastly, the electroweak vacuum should be sufficiently longlived. 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σ region around the best-fit points.

Higgs invisible width
When m 0 H < m h 0 /2, the decay h → H 0 H 0 is allowed and contributes to the invisible decay width of the Higgs boson as 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, BR inv 0.23 [41]. For the numerical analysis we use the log-likelihood function for BR inv provided in [41].

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 χ 2 contribution through where v T ≡ (S −Ŝ, T −T ) and C is the covariance matrix. For U = 0, the electroweak fit performed in [42] gives the valueŝ for a reference Higgs mass of 125 GeV, while the covariance matrix is given by The contributions of the new scalar states to S and T are computed with 2HDMC [101].

LEP-II bounds on the masses of the heavy Z 2 -odd states
The masses of the heavier 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: Limits on the mass of the heavier neutral state (which, we recall, in our case we take to be A 0 ) were extracted in [9] and amount to In the subsequent analysis, we will simply restrict our scanning regions to m A 0 , m H ± 100 GeV in order for these limits to be (conservatively) satisfied, without including them in our global fit analysis.

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 for 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 tol = 0.1-0.001, and an enlargement factor between 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) χ 2 distribution.

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 Sect. 3 except the GCE spectrum and the two unconfirmed dwarfs described in Sect. 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.

Global dark matter fit
The results of our global fit are presented in Fig. 1, where we show projections of the parameter space defined in Table 2 Fig. 1 Global fit results for the IDM free parameters m H 0 , m 0 A , m ± H , λ L and the DM fraction R. The brown, red, orange and yellow points lie within 1-, 2-, 3-and 4σ away from the best-fit point (denoted by a white dot), respectively. Here we take into account the log-likelihood contributions from all observables described in Sect. 3, except the GCE spectrum and unconfirmed dwarfs onto two-dimensional planes of all combinations of the involved parameters. As log(J/J nom ) only concerns the GCE it is not included in the plot. Furthermore, we do not show λ 2 as it is virtually featureless, since the considered observables do not have any dependence on λ 2 at leading order in perturbation theory. We highlight 1-, 2-, 3-and 4σ 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 Fig.  1 is that the masses of the Z 2 -odd scalars cannot vary arbitrarily with respect to each other. Concretely, the mass split-ting between the lightest (H 0 in our case) and the next-tolightest Z 2 -odd particle cannot exceed a few hundred GeV: it can reach maximally ∼ 500 GeV for m H 0 = 100 GeV and decreases to ∼ 300 GeV for m H 0 = 500 GeV. This behavior is a consequence of the perturbativity requirement: from Eqs.
(3) we observe that m 2 H 0 −m 2 A 0 = λ 5 v 2 , implying that indeed large mass splittings between the two scalars can drive the λ 5 coupling to non-perturbative values. At the same time, the mass difference between A 0 and H ± has to be small due to the constraints from the S, T parameters described in Sect. 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 λ 4,5 terms. Assuming all parameters to be real, the symmetry can be restored if and only if λ 4 = λ 5 [105], whereas deviations from this condition amount to contributions to the oblique parameters. At the same time, from Eqs. (3) we see that m 2 A 0 − m 2 H ± = (λ 4 − λ 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 Fig. 1 we see that the IDM can account for the entire DM abundance in the Universe (R ∼ 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 500 GeV. The general reasons for this behavior have been analyzed in the literature [17,[22][23][24][25]: 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 → h 0 → ff /V V * (with V = W ± /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 Z 2 -odd scalar masses described in Sect. 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 ∼ 125 GeV can happen for negative (positive) values of λ 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 λ 2 L . The only regions surviving are those characterized by very small values of λ 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 λ L .
In the high-mass regime, additional effects come into play. As described in detail in [23], once m H 0 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 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 m H 0 − m A 0 /H ± 10 GeV must be satisfied otherwise σ v 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 λ L (of 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 ∼ 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 λ 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 the relic density and the direct detection constraints induce the same relation between the Higgs-portal coupling and the WIMP fraction, λ ∝ R −1/2 . In the IDM this is true only for regions where Higgs-mediated annihilation is clearly dominant. Above ∼ 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 λ 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 Sect. 4.3.
With direct detection constraints being largely inefficient for DM masses above ∼ 75 GeV, the behavior 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 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 counterintuitive, 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 λ L ).

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 Sect. 4.1 we include the GCE likelihood in our global fit as described in Sect. 3.1.3. The result of the fit is shown in Fig. 2 (again, we omit λ 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 Fig. 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 m h 0 /2. In Fig. 2b we zoom into this part of the parameter space, restricting λ L to positive values in order to allow for a logarithmic scaling (we will comment on the asymmetry regarding positive and negative λ 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 and 62.55 GeV (i.e. just below Table 3 Parameters describing the GCE best-fit regions appearing in Fig. 2. The regions 1a and 1b correspond to the Higgs funnel with m H 0 ≈ 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  . 3 Comparison of the gamma-ray spectra predicted by the IDM best-fit points with the GCE spectrum m h 0 /2) and λ L to small values 10 −4 λ L 10 −2 . In this case the IDM can also account for the entire DM abundance in the Universe, R 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, λ L ∼ 10 −2 and the fit favors R 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 σ v today / σ v 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 h 2 | DM, total ∝ 1/(R σ v freeze−out ) whereas the gamma-ray flux is proportional to σ v today × R 2 . Now, the likelihood function is minimized for h 2 | DM, total ∼ h 2 | Planck , which in turn implies σ v freeze−out × R ∼ σ v 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 σ v thermal ∼ σ v today × R 2 is required. Hence, our fit globally prefers parameter space regions in which As R ≤ 1, this condition can only be satisfied if σ v today ≥ σ v freeze−out , a situation occurring for m H 0 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 σ v today / σ v freeze-out 1 and, hence, R = 1. As m H 0 approaches m h 0 /2 from below, σ v today / σ v 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 σ v today / σ v freeze-out is again large but gradually starts to decrease. However, this region does not extent up to masses for which σ v today / σ v freeze-out 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 bb) tends to yield a worse fit of the GCE. On the other hand, the larger λ 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 → h → V V ( * ) , which introduces some additional features.
In region 1a λ L becomes as small as 10 −4 for R 1. For such small values of λ 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) λ L in the annihilation cross section into vector bosons. 5 This leads to a large variation of the W W * contribution to the total gamma-ray spectrum, ranging from a negligible fraction for small positive λ L to around 50% for small negative λ L . Since the bb annihilation spectrum pro-vides a better fit of the GCE than the W W * one, the fit tends to favor a small positive λ L 10 −4 coupling, where the W W * contribution is below the percent level. This is the origin of the asymmetry in λ L and furthermore explains the preference for large R values in region 1a-due to the strong correlation between λ L and R. For larger values of λ L the Higgsexchange diagram dominates and the W W * contribution becomes comparable to the one in the singlet scalar model.
In region 1b the situation is reversed and we obtain a smaller W W * contribution for negative values of λ L . Although, following the same line of reasoning, the GCE should now favor negative λ 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 λ L only. This causes a strong preference for positive λ L in region 1b, despite the slightly worse GCE fit (given the larger W W * 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 analyzed 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, χ 2 /dof = 33.6/22. However, it globally remains at acceptable levels (within 2σ 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 λ 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 λ 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 W W threshold, this process gradually becomes too efficient which would rather imply R < 1, cf. the discussion in Sect. 4.1. However, the spectrum for annihilation into W W * provides a considerably worse fit to the GCE for masses above m H 0 = 72 GeV. (In fact, for pure W W * annihilation a DM mass around 55 GeV would fit best.) Note that σ v today is somewhat smaller than σ v freeze-out due to the kinematic suppression away from the W W 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 nom ). The gamma-ray spec- 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 ∼ 2σ local) have been found: Reticulum II, Tucana III, Tucana IV and Indus II [39,68]. 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 Fig. 4 where we omit the GCE likelihood contribution in order to allow for 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 600 GeV (that does not fit the excess), is only at the level of 1-2σ .
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 × 10 −26 cm 2 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 bb and W W , which are important in the IDM.

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.  [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) × 10 −46 cm 2 for the same mass range. The three GCE best-fit regions found in Sect. 4.2 are characterized by small λ 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 Sect. 3.1.2, electroweak radiative corrections provide a contribution to the DM-nucleon scattering that is independent of λ 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 × 10 −48 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 λ L , for which the decrease in the signal due to the small value of R is roughly compensated by the larger λ 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 Sect. 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 λ 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 λ L can be made sufficiently small so as to render the Higgs-mediated DM-nucleon scattering negligible. The corresponding cross section (in the limit λ L = 0) is about 1.4 × 10 −47 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 highmass 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 Fig. 5 we project our scan points onto the m H 0 -σ v 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 bb final state. In the two IDM GCE best-fit regions that we have found the dominant DM annihilation channels are either bb or W W * , while for the high-mass region the W W and Z Z channels dominate. Since, however, the current bounds on σ v from dwarf spheroidal galaxies are similar for the bb and W W 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σ 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 Fig. 5 we show the projected limits for 100 h of CTA observation of the Galactic center, taken from [111], assuming annihilation into W W 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 λ L values that characterize 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 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 ± H ± pairs does not depend on λ L , as it only involves gauge couplings and could, at least for sufficiently light A 0 and H ± 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 ± 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 ± − m H 0 ∼ 0.2 GeV regardless of the value of λ 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 λ L values.

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 ∼ m h 0 /2, in which H 0 particles annihilate mostly into bb 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 freezeout 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 W W 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 σ v on the DM velocity is milder, and the GCE can mostly be explained for R ≡ h 2 | IDM / h 2 | DM, Planck ∼ 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 (O(10 −2 ) or less) values of the coupling between H 0 pairs and the standard model Higgs boson and masses m H 0 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 WIMPnucleon 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.