Probing the Inert Doublet Dark Matter Model with Cherenkov Telescopes

We present a detailed study of the annihilation signals of the inert dark matter doublet model in its high mass regime. Concretely, we study the prospects to observe gamma-ray signals of the model in current and projected Cherenkov telescopes taking into account the Sommerfeld effect and including the contribution to the spectrum from gamma-ray lines as well as from internal bremsstrahlung. We show that present observations of the galactic center by the H.E.S.S. instrument are able to exclude regions of the parameter space that give the correct dark matter relic abundance. In particular, models with the charged and the neutral components of the inert doublet nearly degenerate in mass have strong gamma-ray signals. Furthermore, for dark matter particle masses above 1 TeV, we find that the non-observation of the continuum of photons generated by the hadronization of the annihilation products typically gives stronger constraints on the model parameters than the sharp spectral features associated to annihilation into monochromatic photons and the internal bremsstrahlung process. Lastly, we also analyze the interplay between indirect and direct detection searches for this model, concluding that the prospects for the former are more promising. In particular, we find that the upcoming Cherenkov Telescope Array will be able to probe a significant part of the high mass regime of the model.


I. INTRODUCTION
Multiple observations strongly suggest that the Standard Model of particle physics should be extended by at least one additional particle, electrically neutral and colorless, and longlived on cosmological time-scales, dubbed the dark matter (DM) particle [1][2][3][4]. Among the many models that have been constructed over the last decades containing a DM particle, the Inert Doublet Model (IDM) stands out for its simplicity and for its rich phenomenology.
The IDM [5][6][7] postulates the existence of a new scalar field η, with identical gauge quantum numbers as the Standard Model's Brout-Englert-Higgs scalar doublet (Higgs for short), and the invariance of the vacuum under a Z 2 symmetry, under which η is odd while all the Standard Model particles are even. These two simple assumptions have a number of implications. First, the Z 2 symmetry ensures that the doublet η contains an absolutely stable particle which is a DM candidate. Second, the exotic doublet η does not interact at tree level with any of the Standard Model fermions, hence the name "inert". The DM, nonetheless, interacts with the Standard Model via gauge interactions and via the quartic term in the scalar potential |η| 2 |Φ| 2 , with Φ the Higgs doublet. These terms are of utmost importance in the phenomenology of the IDM, since they allow to generate a population of DM particles in the early Universe via thermal freeze-out and they induce potentially observable signals in direct and indirect DM searches [7][8][9][10][11][12][13][14][15][16][17][18], collider searches [7,[19][20][21][22], electroweak precision tests [7,23] and the Higgs diphoton decay rate [24][25][26][27].
Of particular interest in the IDM is the scenario where the DM is entirely constituted by the lightest component of the doublet and where the observed DM abundance Ωh 2 0.12 [28] is generated by thermal freeze-out of this particle. As is well known, this requirement implies for this model a DM particle mass either smaller than the W boson mass, M W 80 GeV, or larger than ∼ 500 GeV (see e.g. [7,8,16,[29][30][31]). In this paper we will concentrate in the latter mass window, and we will investigate the possibility of detecting gamma-ray signals generated by the annihilations of these DM particles in the Milky Way center.
The DM induced gamma-ray flux is comprised of two main components. One component is the prompt gamma-rays, mainly generated by the hadronization of massive gauge and Higgs bosons in the DM annihilation final states, and one lower energy component, consisting of photons of the interstellar radiation field that have been up-scattered to gamma-ray energies due to collisions with the energetic electrons and positrons produced in the annihilations. Since we are interested in the energy spectrum at the highest energies, we will neglect the latter contribution in what follows. Besides, the annihilation also produces sharp gamma-ray spectral features which, if observed, would strongly hint toward an exotic origin of this signal. So far, three different gamma-ray spectral features have been identified in DM scenarios: gamma-ray lines [33][34][35], internal electromagnetic bremsstrahlung [36][37][38][39][40][41] and gamma-ray boxes [42]. Notably, the three spectral features arise in the IDM: the gamma-ray lines arise from annihilations at the quantum loop level into γγ and γZ, the internal bremsstrahlung signal arises from annihilation into W + W − γ through the t-channel exchange of the charged Z 2 odd scalars of the inert doublet, and gamma-ray boxes arise from the annihilation into a pair of Higgs bosons and their subsequent decay in flight into γγ.
Due to the small branching fraction of the process h → γγ, the gamma-ray box produced by the decay in flight of the Higgs boson is fainter than the other two spectral features and will not be considered here.
Since the DM candidate in the IDM possesses a SU (2) L charge, weak gauge bosons could be exchanged between the non-relativistic particles in the initial state of the annihilation process. As argued in [43][44][45], for heavy DM particles the long-range interaction associated to the exchange of the weak gauge bosons can significantly distort the wave function of the initial state particles, therefore the correct description of the annihilation process must include non-perturbative effects, which generically lead to an enhancement of the annihilation cross section. This phenomenon, commonly known as Sommerfeld enhancement, can boost the annihilation signal by many orders of magnitude and has been proved to be pivotal in ruling out some well motivated DM scenarios, such as the Wino DM [46][47][48][49][50][51][52] or the 5-plet minimal DM [53,54] (assuming the DM density follows an Einasto profile in our Galaxy).
The main goal of this paper is to investigate the prospects to observe signals in gamma-rays from the IDM in the high mass regime, including the Sommerfeld enhancement and including not only the channels generating a continuum of gamma-rays, but also those generating sharp spectral features in the energy spectrum.
The paper is organized as follows. In Section II we present a brief overview of the Inert Doublet Model. In Section III we discuss the process of annihilation in the non-relativistic limit and we describe our non-perturbative approach to calculate the cross section. In Section IV we calculate the expected gamma-ray flux and we confront the predictions of the model to limits on continuum gamma-ray fluxes and on sharp gamma-ray spectral features.
In Sections V and VI, we discuss the complementarity between direct detection experiments and gamma-ray instruments in probing the parameter space of the IDM and, lastly, in Section VII we present our conclusions. We also include three appendices discussing various theoretical and experimental constraints on the IDM, technical details of the Sommerfeld enhancement and an estimation of its effect on relic density calculations in the early Universe.

II. DARK MATTER AS AN INERT SCALAR
The IDM is an extension of the Standard Model by one complex scalar field η, which is a singlet under SU (3) C , doublet under SU (2) L and has hypercharge 1/2. Furthermore, the model postulates a discrete Z 2 symmetry, preserved also in the electroweak vacuum, under which the Standard Model particles are even while the extra scalar doublet η is odd.
With this particle content, the Lagrangian can be cast as L = L SM + L η , where L SM is the Standard Model Lagrangian including a potential for the Higgs doublet Φ and L η is the most general Z 2 invariant Lagrangian involving the scalar doublet η where D µ denotes the covariant derivative.
Due to the postulate that the Z 2 symmetry remains unbroken, only the Higgs doublet acquires an expectation value, therefore the doublets can be cast as where v h ≡ −m 2 1 /λ 1 ≈ 246 GeV, G 0 and G + provide the longitudinal components of the of the Z and W + bosons through the Brout-Englert-Higgs mechanism and h is the Standard Model Higgs. On the other hand, the inert sector consists of two charged states H ± , one CP-even neutral state H 0 and one CP-odd neutral state A 0 . Furthermore, the preserved Z 2 symmetry ensures that the lightest particle in the inert doublet is absolutely stable and, if it is neutral, it constitutes a DM candidate; we will assume in what follows that this is the case for the CP-even neutral state H 0 .
The seven parameters in the scalar potential of the model can be recast in terms of the Higgs boson mass M h ≈ 125 GeV, the vacuum expectation value of the Higgs field v h and the DM mass M H 0 , together with the quartic couplings λ 2 , λ 3 , λ 4 and λ 5 . The masses of the remaining inert scalars are given in terms of these parameters by These masses receive corrections at the quantum level. In particular, gauge interactions induce a splitting between the neutral and charged scalar masses which is approximately 356 MeV [29]. However, this contribution can be compensated by an appropriate renormalization of the quartic couplings, resulting in a mass difference which can, in principle, be arbitrarily small. Therefore, in this paper we will take the mass differences among the inert scalars as free parameters, only constrained by, e.g., the perturbative condition |λ i | < 4π in Eq. (4).
The parameters of the IDM are further constrained by the stability of the vacuum [55,56] and by the unitarity of the S-matrix [57,58]. We summarize the constrains we impose in Appendix A, along with various experimental bounds coming from electroweak precision observables and collider searches.
The inert scalar H 0 has the characteristics of a Weakly Interacting Massive Particle (WIMP) because its gauge interactions with the electroweak gauge bosons and its quartic coupling to the Higgs particle are of the required size to thermally produce H 0 of the right amount, via the freeze-out mechanism, to match the observed DM content of our Universe Ωh 2 0.12 [28].
There are two allowed DM mass regimes for H 0 (see e.g., [7,8,16,30,31] there is no longer room to avoid this conclusion by invoking destructive interference effects; see [16] and then, e.g., [59]). Nonetheless, for M H 0 535 GeV, and vanishing quartic couplings, the annihilation rate into gauge bosons is sufficiently small to reproduce the observed DM density. For masses above 535 GeV, the correct relic density can also be obtained if the quartic couplings are appropriately chosen, because their effect is to increase the annihilation cross section. This forms the second, so called, high DM mass regime of the IDM.
The requirement of correct H 0 abundance thus implies larger and larger couplings as the DM mass increases. In fact, an upper limit on the DM mass can be derived by imposing perturbativity on the couplings. With the bounds of Appendix A, and from the relic abundance calculation in Appendix C, we find an upper limit of M H 0 20 TeV. In this paper we will investigate this high mass regime of the IDM and in particular the possible signals in gamma-ray signals from annihilation of H 0 particles in the galactic center. To this end, 1 Three body annihilations of the type H 0 H 0 → W W * → W ff are also important in some regions of the parameter space [15].
we discuss the corresponding annihilation cross sections in the following section.

III. ANNIHILATION CROSS SECTION INTO GAMMA-RAYS
Various annihilation channels contribute to the gamma-ray flux in the IDM. The processes with the largest cross section are the tree-level two-body annihilations into W + W − , ZZ and hh, which generate a gamma-ray flux with a featureless energy spectrum. Processes arising at higher order in perturbation theory have a smaller cross section, however they can contribute significantly to the gamma-ray flux at energies close to the kinematical endpoint of the annihilation and produce a sharp spectral feature in the energy spectrum. This is the case of the tree-level three-body annihilation into W + W − γ, as well as the one-loop annihilations into γγ and γZ. The former was studied in [18] and leads to a bump close to the kinematical end-point of the gamma-ray energy spectrum. This process is sizable especially when H 0 and H ± are relatively close in mass, which is typical in the high mass regime (see Eq. (4)). The latter, on the other hand, were studied in [10] in the low mass regime.
The perturbative approach pursued in that paper, however, cannot be applied to the high mass regime since for very large DM masses the predicted annihilation rates in the galactic center into γ γ and γ Z exceed the upper bound set by unitarity. In fact, for non-realtivistic velocities, the one-loop annihilation cross sections into photons are not suppressed by the DM mass but rather by the W boson mass (see Ref. [60] for a detailed discussion). Similar shortcomings of the perturbative calculation have been pointed out for neutralino DM in the MSSM [61,62], which were solved in Ref. [45] by pursuing a non-perturbative approach.
To calculate the annihilation cross section into the various final states in the high mass regime of the IDM, we thus follow closely the formalism introduced in [43,45,63]. There they use the framework of non-relativistic field theory, which is well motivated by the fact that DM particles move slowly in our Galaxy. 2 The non-relativistic action is obtained by taking the non-relativistic limit description of the components of the inert doublet, which are assumed to be quasi-degenerate in mass, and by integrating out the light particles, namely the Higgs boson and the gauge bosons. In this formalism, it is convenient to introduce 2 A detailed description of the non-perturbative calculation of the annihilaton rate in the IDM can also be found in Ref. [60]. auxiliary fields for the two-body states where s i describes the wave functions of any of the pairs i = (H 0 , H 0 ), (A 0 , A 0 ) and (H − , H + ). Here x is the position of the center of mass of the system and r is the relative position vector for the pair of particles. In terms of these auxiliary fields, the two-body state effective action reads The matrix V (r) represents a central potential consisting of three terms: one specifying the mass splittings among the pairs of particles and the other two describing Yukawa potentials induced by the non-relativistic exchange of gauge bosons (as shown in Fig. 1) and light scalars (as shown in Fig. 2). Thus with Here Moreover, the action contains an absorptive -or imaginary -term, which takes into account that the two body states can annihilate into a Higgs pair or into two gauge bosons.
This term is proportional to the matrix Γ = f Γ (f ) , where f is any final state in which the pairs in the auxiliary field s( x, r) can annihilate into. More concretely, where P i and P f are the total 4-momenta of the initial and final states, the N i are symmetry factors for the initial state particles with where we introduce for convenience the following matrix For the internal bremsstrahlung process, corresponding to the final states W W γ, we directly use Eq. (11), before integrating on the photon energy (see the Appendix of Ref. [18] for details).
These matrices are of interest here because they allow to calculate the annihilation s-wave cross section in the final state f by means of the formula where the matrix d are the Sommerfeld enhancement factors that can be calculated by solving the Schrödinger equation associated to the potential (7), as described in detail in For DM annihilation, Eq. (17) can be cast as The quantities d 11 , d 12 and d 13 are therefore interpreted as non-perturbative enhancement factors that account for the long range interactions between the annihilating DM particles due to the exchange of gauge and Higgs bosons in the non-relativistic limit. As an example, we show in Fig. 3 the absolute value of d 11 , d 12 and d 13 as a function of the DM mass, for the case λ 3 = λ 5 = 0 and λ 4 chosen so that the mass splitting between the charged and the neutral component is 1 GeV (left panel) and 10 GeV (right panel) 3 . We find that, for masses below approximately 2 TeV, the inclusion of these factors in the calculation is irrelevant; however, once the DM mass increases, the enhancement factors dramatically affect the annihilation cross sections in Eq. (18). Furthermore, we find a resonance, which moves to higher DM masses as the mass splitting between the charged and the neutral component is increased. This is in agreement with what was found in Ref. [45] for neutralino DM.
In order to study the impact of the Sommerfeld enhancement in the IDM, we performed a scan over the five dimensional parameter space of the DM sector following the procedure described below. Considering, as observed in Fig. 3, that scalar mass splittings are essential quantities, instead of taking λ 4 and λ 5 as parameters of the scan, we let the relative mass vary logarithmically in between 10 −5 and 1. Then using Eq. (4), we solve for λ 4 and λ 5 and we discard points whose magnitude is greater than 4π. In contrast, the quartic couplings λ 2 and λ 3 , which do not lead to any mass splitting, are randomly sampled on a linear scale. Finally, we take M H 0 in-between 0.5 and 20 TeV and impose the remaining constraints of Appendix A. We then calculate the relic abundance of H 0 for each model by means of micrOMEGAs 3.1 [64] and require that it is in agreement with the observed value Ω DM h 2 = 0.1199 ± 0.0027 within a 30% range.
As discussed in Appendix C, 30% is the error in the relic density calculation that can be expected from our approximation of not accounting for the Sommerfeld effect in the early Universe. The enhancement of the cross section in the channels producing gamma-rays in the final state has important implications for the indirect searches of the inert doublet DM at gammaray telescopes, as we discuss in the next section.

IV. GAMMA-RAY SIGNALS OF THE IDM
The DM induced gamma-ray signal from a given sky region is The astrophysicalJ-factor is here the line-of-sight integral over the squared DM density ρ H 0 , averaged over the observed solid angle ∆Ω, The d(σv) γ /dE γ is comprised of three different components: first, the fairly featureless spectrum generated in the decay and fragmentation of the gauge and Higgs bosons, second, the gamma-ray lines produced by the monochromatic photons emitted in the processes H 0 H 0 → γγ and H 0 H 0 → γZ, and third, the virtual internal bremsstrahlung signal from Concretely, the inclusive differential cross section into photon is where dN f γ /dE γ is the photon multiplicity associated to the two body final states f . When f is a electroweak or Higgs boson pair, we use the parametrization dN γ /dE γ = dN frag γ /dE γ = 0.73 The relative strength of each of these components in the total gamma-ray flux strongly depends on the concrete choice of the parameters of the model. In order to assess the prospects to observe annihilation signals of the IDM, we have calculated the predicted gamma-ray flux for all the viable models of our parameter scan from section III. We include the gamma-ray contributions from the final states W + W − , ZZ, hh, γγ, γZ and W + W − γ and take into account the Sommefeld enhancement for each of these channels 4 , as described in Section III.
In order to better illustrate results, we also selected six benchmark model points (BMPs) displaying qualitatively different spectra. The parameters corresponding to each of these points as well as their predicted gamma-ray energy spectra are shown in Table I. Benchmark points BMP1 and BMP2 produce a very intense gamma-ray line, BMP3 and BMP4 produce a significant virtual bremsstrahlung signal, while BMP5 and BMP6 produce an intense continuum. In the plots, the contributions of the virtual internal bremsstrahlung (VIB), the continuum part, and the γγ and γZ monochromatic lines are shown, respectively, in blue, green, magenta and pink. Considering that the total gamma-ray flux observed by H.E.S.S. telescope falls roughly as E −2.7 [66], our spectra have been multiplied by E 2.7 in order to better appreciate their features at the highest energies.
For all viable points from our scan, we find fairly large annihilation cross sections for both the channels producing continuum gamma-ray emission and those producing sharp gamma-ray spectral features. These signals could therefore be in reach by present and upcoming gamma-ray telescopes. In Fig. 5   Einasto DM density profile using a local DM density of ρ = 0.39 GeV/cm 3 and our various data sets are from the inner Galactic center sky region specified in [66,67] unless otherwise stated. In the same plot we also include limits from the analysis of Ref. [68] using preliminary measurements of the cosmic antiproton-to-proton fraction by the AMS-02 collaboration [69] (solid red line). The two plots in the lower panel instead includes the limits derived by the H.E.S.S. collaboration on the channels γγ and W + W − γ [66], respectively. Furthermore, the Fermi-LAT Collaboration searches for gamma-ray signals from dwarf spheroidal galaxies provide relevant limits [70]. However, the DM annihilation cross sections predictions can be somewhat different in these galaxies, because the DM velocity dispersion is lower there than in the Galactic center region, and consequently we do not include them in our analysis.
To examine the expected reach of the upcoming CTA telescope, we show in the upper plot the projected limits on annihilations into W + W − as derived in [71], assuming 100 hours of observation of a Milky Way center region, and in the lower right plot the limit predictions on monochromatic photons from [72] (after a proper rescaling of their limits on narrow boxed shaped spectra), assuming an observation time of 112 hours of the Galactic center region given in [66,67].
These estimates indicate that present instruments are already sensitive to large regions of the viable parameter space of the IDM and that CTA has good prospects to observe a signals from this model by the observation of a broader continuum excess in the gamma-rays spectrum. Furthermore, the continuum signal flux might be complemented by a simultaneous univocal DM signal in the form of a sharp gamma-ray spectral feature.

V. COMPLEMENTARITY AMONG GAMMA-RAY SIGNALS IN THE IDM
In order to more carefully asses the prospects to observe signals of the IDM, we derive dedicated limits on each of the IDMs from our parameter scan. For each model's cross section limit, we then define a maximal boost factor (BF) that corresponds to how much the model's predicted gamma-ray signal can be increased before it saturates its derived limit.
To derive signal limits on each IDM induced continuum signal, we use the data collected by the H.E.S.S. instrument and closely follow the method pursued in [67], which compares the gamma-ray fluxes measured in a "search region" and in a "background region" around the in the colors (to be used also for future references) red, green, cyan and blue, respectively.
Among the viable models, the six benchmark points of Table I are highlighted in the plot with their corresponding tag. Notably, there are many points, especially with mass above ∼ 2 TeV which are already excluded by observations with the H.E.S.S. instrument. Furthermore, for most of the points the BF value is constrained to be smaller than ∼ 10. Therefore, an improvement in sensitivity of gamma-ray telescopes to this type of exotic continuum flux by a factor of ∼ 10, which seems to be feasible with the upcoming CTA (see e.g. Ref. [71] and our Fig. 5), could suffice to cover all the signal predictions from the IDM, assuming that the DM halo distribution follows the Einasto profile. For a Navarro-Frenk-White profile of the DM distribution, theJ factor in the target region is a factor of two smaller [67], hence the annihilation signal would in this case be a factor of two fainter and the prospects for detection, somewhat poorer.
We have also calculated the maximal boost factor BF from the non-observation of the sharp gamma-ray spectral features produced by the final states γγ, γZ and W + W − γ. These limits on the IDM were derived following the procedure pursued by the H.E.S.S. collaboration in [66], which adopts a phenomenological background model defined by seven parameters.
The result of this procedure is illustrated in Table II Fig. 6. We find again models which are already ruled out by present observations, especially at large DM masses. 5 5 The strengthening of the limits at M H 0 ∼ 600 GeV is due to a dip around E γ ∼ 700 GeV in the gamma- Furthermore, with an increase in sensitivity by a factor ∼ 10, which is likely to be achieved with the upcoming CTA (see e.g. [53,72]), a significant part of the IDM parameter space will be probed, thus opening the exciting possibility of observing unambiguous signals from DM annihilation at future gamma-ray telescopes. Unfortunately, to guarantee the observation of a sharp feature in the gamma-ray spectrum a larger increase in sensitivity is necessary, concretely by a factor ∼ 100, assuming the Einasto profile.
From the above discussions it apparently becomes relevant to investigate the potential complementarity between the searches for a continuum exotic flux and a sharp spectral feature. In Fig. 7 we illustrate this complementary. For each model in our scan, the required ray flux measured by the H.E.S.S. collaboration, and which is possibly due to a downward statistical fluctuation.
boost factor for a model to become excluded by the conti-nuum spectrum constraints (continuum BF) is confronted to the required boost factor to become excluded by the sharp spectral feature (feature BF). It follows from the figure that, for most of the points, the former provides a stronger constraint than the latter, with the difference between the BFs being more notable as the DM mass increases. This behaviour is a consequence of the requirement on the parameters to correctly reproduce the DM density via thermal freeze-out. As argued in Section II, larger and larger quartic couplings are required when the DM mass increases.
As a result, the annihilation rates into W + W − , ZZ and hh, which can be induced by quartic coupling interactions (as follows from Eqs. (13), (14) and (15)), are enhanced compared to γγ and γZ, which are induced only by gauge interactions (as follows from Eq. (12)).
We would like to briefly comment on the difference between the IDM and two other minimal DM scenarios with a DM candidate in a 5-plet and 7-plet representations of SU (2) L [53,54,73]. In the these papers it was concluded that if the DM candidates account for all the DM then they are excluded if their masses are below 20 TeV. Although we cover similar masses in the IDM, not all them are excluded, as it is shown in Fig. 7. The underlying reason for this is related to the mass splittings between the charged and the neutral components of the inert states. In the minimal 5-plet and 7-plet DM scenarios, the mass splitting is set by radiative corrections and is fixed to a constant value. In the IDM there is more freedom, and the mass splitting can be larger than the quantum effect. In fact, as already discussed, for large DM masses, large quartic couplings are typically needed to achieve the right relic abundance. Unless a cancellation between λ 4 and λ 5 takes place, the mass difference between H + and H 0 becomes relatively larger. This leads to relatively smaller Sommerfeld effects in comparison to these minimal DM models (even if it is still large for the heaviest DM masses in IDM). Another important reason is that larger SU (2) L multiplets will contain particles with larger electric charges. This leads to larger annihilation cross sections for all gauge mediated annihilation channels and, in particular, increases monochromatic gamma-ray signals.

VI. COMPLEMENTARITY WITH DIRECT DETECTION
A complementary avenue to probe the high mass regime of the IDM is direct detection.
The spin-independent scattering cross section of DM particles with nuclei receives in this model two different contributions. The first one is induced by the t-channel exchange of a Higgs boson with the nucleon n [7]: where f ≈ 0.3 is a form factor with its precise value taken from micrOMEGAs 3.1 [64] and M n = 0.939 GeV is the nucleon mass. This cross section is suppressed for large DM masses and for small DM-Higgs coupling, which corresponds to |λ 3 + λ 4 + λ 5 |. The second contribution is induced by one-loop exchange of gauge bosons [17,29,31], which is independent of the quartic couplings and which sets a lower limit on the interaction cross section of σ SI 2.6 × 10 −46 cm 2 independently on the DM mass. We show in Fig. 8 the predicted spin-independent scattering cross section with protons for our sample of viable points and compare them to the limit from the LUX experiment as well as to the projected reach of XENON1T [74] and LZ [76]. As apparent from the plot, the current data from the LUX experiment barely constrain the viable parameter space of the IDM. On the other hand, and due to the above-mentioned lower limit on the interaction cross section induced by the oneloop exchange of gauge bosons, the upcoming XENON1T (LZ) experiment should observe a signal of the IDM for thermally produced DM particles with M H 0 1.6 TeV (13 TeV). We note that many of the points of our scan have an interaction cross section σ SI 5×10 −45 cm 2 , which is well within the reach of XENON1T.
We finally confront the constraints from direct detection DM searches to those from indirect detection searches in gamma-rays. In Fig. 9, the left plot (right plot) shows the viable DM models in the plane of the ratio σ LUX /σ against the required boost factor of the broad DM gamma-ray signal (sharp spectral feature) to reach the current H.E.S.S.
constraints. Here σ is the predicted spin-independent cross section with a proton and σ LUX is the upper limit from LUX for that model. Points with σ LUX /σ < 1 are then ruled out by the LUX experiment and points with BF< 1 are ruled out by the H.E.S.S. instrument.
The values of BF and σ LUX /σ for our six benchmark points are indicated in the figure and displayed for reference in Table III

VII. CONCLUSIONS
We have investigated the gamma-ray spectrum produced in DM annihilations in the center of the galaxy for the high mass regime of the Inert Doublet Model (IDM). We have found that, in order to satisfy the requirements of unitarity on the annihilation cross sections, it is necessary to account for the so-called Sommerfeld enhancement. This is a non-perturbative effect arising from the exchange of gauge and Higgs bosons between non-relativistic annihilating DM particles. In the mass regime under consideration, such exchange induces long range interactions that lead to a significant modification of the annihilation cross sections.
We have argued that including that effect is crucial for phenomenological studies of indirect DM signals from the galactic center, specially for masses much larger than 1 TeV and small mass splittings between the charged and the neutral scalars. We have also showed that such effect is much less important for the DM production in the early Universe (see Appendix C).
We have calculated the impact of the Sommerfeld enhancement in the framework of the effective field theory resulting from the non-relativistic limit of the inert scalar particles. The main ingredients to consider are the potential matrices of Eq. (7), encoding the long-range effects, and the matrices of Eqs. (12), (14), (13) and (15), which describe the annihilation processes. In appendix B we have described succinctly how to use these matrices in order to the obtain the annihilation cross sections. Using this formalism, we have been able for the first time to reliably calculate the gamma-ray spectrum. It receives contributions from (i) a featureless soft part arising from annihilations into W + W − , ZZ and hh pairs, Subsequently, we have analyzed the interplay between this indirect search and the direct DM searches with the LUX experiment. The corresponding results are shown in Fig. 9, and for our benchmark points, on table III. Current direct DM searches are not sensitive enough to detect signals of the IDM. Nevertheless, the upcoming XENON1T experiment and the projected LZ experiment will be able to close in on the viable parameter space of the model.
Finally, we would like to comment on the sensitivity of the upcoming Cerenkov Telescope Array to the IDM. As shown in Fig. 5, a significant part of the viable models of our scan can be potentially probed by CTA since most of them produce a continuum spectrum that is within a factor ten from current experimental sensitivity.

Note Added
During the last stages of this work, we learned of the analysis of Ref. [77], where the gamma-ray signals of the IDM in its high mass regime are also studied. In that paper, however, the Sommerfeld effect was neglected and the contribution to the spectrum from gamma-ray lines and from virtual internal bremsstrahlung was not considered.
• Unitarity of the S-matrix on scalar to scalar, gauge boson to gauge boson and scalar to gauge boson scatterings at the perturbative level furthermore requires that [57,58] • Electroweak precision observables constraints the contributions to the Peskin-Takeuchi S, T, U parameters to remain in the region ∆S = 0.06 ± 2 × 0.09, ∆T = 0.1 ± 2 × 0.07 with a correlation coefficient of +0.091 (when ∆U is fixed to zero, which is appropriate for the IDM). The contribution from the IDM can be calculated as in, e.g., [24]. This typically prohibit large mass splittings among inert states, but for DM masses with M H 0 500 GeV relatively small splittings are already required, especially when combined with the relic density constraint [31].
Besides these theoretical constraints on the parameters of the IDM, there are also limits from experimental searches of the inert scalars.
• LEP bound comes first from that the decay channels Z → A 0 H 0 , Z → H + H − , W ± → A 0 H ± and W ± → H 0 H ± would alter the gauge bosons measured mass widths.
As a good approximations, these implies that  [24][25][26][27]. These bounds are typically only relevant for masses below M h /2, and thus play a little role for inert scalar particle masses well above. Direct di-lepton searches have also been shown to restrict the inert scalar masses in the region of M H 0 60 GeV and M A 150 GeV [80].
From these constraints it is clear that the IDM is strongly restricted if the new states are at sub 100 GeV masses and not so constrained for masses above 500 GeV. In addition there are various astrophysical constraints (see, e.g., the review in [56]), and direct detection constraints. Due to the relevance of the latter for this work, they are discussed in the main text.
Appendix B: Algorithm for the Sommerfeld Enhancement In this appendix we briefly describe the algorithm to calculate the Sommerfeld enhancement factors, given the potential of Eq. (7) and the annihilation matrices of Eqs. (12), (14), (13) and (15).
As shown in [45], the Sommerfeld effect is encoded, in the basis of the inert pair states (H 0 , H 0 ), (A 0 , A 0 ) and (H − , H + ), in a 3×3 matrix g(r), which satisfies the following secondorder differential equation where v is the relative velocity of the initial state particles and V (r) is given in Eq. (7). The boundary conditions to solve the differential equation can be determined by analyzing the behavior of the solution g(r) at r = 0 ant r → ∞. At the origin, On the other hand, for large values of r, the matrix g(r) depends on the mass splitting between the pair states. If the mass splitting δm ij associated to the inert pairs i and j is smaller than the initial kinetic energy 1 4 M H 0 v 2 , then there is enough energy to produce on-shell states of the corresponding pair and, therefore, the matrix element g ij (r) at infinity behaves as an out-going wave with momentum given, according to Eq. (B1), by The corresponding boundary condition is In the opposite case, namely when δm ii > 1 4 M H 0 v 2 , there is not enough energy to produce on-shell states of the corresponding pair, and therefore the matrix elements g ij (r) decay exponentially at infinity. Hence The boundary conditions at r = 0 and r → ∞ then allow to find g(r) by solving Eq. (B1).
Once the solution is obtained, the oscillating phases of g(r) at large values of r can be factorized by casting the matrix as where d is a 3 × 3 matrix which contains the non-perturbative enhancement factors due to the exchange of scalar and gauge bosons in the initial state. 6 Finally, the s-wave cross section for the annihilation of the pair i into a final state f can be determined from where N i = 1/ √ 2 for initial states with identical particles, and N i = 1 otherwise. Notice that when the potential in Eq. (B1) is negligible then d = 1 1, and therefore Eqs. (18) and (B7) reduce to the standard expressions for calculating the cross section in the s-wave limit.
Let us first consider pairs of the co-annihilating species H 0 , A 0 , H + , H − . The subspace generated by such pairs can be decomposed into one self-conjugate isospin singlet |m I = 0, I = 0 Y =0 , one self-conjugate triplet |m I , I = 1 Y =0 and one isospin triplet |m I , I = 1 Y =1 and its corresponding complex conjugate |m I , I = 1 Y =−1 . In terms of these states, the co-annihilating pairs with charge Q = 0, Q = 1 and Q = 2 read Q = 2 : H + H + = |m I = 1, I = 1 Y =1 (C1) Q = 1 : Because the isospin is a good quantum number, the non-relativistic potential is proportional to the identity in subspaces with definite isospin and hypercharge. Concretely, This equation and the previous transformation matrices allow to express the potential in the basis of co-annihilating pairs. Comparing the resulting potential for Q = 0 with Eq. (7), it is possible to solve for the potential V IY , after neglecting the vev triggering the electroweak symmetry breaking (and therefore the mass splittings, the gauge boson masses and the scalar potential, which are all proportional to v h ). The result reads: An analagous procedure can be applied to the annihilation matrices. That is Γ |m I , I Y = Γ IY |m I , I Y (C7) Once again by employing this and Eqs. (C1), (C2) and (C3), we write the annihilation matrices in the basis of co-annihilating pairs. Furthermore by comparing the result with the addition of Eqs. (12), (14), (13) and (15), we solve for the annihilation matrix In each subspace with definite hypercharge and isospin, the potential and the annihilation matrices take, neglecting the gauge boson masses, the simple form of a Coulomb potential, thus allowing to analytically estimate the Sommerfeld enhancement in each subspace. Concretely, the total annihilation cross sections reads where the function z/(e z − 1), with z = πα IY /v, is the well known Sommerfeld enhancement factor associated to a Coulomb potential [63,82]. Besides, the factor of 2 in the numerator comes from the fact that the DM in the IDM is its own antiparticle, the symmetry factor 4 2 is the total number of pairs that can be constructed from the co-annihilating species, and the factor (2I + 1) is the multiplicity of each isospin state Lastly, taking the thermal average of this expression and using the instantaneous freezeout approximation, one can estimate the DM relic density as [31] Ωh 2 = 1.07 × 10 9 GeV −1 where the inverse freeze-out temperature x f = M H 0 /T f can be found by solving Remarkably, Eq. (C12) can be used to calculate the relic density even where there is no Sommerfeld enhancement, by taking α IY → 0. In this limit, Eq. (C11) reduces to (σv) ef f 2.4 + 5.2 [(2λ 3 + λ 4 ) 2 + 3λ 2 4 + 6λ 2 5 ] (M H 0 /530 GeV) 2 × 10 −26 cm 3 /s.
From this, we find a deviation of at most 10% from the result of the scan of Section III, which was derived using micrOMEGAs 3.1 [64]. We also find that the difference becomes stronger when the mass splitting is very large, as expected from the fact that in this regime the approximation of taking the SU (2) L symmetric limit is worst. When the Sommerfeld enhancement is taken into account, the disagreement with respect to perturbative result is at most of 30%, in remarkable agreement with what has been found for Higgsino DM [51,83], another SU (2) L doublet candidate. We thus conclude the the Sommerfeld effect in the early Universe is less dramatic than in the galactic center.
For a given set of quartic couplings, using Eq. (C14) and Ω DM