Anatomy of the Inert Two Higgs Doublet Model in the light of the LHC and non-LHC Dark Matter Searches

The inert Two Higgs Doublet Model (i2HDM) is a theoretically well-motivated example of a minimal consistent Dark Matter(DM) model which provides mono-jet, mono-Z, mono-Higgs and Vector-Boson-Fusion+Missing Transverse Momentum signatures at the LHC, complemented by signals in direct and indirect DM search experiments. In this paper we have performed a detailed analysis of the constraints in the full 5D parameter space of the i2HDM, coming from perturbativity, unitarity, electroweak precision data, Higgs data from LHC, DM relic density, direct/indirect DM detection and LHC mono-jet analysis, as well as implications of experimental LHC studies on disappearing charged tracks relevant to high DM mass region. We demonstrate the complementarity of the above constraints and present projections for future LHC data and direct DM detection experiments to probe further i2HDM parameter space. The model is implemented into the CalcHEP and micrOMEGAs packages, which are publicly available at the HEPMDB database, and is ready for a further exploration in the context of the LHC, relic density and DM direct detection.


Introduction
The evidence for dark matter (DM) is well-established from several independent cosmological observations, including galactic rotation curves, cosmic microwave background fits of the WMAP and PLANCK data, gravitational lensing, large scale structure of the Universe, as well as interacting galaxy clusters such as the Bullet Cluster. Despite these large-scale evidences, the microscopic nature of the DM particles remains unknown, since no experiment so far has been able to claim their detection in the laboratory and probe their properties. Potentially, DM can be produced at the LHC and probed in the DM direct detection (DD) underground experiments. The fundamental importance and vast experimental opportunities make the search for and investigation of DM one of the key goals in astroparticle physics and high energy physics (HEP), worthy of the intense efforts undertaken by the physics community.
At the other end of the length scale, the Standard Model (SM) of particle physics recently demonstrated its vitality once again. The scalar boson with mass m H ≈ 125 GeV found at the LHC [1,2] closely resembles, in all its manifestations, the SM Higgs boson. Since the SM cannot be the ultimate theory, many constructions beyond the SM (BSM) have been put forth, at different levels of sophistication. Yet, without direct experimental confirmation, none of them can be named the true theory beyond the SM.
One way the particle theory community can respond to this situation is to propose simple, fully calculable, renormalizable BSM models with viable DM candidates. We do not know yet which of these models (if any) corresponds to reality, but all models of this kind offer an excellent opportunity to gain insight into the intricate interplay among various astrophysical and collider constraints. We call here these models as Minimal Consistent Dark Matter (MCDM) models. MCDM models which can be viewed as toy models, are self-consistent and can be easily be incorporated into a bigger BSM model. Because of these attractive features, MCDM models can be considered as the next step beyond DM Effective Filed Theory (EFT) (see e.g. [3,4,5,6,7,8,9,10,11,12,13,14,15]) and simplified DM models (see e.g. [16,17,18,19,20,21,22,23]) models.
In this paper, we explore, in the light of the recent collider, astroparticle and DD DM experimental data, the inert Two-Higgs Doublet Model (i2HDM), also known as the Inert doublet model. This model is easily doable with analytic calculations, its parameter space is relatively small and can be strongly constrained by the present and projected future data. The model leads to a variety of collider signatures, and, in spite of many years of investigation, not all of them have yet been fully and properly explored. It is the goal of the present paper to investigate in fine detail the present constraints and the impact of the future LHC and DD DM data on the parameter space of this model.
The i2HDM [24,25,26,27] is a minimalistic extension of the SM with a second scalar doublet φ 2 possessing the same quantum numbers as the SM Higgs doublet φ 1 but with no direct coupling to fermions (the inert doublet). This construction is protected by the discrete Z 2 symmetry under which φ 2 is odd and all the other fields are even. The scalar Lagrangian is with the potential V containing all scalar interactions compatible with the Z 2 symmetry: All free parameters here are real, 1 which precludes the CP -violation in the scalar sector. There is a large part of the parameter space, in which only the first, SM-like doublet, acquires the vacuum expectation value (vev). In the notation φ 0 i = v i / √ 2, this inert minimum corresponds to v 1 = v, v 2 = 0. In the unitary gauge, the doublets are expanded near the minimum as The Z 2 symmetry is still conserved by the vacuum state, which forbids direct coupling of any single inert field to the SM fields and it stabilizes the lightest inert boson against decay. Pairwise interactions of the inert scalars with the gauge-bosons and with the SM-like Higgs H are still possible, which gives rise to various i2HDM signatures at colliders and in the DM detection experiments.
The idea that the symmetry-protected second Higgs doublet naturally produces a scalar dark matter candidate was first mentioned more that 30 years ago [24]. However, the real interest in phenomenological consequences of the i2HDM woke up in mid-2000 and intensified in the last few years. Its simplicity, predictive power, rich yet manageable parameter space, makes it an ideal playground for checking its compatibility with the DM relic density, with the results of the direct and indirect DM searches, and with collider searches of various BSM signals.
Assuming that the lightest inert scalar is the only DM candidate, one typically finds that the low-mass region, below about 50 GeV, is excluded by the relic density constraints coupled with the LHC constraints on the invisible Higgs decay [28,29,30]. The funnel region, with the DM mass close to M H /2, the intermediate, 100-500 GeV, and the high mass regions are still compatible with data and lead to interesting predictions at colliders. Additional theoretical constraints on the parameter space and DM candidate properties can be deduced from assumptions of full stability of the i2HDM up to the Planck scale [31,32] or of multi-doublet Higgs inflation [33]. The i2HDM can also produce signals for direct [34] and indirect DM search experiments via heavy inert scalar annihilation, which can be detectable via γ-rays [35,36,37] or via its neutrino [38,39] and cosmicray signals [40].
The i2HDM can also have interesting cosmological consequences. Being an example of 2HDM, it possesses a rich vacuum structure, which evolves at high temperatures [41,42,43]. This opens up the possibility within i2HDM that the early Universe, while cooling down, went through a sequence of phase transitions including a strong first-order phase transitions [44,45,46,47,48,49,50]. Such analyses are capable of restricting the parameter space; for example, the recent study [50] showed that combining the strong first-order phase transition with other astroparticle and collider constraints gives preference to the funnel region.
There has also been a number of studies on collider signatures of the i2HDM. They focus either on specific processes such as SM-like Higgs decays to γγ and γZ [51,52,28,53], multi-lepton plus missing E T [54,55,56] with as many as five leptons [57], dijet plus missing E T [58] and dileptons accompanied with dijets [56]. Other works present combined analyses of astroparticle and collider constraints [59,60,29,61,57]. Comparing the i2HDM predictions with the electroweak precision data, the measured SM-like Higgs properties, the non-observation of long-lived charged particles and other exotic signals, and finally the astroparticle observations, allows one to significantly restrict the i2HDM parameter space. Ref. [59] takes into account one-loop corrections to the masses and, for a part of the numerical scans, includes the additional theoretical constraint that the perturbativity and stability be satisfied up to large scale Λ. Another recent analysis [29] gave a detailed account of these constraints. For specific benchmark points or benchmark planes in the surviving parameter space, it predicted the cross section of pair production of inert scalars followed by various modes of their decay. As for the specific signatures of the i2HDM at the LHC, dileptons and mono-Z signals were mentioned. An earlier analysis [60] investigated multilepton, multijet, mono-Z, and several channels for the mono-jet with large missing E T . The version of i2HDM equipped with Peccei-Quinn U (1) symmetry spontaneously broken to Z 2 was investigated in [61]. Here, dark matter acquires a second component, the axion, which changes the DM phenomenology. It is also possible to hunt for i2HDM at the future colliders, via searching for new scalars and reconstructing the potential [62] or by accurately measuring the SM-like Higgs couplings and deducing patterns of the deviations from the SM [63].
In the present work, to these many studies on the i2HDM, we add: • detailed combined analysis of the i2HDM model in its full five-dimensional (5D) parameter space, taking into account perturbativity and unitarity, LEP and electroweak precision data, Higgs data from the LHC, DM relic density, direct/indirect DM detection complemented by realistic (beyond-the-parton-level) LHC mono-jet analysis at the LHC; • quantitative exploration of the surviving regions of parameters, including very fine details and qualitatively new region not seen in previous studies, which is enabled by our extensive numerical scans; • a combination of different processes giving the LHC mono-jet signatures: those with direct DM pair production and those with coproduction of DM and another scalar with a close mass from the inert multiplet; • implication of experimental LHC studies on disappearing charged tracks relevant to high ( 500 GeV) DM mass region; • separate, equally detailed analyses for the assumptions of the DM relic density being fitted to the Planck results or under-abundant, allowing thus for additional allowed regions of the parameter space.
All these points above are in close focus of the present paper where we have performed a comprehensive scan and study of the full parameter space of the i2HDM model. In addition we have performed an independent implementation and validation of the model in two gauges including Higgs-gluon-gluon and Higgs-photon-photon effective couplings, and we made it public together with the LanHEP model source. The paper is organised as follows. In Sect. 2 we discuss the i2HDM model parameter space, implementation, theoretical constraints as well as constraints from LEP and electroweak precision data. In Sect. 3 we discuss results of a comprehensive scan of the i2HDM parameter space and combined constraints considering both cases, when the relic density is "just right" and in agreement with the Planck data as well as the case of under-abundant DM relic density. In this section we also present the reach of LHC studies in the high DM mass region using results on disappearing charged tracks. In Sect. 4 we present results on future projections of the LHC and DM DD experiments in combination with all previous constraints. Finally, in Sect. 5 we draw our conclusions.

Constraints from the Higgs potential
In order to represent a viable model, the potential (2) must be bounded from below and must have a neutral, non charge-breaking vacuum. The former requirement leads to the well-known restrictions on the free parameters of the model: The absence of the charge-breaking vacuum is guaranteed if one assumes This is a sufficient but not necessary condition for the vacuum to be neutral. A neutral vacuum can also be achieved for positive λ 4 − |λ 5 | with appropriate m 2 1 and m 2 2 . However in this case the lightest DM candidate will be the charged scalar. Condition (5) avoids this situation as well.
Once these restrictions are applied, the vacuum is neutral, and one can calculate the masses of the physical Higgs bosons. In addition to the SM-like scalar H, one gets charged h ± and neutral h 1 , h 2 scalars. It is well known that the two neutral scalars of the i2HDM have opposite CPparities, but it is impossible to unambigously assign which of them is CP -even and which is CPodd. In the absence of any suitable vertex, the model has two CP -symmetries, h 1 → h 1 , h 2 → −h 2 and h 1 → −h 1 , h 2 → h 2 , which get interchanged upon basis change φ 2 → iφ 2 . Either can be used as "the CP -symmetry" of the model, making the specification of the CP properties of h 1 and h 2 a basis dependent statement. Therefore, we denote the two neutral inert scalar masses as M h 1 < M h 2 , without specifying which is scalar and pseudoscalar. The masses of the physical scalars are The mass differences, written as highlight the role of parameters λ 4 and λ 5 and are consistent with (5). It should also be stressed that the parameters λ 1 and m 2 1 correspond to the Higgs potential in the SM, and can thus be fixed by the values of the VEV and Higgs mass.
One also notices that the sign of λ 5 is phenomenologically irrelevant: flipping the sign of λ 5 would only lead to swapping the CP -parities of the inert neutral scalars, which are unobservable anyway. In order to eliminate double-counting, we make the standard choice of λ 5 < 0, and introduce the shorthand notation λ 345 = λ 3 + λ 4 + λ 5 . The latter parameter plays an important phenomenological role, as it governs the Higgs-DM interaction vertex Hh 1 h 1 . With all these conventions, we describe the five dimensional parameter space of i2HDM with the following phenomenologically relevant variables: Another set of theoretical constraints comes from the symmetry breaking patterns in i2HDM [24] and from the fact that the potential can have two minima at different depths. Following [44], we introduce R = λ 345 /2 √ λ 1 λ 2 , which satisfies R > −1. Requiring that the inert vacuum corresponds to the global minimum leads to the following conditions on the parameters of the potential, apart from m 2 1 > 0: In Fig. 1 we visualise these restrictions on the (m 2 1 , m 2 2 ) plane for the three choices of R. The inert, v 1 = v, v 2 = 0, and pseudoinert, v 1 = 0, v 2 = v, vacua can coexist only when R > 1, which is shown by the dashed region in Fig. 1 (a). For R > 1, the second line in Eq. (9) is a stronger condition than the first line and it guarantees that the inert minimum is the deepest one. This condition is shown in Fig. 1 (a) by the solid black line. 1 , m 2 2 ) plane coming from the requirement that the inert vacuum is the deepest minimum of the potential. The three cases correspond to (a) R > 1, (b) 0 < R < 1, (c) −1 < R < 0. Light and dark grey correspond to models with an inert v 1 = v, v 2 = 0 and a pseudoinert v 1 = 0, v 2 = v vacuum, respectively, while the blue region in between corresponds to the mixed vacuum, when both v 1 and v 2 are non-zero. The dashed region in the left plot shows the region when the inert and pseudoinert minima coexist but have different depths.
Rewriting conditions (9) for the physical parameters we get the constraint on the Higgs potential in the following compact final form: and where λ 1 ≈ 0.129 is fixed in the Standard Model by the Higgs mass (6). The latter condition places an upper bound on λ 345 for a given DM mass M h 1 .

Model implementation
We have implemented the i2HDM into the CalcHEP package [64] with the help of the LanHEP program [65,66] for automatic Feynman rules derivation. The effective Hgg and Hγγ vertices were included and the model was cross-checked in two different gauges to ensure a correct, gauge invariant implementation. It is publicly available at the High Energy Physics Model Data-Base (HEPMDB) [67] at http://hepmdb.soton.ac.uk/hepmdb:0715.0187 together with the LanHEP source of the moedl. The model is implemented in terms of the five independent parameters defined in Eq. (8), consisting of three physical masses and two couplings. We found this choice the most convenient for exploration of i2HDM phenomenology and constraints of its parameter space. We should stress that the M h 1 and M h 2 parameters conveniently define the mass order of the two neutral inert states independently of their CP properties. This choice is especially convenient and relevant for collider phenomenology since, as we discussed above, one can not assign (or determine) the CP parity of each neutral inert scalar.
To explore the phenomenology of the i2HDM we need to consider other constraints on its parameter space in addition to those coming from vacuum stability which we discussed above.

Constraints from perturbativity and unitarity
The first requirement we impose on the quartic couplings in (2) is that their value is such that perturbative calculations can be trusted in the model. The most effective way is to impose perturbative unitarity on all the scattering processes involving the scalars. Following [51], we impose this condition on the full scattering matrix, which leads to the following bounds on combinations of couplings e i : where e 1,2 = λ 3 ± λ 4 , e 3,4 = λ 3 ± λ 5 , e 5,6 = λ 3 + 2λ 4 ± 3λ 5 , e 7, The parameter λ 1 is fixed by SM-Higgs mass and the vacuum expectation value. One can verify that the constraints given by Eq. (12) imply that all quartic couplings in (2) are bound to be smaller than 8π, thus within the perturbative regime. The perturbativity constraints can also be used to find upper bounds on the two input couplings we defined in the previous section, i.e. λ 2 and λ 345 . From e 10 one finds: where λ max 2 is a function of model parameters, while from e 5 = 3λ 345 − (2λ 3 + λ 4 ), combined with e 10 in the limit λ 2 = 0, we obtain an upper bound for λ 345 : where we expanded at leading order in the small coupling λ 1 , and the lower bound comes from the stability of the potential. This limit, derived from the constraints on e 5 and e 10 is not actually the most stringent one: in the limit of λ 2 → 0 we have found that the biggest value for λ 345 is realised in the |λ 4,5 | → 0 limit when λ 3 8π 2 = 4π and respectively λ 345 4π. After expansion in the small coupling λ 1 , the upper limit on λ 345 in the small λ 2 limit reads as while for finite λ 2 the limit can be found numerically.
One should also stress that the vacuum stability condition given by Eq.(11) sets an important constraint on the maximum value of λ 345 in the small M h 1 region (which is the region of our special interest because of collider phenomenology constraints as we discuss below). This can be seen from Eq.(11) which can be written as: In Fig.2, we present viable parameter space in the (λ 345 , λ 2 ) plane after constraints from Eq. (12) as well constraints from scalar potential, given by Eq.s (10)(11)16). To produce this plot  we have performed the wide random scan to cover the full five-dimensional parameter space of the model, with the following chosen range for the model parameters: 10 GeV < M h1,h2,h + < 1000 GeV The colour map in Fig

Constraints from LEP and electroweak precision data
Very strong constraints on the i2HDM arise from precision data and searches from LEP experiments. First of all, the model should respect the precise measurements of the W and Z widths which lead to the following lower limit on the odd scalar masses: to make sure that Γ(W + → h 1 h + , h 2 h + ) and Γ(Z → h 1 h 2 , h + h − ) decay channels are kinematically forbidden.
While studying the phenomenology of the i2HDM, we should also make sure that Electroweak Precision Test (EWPT) data is respected. As we know, EWPT can be expressed in terms of three measurable quantities, called S, T, and U, that parameterise contributions from beyond standard model physics to electroweak radiative corrections [68]. The contribution to the S and T parameters [26] can be written as where where the function f c (x, y) is defined by We have written the contributions to S and T in a form which demonstrates explicitly their symmetry with respect to swapping h 1 ↔ h 2 , pointing again to the fact that one can not distinguish their CP properties. With U fixed to be zero, the central values of S and T , assuming a SM Higgs boson mass of m h = 125 GeV, are given by [69]: with correlation coefficient +0.91. The effect of the constraints on S and T is presented in Fig.3, where panels a) and b) present the colour map of the S and T parameters respectively in the (M h + , M h 2 ) plane. One can see that the T variable is more sensitive than S to this mass split, thus only modest splits are allowed by EWPT data. Finally, Fig.3 c) presents the colour map of the (M h + − M h 2 ) split in the (S, T ) plane together with the 65% and 95% exclusion contours, based on a χ 2 with two degrees of freedom. One can see that EWPT data prefer a modest positive (M h + − M h 2 ) mass split below about 100 GeV, which is mainly defined by the T parameter, while the role and the respective range of variation of S is milder. One should stress that it is crucial to take into account the correlation between S and T and combine limits from these two parameters. This combination gives a much stronger limit on the parameter space, in particular We also excluded the region defined by the intersection of the conditions below: which is excluded by LEP data since it would lead to a visible di-jet or di-lepton signal as demon-strated in [70] where a reinterpretation in the i2HDM of a LEP-II limit of the second neutralino production in the MSSM is presented. A more detailed analysis of this specific region of the parameter space -low M h 1 and M h 2 with large enough mass gap -was studied recently [71]. One should also mention that e + e − → h + h − production at LEP2 sets as found in [72] as a result of the re-interpretation of LEP-II limits on charginos.

Constraints from LHC Higgs data
The LHC Higgs data also poses an important constraint on the i2HDM parameters space in the form of constraints on the couplings of the SM-like Higgs boson. A collection of combined fits from the Run-I data, for both ATLAS and CMS, can be found in [73]. In the i2HDM, the leading effects is encoded in two observables: the decays of the Higgs into two Dark Matter scalars, which is kinematically open when m h 1 < M H /2; and the contribution of the charged Higgs loops to the H → γγ rates. In principle, therefore, we would need to do a two-parameter fit of the available Higgs data. None of the fits presented in [73] can therefore be directly applied in our case. A simpler possibility is, instead, to consider the best possible bound from the available fits on the two parameters, in turn. We follow this simpler procedure, confident that it will lead to a somewhat more conservative estimation of the bounds. For the invisible Higgs branching ratio, we consider the bound coming from the dedicated ATLAS search [74], requiring at the 95% CL, which is comparable with a 36% limit from the combined CMS analysis [75] 2 . For the second observable, i.e. the di-photon rate, we consider the result from the combined fit on the signal strength in the di-photon channel [73]: where the central value variation is given at the 95% CL 3 . A sufficiently light charged Higgs with sufficiently large λ 3 coupling to the SM Higgs boson which deviate radiative H → γγ decay beyond the quoted limit, is excluded. It should be noted that we would expect a proper 2-parameter fit to lead to stronger constraints that the ones we use, however the qualitative impact of the constraints should be unchanged. For example, the partial decay width of the Higgs into DM which is defined by where g W , is the weak coupling constant, provides the following bound on λ 345 : where Br max invis = 0.28 is the current bound on the maximal value of branching ratio of the Higgs boson decay into invisible mode. The above limit on λ 345 is M h 1 dependent: for M h 1 /M H 1 it is about 0.019, while for M h 1 closer to the threshold, e.g. 60 GeV, the limit on λ 345 increases almost by a factor of two and reaches a value of 0.036. In addition we have included the limit from H → h 2 h 2 when h 2 is close in mass to h 1 , which can be trivally done, taking into account that Hh 2 h 2 coupling takes the following form: We discuss these limits in more details below, together with the Dark Matter (DM) constraints.

Dark Matter relic density and direct/indirect detection
The results from Planck [76,77] (see also WMAP [78]) have further decreased the error on the already quite precise measurement of the dark matter relic density, Ω DM h 2 : In the i2HDM model, the lightest inert scalar, h 1 , is stable and will contribute to this relic density. In our study we take the upper limit on Ω DM h 2 as a really hard one, since it would lead to DM overabundance. At the same time we do not exclude i2HDM parameter space where h 1 underabundance takes place, since there could be other sources of DM coming from an additional new physics sector. We have evaluated Ω DM h 2 values using the micrOMEGAs 2.4.1 package [79,80,81] since it directly reads the model files in CalcHEP format. Fig. 4(a) shows the relic density in case with quasi-degenerate h 1 , h 2 and h + masses, GeV. This is qualitatively different from the case with a non-negligible mass splitting as illustrated in Fig. 4 One should also note that scenarios with positive or negative λ 345 values of the same magnitude are qualitatively similar, except for the effect of interference (see dashed versus solid curves in Fig. 4). One can observe the following effects and features of the model in Fig. 4: • The red-shaded region in Fig. 4(a) is excluded by LEP data, since in this region W and Z boson decay channels to the light inert scalar would take place. Respectively, the effect of resonant co-annihilation, h 1 h 2 → Z and h 1 h + → W + , can be seen in this region in the first two dips for M h 1 ∼ 40 and 45 GeV.
• In the case of larger M h 1 − M h 2 mass split ( Fig. 4(b)), this effect disappears since the sum of the masses • The sharpest dip in the Ω DM h 2 distribution around 65 GeV corresponds to the DM annihilation through the Higgs boson, h 1 h 1 → H and it is present in both cases.
• At higher masses, we observe a wider and more shallow dip at around 80-90 GeV from h 1 h 1 → W + W − and h 1 h 1 → ZZ channels which are merged together.
• Finally the last dip around 125 GeV, corresponds to the reduction of the DM relic density due to the opening of the h 1 h 1 → HH annihilation channel. This dip takes place only for large values of λ 345 , which provide a high enough rate for the h 1 h 1 → HH process via an s-channel Higgs boson.
• The pattern of these last three dips is the same for larger mass split scenario, presented in Fig. 4(b). One should also note the interference effect in both scenarios, which is sensitive to the sign of λ 345 and appears in this region as a result of the positive or negative interference of the s-channel Higgs boson exchange diagram and the rest of annihilation diagrams.
• One can also observe the qualitative difference in the asymptotic behaviour of the DM relic density for small and large M h 1 values for the different values of ∆M . In the ∆M = 1 GeV case with M h 1 below 65 GeV, the effective co-annihilation of the inert scalars keeps the DM density always below the Planck limit, while in the case of ∆M = 100 GeV, DM coannihilation is suppressed and the relic density is equal or below the experimental limit only for large values of λ 345 (λ 345 0.3) which are excluded by LHC limits on the invisible Higgs decay.
• For M h 1 well above 65 GeV, co-annihilation effects become less important in comparison with h 1 h 1 annihilation into vector bosons, which opens in this region. For this annihilation process the quartic couplings of DM with longitudinal vector bosons, , the respectively small values of the h 1 h 1 V L V L quartic couplings generate a low h 1 h 1 annihilation cross section, σv , which decreases with increasing M h 1 . Eventually this leads to comparatively high values of Ω DM h 2 (which increases with M h 1 both due to the decrease of σv as well as the increase of the DM mass) which reaches the PLANCK limit for large enough M h 1 as one can see from Fig. 4(a). On the contrary, for large values of ∆M c and/or ∆M 2 , the mass splitting generate a high rate for h 1 h 1 annihilation into vector bosons, which increases with increasing M h 1 . This generates a DM density below the experimental limit even for large values of M h 1 . In this scenario the potential increase of Ωh 2 due to an increase of the DM mass is compensated by the respective increase of σv and leads to an approximately flat behaviour of Ωh 2 versus M h 1 in the 100-1000 GeV range. This makes the asymptotic behaviour of the DM density versus M h 1 qualitatively different for ∆M = 100 GeV compared to ∆M = 1 GeV, as one can see from Fig. 4(b). These two scenarios, with a large and small ∆M c , ∆M 2 mass split respectively, qualitatively cover the whole parameter space of the i2HDM. We have also checked whether the i2HDM parameter space is consistent with the limits from DM direct detection (DD) experiments. We have evaluated the spin-independent cross section of DM scattering off the proton, σ SI , also using the micrOMEGAs package. In Fig. 5 limits from LUX100 are shown by the shaded green area where the left and right frames illustrate the small and large ∆M scenarios in analogy to the previous Fig. 4. To present the results in Fig. 5 we use the re-scaled DD cross section,σ SI = R Ω × σ SI , where the scaling factor R Ω = Ω DM /Ω Planck DM allows us to take into account the case when h 1 only partly contributes to the total budget of DM relic density, thus allowing a convenient way to compare the model prediction with the limits from LUX [82].
One can note the flat asymptotic ofσ SI in Fig. 5(a) for large values of M h 1 , which means that the decrease of the proton-DM scattering cross section σ SI with increasing M h 1 , is compensated by the respective increase of the relic density which one can observe in Fig. 4(a). In Fig. 5(b), on the other hand,σ SI drops with large and increasing values of M h 1 . This can be understood  Figure 5: Rescaled spin independent direct detection ratesσ SI versus M h 1 and the LUX100 constraint. The red-shaded region in the left frame is excluded by LEP data.
by observing from Fig. 4(b) that in this region R Ω = Ω DM /Ω Planck DM constant, and therefore the asymptotic behaviour ofσ SI is the same as for σ SI , which is to decrease as M h 1 increases due to the reducing σv .
A related question is whether the model can be better probed by indirect detection (ID) experiments, i.e. the detection of energetic cosmic rays like e + , γ, p orp, which may be created by the annihilation of h 1 pairs. We have checked that the strongest bounds on the i2HDM parameter space coming from such experiments are set by gamma ray telescopes: both the Fermi-LAT gamma-ray space telescope [83] as well as ground based telescopes. Fermi-LAT is sensitive to gamma rays particularly in the low mass range up to O(100 GeV), but the bounds are not competitive with those coming from DD. This conclusion is also confirmed by studies in Ref. [60].
3 Numerical scan of the parameter space 3

.1 Results for the generic scan
To have a complete picture of the properties of i2HDM in the whole parameter space, we have performed a five-dimensional random scan of the model parameter space with about 10 8 points, evaluating all relevant observables and limits mentioned above. The range for the model parameters of the scan was chosen according to the Eq. 17.
To better delineate the impact of each constraint, we have imposed different cuts sequentially on the parameter space, following the classification below: Cut-1: theoretical constraints on the potential from vacuum stability [Eq. (10)(11) and (16) The results of the scan are presented in Fig. 6 in the form of a colour map of DM relic density, projected on two-dimensional planes: (M h 1 , λ 345 ) in the first, (M h 1 , M h 2 ) in the second, and (M h 2 , M h + ) in the third column, respectively. The four rows reproduce the effect of the progressive application of the 4 Cuts defined above. In Fig. 7 we also present, in the same format, the results of a finer scan, zoomed to the region of low masses, where the range has been restricted to  GeV for the three masses M h 1 , M h 2 and M h + . The latter is the most relevant corner of parameter space for the LHC phenomenology that we will discuss in the next section. Note that the lower bound of λ 345 presented in these plots corresponds to the lowest limits allowed by unitarity, perturbativity and scalar potential constraints (see Fig. 2).
One can see from Figs. 6-7(a) that λ 345 is limited from above, and the dependence which defines the shape of this limit as a function of M h 1 comes from the vacuum stability condition given by Eq. (15). One can also see from Figs. 6-7(a), (b), and analogous figures in the rows below, that the relic density is too high for small M h 1 values and small λ 345 . Therefore, the relic density constraint combined with the LHC Higgs data constraints (limiting the invisible decays of the Higgs) restricts M h 1 to be above 45 GeV, as it can be clearly seen from Figs. 6-7(g) and (h). For example, the range 45 GeV< M h 1 <50 GeV is allowed but it requires h 1 − h 2 co-annihilation and respective mass degeneracy, as one can see from Figs. 6-7(h) and (k). From Figs. 6-7(a), (b) and analogous ones in the rows below, one can see a clear vertical blue pattern of low relic density corresponding to the h 1 h 1 → H resonant annihilation. For M h 1 > M H /2 the pattern of DM relic density follows the pattern of W W , ZZ and HH thresholds presented earlier in Fig. 4.
One can also observe that the effect of Cut-1 plus Cut-2 is quite dramatic: a) Br(H → h 1 h 1 ) < 0.28 and µ γγ = 1.14 +38 Actually the region with M h 1 < M Z 2 is excluded due to the interplay of several constraints. In the M h 1 < M H 2 region with low |λ 345 | (|λ 345 | 0.02) as required by LHC Higgs data, the only possibility for relic density of h 1 to be sufficiently low to satisfy Planck constraints is the h 1 h 2 co-annihilation channel: potentially this co-annihilation could provide low enough relic density HH channel is open for negative λ 345 , the relic density drops substantially below the Planck limit, which makes the rescaling factor low enough to avoid limits from LUX searches. The difference between the positive and negative λ 345 cases is related to the respective positive and negative interference of h 1 h 1 → H → XX channel with non-Higgs-exchange diagrams. This asymmetry between positive and negative λ 345 cases was seen initially in Fig. 4, where the h 1 relic density was presented as a function of M h 1 for different λ 345 values.
In summary, after all constraints given by Cut-1-Cut-4, we found that the parameter space with is completely excluded. Our results agree with the results of previous studies on the i2HDM (see, e.g., [60,29]), however we would like to stress that the general exclusion given by Eq.30 is established here for the first time, to the best of our knowledge. In [29], for example, the authors demonstrate (see Fig.6 and Eq. 18 in [29]) that M h + above M H is excluded from a specific scan.
Here we find that M h + as light as 70 GeV is allowed by all present constraints, while M h 1 and M h 2 are generically allowed to be as light as 45 GeV. One should note that specific regions of the parameter space can be excluded using di-lepton and missing transverse momentum signatures: for example, in a recent study [71] the authors showed that values of the masses below M h 1 50 GeV and M h 2 140 GeV can be excluded using this signature, provided that the mass gap between M h 2 and M h 1 is large enough. However, we find that this parameter space region is already excluded by the upper cut on the relic density (Cut-3), as one can see from Fig. 7(h): for M h 2 > 100 GeV, the entire region M h 1 50 GeV is excluded by Cut-3 combined with previous cuts (including LEPII limits).
We would also like to point to some features of the scan for the region of M h 1 , M h 2 and M h + above 200 GeV, presented in Fig. 6. From Figs. 6(f),(i),(l), one can see that EWPT constraints require a very modest mass split between M h 2 and M h + since this mass split is directly related to values of the M h 2 and M h + couplings to the SM Higgs as well as to the couplings to longitudinal components of the W and Z-bosons. Therefore constraints from S and T parameters leave only a rather narrow corridor in the (M h + , M h 2 ) plane.

Fitting the relic density
In our analysis, we generically allow the DM relic density to be equal or below the Planck constraints, Eq. (29). This is the concept of our approach, following which we assume that, in the case of under-abundance, there should be either additional sources of DM or mechanisms other than the thermal freeze-out that compensate for the DM deficit, such as DM freeze-in scenarios [84].
Keeping this in mind, in our analysis we exclude parameter space only where the relic density exceeds the Planck constraint. It is however interesting to distinguish the parameter space where both the upper and lower Planck limits are satisfied. This parameter space is presented in Fig. 8 for the "full" scan ( M h 1 , M h 2 , M h + ⊂[10 GeV-1000 GeV]) and in Fig. 9 for the "zoomed" scan (M h 1 , M h 2 , M h + ⊂[10 GeV-200 GeV]), respectively. Many interesting features of the i2HDM parameter space arise from Figs. 8 and 9 once the "correct" amount of DM relic density is required. From the (λ 345 , M h 1 ) plane, shown in Fig. 8(a), one can clearly see only two very distinctive regions of M h 1 where the relic density satisfies the Planck limit within 2 standard deviations. The first region at low mass, realised for 53 GeV M h 1 76 GeV, is presented in detail in Fig. 9. It clearly show the presence of various parts with specific physical properties: a) The region for M h 1 < M H /2 consists of 3 thin strips: the two symmetric wings for positive and negative λ 345 , clearly visible in Fig. 9(a), consists of DM annihilation via the Higgs boson exchange; the thin horizontal line for very small λ 345 consists of the h 1 − h 2 co-annihilation region. The latter can also be seen in Fig. 9(b), where it appears as a thin diagonal strip in the low M h 2 part, and one can also see that it actually starts at 54 GeV and extends beyond M H /2 up to about 73 GeV. The width of this strip is defined by the maximum allowed value of ∆M = M h 2 − M h 1 = 8 GeV, above which the parameter space is excluded by LEP di-lepton searches until M h 2 > 100 GeV (see Eq. 22). For ∆M < 8 GeV and M h 1 < 54 GeV on the other hand, the Ωh 2 is below the Planck limit in this region. The upper edge at 73 GeV is defined by the rapid increase of the h 1 h 1 → W W * contribution, which does not require co-annihilation above this mass. The typical M h 2 − M h 1 mass split in the co-annihilation region is 6-8 GeV, required to make the relic density consistent with the Plank limit. This small, but important region of the parameter space, which is realised for λ 345 0 and consistent with all present constraints has been missed in the previous studies.

The high mass region and the LHC sensitivity
We would like to discuss here separately the second region at high mass, where relic density is "just right", is realised for M h 1 490 GeV, as shown in Fig. 8. The most peculiar feature of this region is the high degree of degeneracy of the three inert Higgs bosons, as one can see from Fig. 8(b),(c),(d). Numerically we find that the maximal mass difference among h 1 , h 2 and h + , that we call ∆M max , does not exceed a few GeV.
Remarkably, the mass split is required to be large enough, so that the relic density can reach the lower value of the Planck limit: the increase of the mass split is, in fact, correlated with the One should also note that with increasing DM mass, the required split between h 1 , h 2 and h + increases. At about 20 TeV for M h 1 , the DM relic density constraint together with requirement of unitarity and perturbativity which are saturated by ∆M max 10 GeV, close the i2HDM parameter space.

Dark Matter signatures and rates
The i2HDM exhibits various signatures that are potentially accessible at the LHC. First of all, it provides a mono-jet signature from the pp → h 1 h 1 j process, the Feynman diagrams for which are presented in Fig. 10. For this process, the relevant non-trivial parameter space is one dimensional: it is just the DM mass, M h 1 , since the second parameter, λ 345 , simply scales the production cross section which is proportional to (λ 345 ) 2 for M h 1 > M H /2. One should note that the mediator mass for this signature is the Higgs mass, M H = 125 GeV, thus the Effective Field Theory (EFT) approach is not applicable for this process. Also, the recent limits by ATLAS [86] and CMS [87,88] collaborations are not directly applicable for this process since they have been obtained for a different spin of the mediator and different spin of DM. There is one more process, namely qq → h 1 h 2 g (gq → h 1 h 2 q) (see diagrams in Fig. 11), that can contribute to a mono-jet signature in the special case of a small mass split between h 1 and h 2 . In this scenario h 2 decays to h 1 plus soft jets and/or leptons. The essential parameter space Figure 11: Feynman diagrams for qq → h 1 h 2 g (gq → h 1 h 2 q) process contributing to mono-jet signature.
for this process is the two-dimensional (M h 1 , M h 2 ) plane which fixes its cross section. Besides mono-jets, the i2HDM gives rise to a mono-Z signature, the diagrams for which are presented in Fig. 12. The first diagram scales with λ 345 while the other two are fixed by electroweak interactions. 4 In general, non-trivial interference takes place between the three different topologies represented by each of three diagrams, so this process can not be approximated by a simplified model. However, when |λ 345 | 0.02 with M h 1 < M H /2 the Higgs boson threshold) or |λ 345 | 1 with M h 1 > M H /2 (above the Higgs boson threshold,) we have found that the first diagram is dominant and defines the respective event kinematics. So for these values of λ 345 and M h 1 one can use just a simplified model with the Higgs boson as the mediator to set LHC limits.
One should also note that for values of |λ 345 | below 0.02 the contribution from diagrams which scales with |λ 345 | drops below 1%, and in this case the Z boson will be the only mediator to probe the i2HDM model at the LHC, with the mono-Z processes being the leading signature for this purpose (not only complementary to the mono-jet signature). This signature will be especially pronounced if M h 2 − M h 1 > M Z , such that the cross section of the mono-Z signature is essentially defined by the cross section of the 2 → 2 process, pp → h 1 h 2 → h 1 h 1 Z. The parameter spaces for this process is the two-dimensional (M h 1 , M h 2 ) plane.
The i2HDM could also provide a mono-Higgs signature via gg → h 1 h 1 H and qq → h 1 h 2 H, whose diagrams are presented in Fig. 13 and Fig. 14 respectively. The only mediator for gg → h 1 h 1 H is the Higgs boson, and the respective cross section scales as (λ 345 ) 2 for small values of λ 345 and (λ 345 ) 4 for large values of λ 345 because of the second diagram. On the other hand, the qq → h 1 h 2 H process takes place via either a Z-boson or a h 2 mediator: the first diagram does not scale with λ 345 , while the last two do. Therefore for large λ 345 the (λ 345 ) 2 scaling takes place for qq → h 1 h 2 H process. In fact, the contribution from the second and the third diagrams of Figure 12: Feynman diagrams for qq → h 1 h 1 Z process contributing to mono-Z signature.  Figure 13: Feynman diagrams for gg → h 1 h 1 H process contributing to mono-Higgs signature.
qq → h 1 h 2 H to the total cross section drops below 1% only for a λ 345 < 0.002, below which the process kinematics and the cross section are determined by the first diagram with two Z-boson propagators. Finally, one should mention the production of DM via vector boson fusion, pp → h 1 h 1 jj, the diagrams for which are presented in Fig. 15. Similarly to the mono-Z process, there are three different diagrams with different topologies and mediators which contribute to this process: thus, it cannot be described by just one simplified model. The first and second diagrams, which scale Figure 14: Feynman diagrams for qq → h 1 h 2 H process contributing to mono-Higgs signature. Figure 15: Diagrams for qq → q ( ) q ( ) h 1 h 1 DM production in vector boson fusion process.
with λ 345 (in the second diagram, the Z L Z L h 1 h 1 coupling is proportional to λ 3 + λ 4 − λ 5 λ 3 + λ 4 + λ 5 = λ 345 for small M h 1 − M h 2 which occurs when λ 5 is small), gives the dominant contribution to the pp → h 1 h 1 jj process for λ 345 1, which decouples only for λ 345 1, when the contribution from the first two diagrams drops below the per-cent level. On the other hand the VBF process can be also enhanced by Z L Z L h 1 h 1 , W + L W − L h 1 h 1 couplings in case of large h 1 − h 2 and h 1 −h + splittings respectively. This opens a new perspective for the exploration of the i2HDM model which we plan to perform in the near future.
In Fig. 16 we present a summary of the cross sections versus Dark Matter mass, M h 1 , for all the processes mentioned above which contribute to mono-jet, mono-Z, mono-Higgs and VBF signatures for the LHC@8 TeV and LHC@13 TeV. We have used the following setup for the process evaluation: • the QCD renormalisation and factorisation scales, Q were chosen to be equal to the transverse momentum of the pair of DM particles, i.e. missing transverse momentum, E miss T for all processes.
• the PDF and the strong coupling constant are as provided by the NNPDF23LO (as_0119_qed) PDF set [89].
• for all processes a cut on the minimal value of missing transverse momentum of 100 GeV was applied.
• the VBF cross section has been evaluated with the following additional cuts: • below we present plots and numbers for cross sections (in the text and  is around 20-30% for the tree-level cross sections presented, dominating PDF uncertainties which are below 10%. Presentation and a detailed discussion of these uncertainties is out of the scope of this paper. The cross section for the pp → h 1 h 1 j processes which scales with (λ 345 ) 2 is represented by a solid grey line in Fig. 16 and has been evaluated for λ 345 = 1. The cross section scaling for the gg → h 1 h 1 H processes is a bit more subtle: all diagrams except the second one in Fig. 13 scale as λ 345 , while the second one scales as λ 2 345 . So for large λ 345 the cross section scales as (λ 345 ) 4 while for low values it scales as (λ 345 ) 2 . Moreover one can expect non-trivial interference between these diagrams for negative values of λ 345 . The cross section for this process is presented by solid red line in Fig. 16 and has been evaluated for λ 345 = 1.
The cross section for the GeV, both of which are independent of λ 345 are presented by the grey and blue dashed lines.
When M h 2 − M h 1 = 1 GeV, the processes qq → h 1 h 1 Z (with h 2 being off-shell), qq → h 1 h 2 H and qq → jjh 1 h 1 (VBF) scale partly with λ 345 as discussed above. For λ 345 = 1, a value chosen such that the contribution from diagrams which scale with λ 345 is dominant, the cross section for these processes is presented by blue, red and green dot-dashed lines respectively, while for small values of λ 345 (0.002, 0.02 and 0.02 respectively) below which the contribution from such diagrams is not more than 1%, the cross sections are represented by blue, red and green dotted lines. Using information presented in Fig. 16 one can easily estimate cross sections for other values of λ 345 by recalling those cross sections which are proportional to λ 2 345 . One should stress an importance of the pp → h 1 h 2 j process involving h 2 in the final state and the crucial dependence of the physics, i.e. the signatures and cross sections, on the mass of h 2 : a) in case of a small h 2 − h 1 mass split, this process contributes to a mono-jet signature, which actually has a cross section larger than that of gg → h 1 h 1 j even for λ 345 1 when M h 1 > 75 GeV; b) when the mass split is increased to a few dozens of GeV, the cross section of this process drops, however in this case it produces a SUSY-like signature with soft-leptons and E miss T as studied previously and discuss above; c) when M h 2 − M h 1 > M Z , the cross section of this process is enhanced by about two orders of magnitude providing a mono-Z signature with a rate close to that of the pp → h 1 h 2 g process.
In Table 1 we present six benchmarks (BM) from the i2HDM parameter space together with corresponding observables: DM relic density (Ωh 2 ), spin-independent DM scattering rate on the proton (σ p SI ) accompanied with its ratio to the experimental limit from LUX following re-scaling with the relic density: We also present LHC cross sections for the mono-jet, mono-Z and mono-H signatures discussed above with a E miss T > 100 GeV cut applied. All of these benchmarks are allowed by the present experimental data.
The first two benchmarks have small and medium values of λ 345 and correspond to the scenario when M h 1 is below M H /2, and the mass split ∆M = M h 2 − M h 1 is small. BM1 has a very small value of λ 345 = 10 −4 and is therefore characterised by having a small Br(H → h 1 h 1 ) value and a very low DM direct detection rate, σ p SI , whilst the relic density is consistent with the Planck limit due to co-annihilation. The h 1 h 1 j mono-jet signature rate at the LHC scales with (λ 345 ) 2 and is therefore very low, while the λ 345 -independent h 1 h 2 j signature cross section is about 36.7 fb and 92.4 fb for LHC@8 TeV and LHC@13 TeV respectively.
BM2 differs from BM1 only by the value of λ 345 = 0.027, which has been chosen to be the   Table 1: Benchmarks (BM) from the i2HDM parameter space together with corresponding observables: DM relic density (Ωh 2 ), spin-independent DM scattering rate on the proton (σ p SI ) accompanied with its ratio to the experimental limit from LUX following re-scaling with the relic density: . The corresponding LHC cross sections for mono-jet, mono-Z and mono-H signatures with a E miss T > 100 GeV cut applied. maximum value allowed by the Higgs invisible branching ratio. For this value of λ 345 the h 1 h 1 j mono-jet production rates are 288 fb (LHC@8 TeV) and 878 fb (LHC@13 TeV).
BM3 and BM4 correspond to scenarios where ∆M > M Z with M h 1 below and above M H /2 respectively, with the other parameters chosen such that the relic density is consistent with Planck data. In comparison to BM3, BM4 has a very low h 1 h 1 j production cross section because the SM Higgs boson is produced off mass shell. At the same time the h 1 h 1 Z cross section is of the same order for both benchmarks -6.48 fb and 3.90 fb for LHC@8 TeV, and 17.8 fb and 11.1 fb for LHC@13 TeV respectively.
Finally, BM5 and BM6 represent cases with a small (5 GeV) mass split and M h 1 = 100 GeV.
The only difference in input parameters are that they respectively have a large (λ 345 = 1) and small (λ 345 = 0.002) value of λ 345 . For both benchmarks, the DM relic density is well below the PLANCK limit, and therefore an additional source of Dark Matter is required. Even for BM6 which has a small value of λ 345 , the DM relic density is of the order of 10 −3 because the DM effectively annihilate via h 1 h 1 → V V and h 1 h 1 → HH channels which are open for this value of DM mass and are defined essentailly only by the weak coupling and are independent of λ 345 (the contribution from h 1 h 1 → V L V L process is small because of small h 1 − h 2 mass split, while the contribution from co-annihilation is subdominant for this value of mass split). For both of these benchmarks, the h 1 h 2 j channel which has cross-sections of 6.93 fb (LHC@8 TeV) and 19.1 fb (LHC@13 TeV) looks the most promising. From Table 1 one can see that different mono-X signatures are very complementary for these suggested benchmarks, especially the h 1 h 1 j and h 1 h 2 j processes which are the main focus of the collider study presented below.

Limits from LHC@8TeV and projections for LHC@13TeV
In this section, we examine limits from current 8 TeV and projected 13 TeV LHC data. We concentrate on limits from mono-jet processes, as these are the mono-object signatures with the highest cross sections, as shown in Fig. 16.
For mono-jet signals we consider 2 different process: the first one is pp → h 1 h 1 j in which the only two parameters that contribute to the total cross section are the dark matter mass M h 1 and λ 345 . The second one is pp → h 1 h 2 j in which all the vertices depend only on the gauge constants, and therefore the only two parameters that contribute to the total cross section are the dark matter masses M h 1 and M h 2 , or equivalently M h 1 and ∆M = M h 2 − M h 1 .
In order to calculate the limits from the LHC at 8 TeV, we used the CheckMATE [90,91,92,93,94,95,96,97,98] framework, which allows easy application of implemented search analyses. This tool takes a given sample of Monte Carlo events in the HEP or HepMC format after parton showering and hadronisation (for which we used Pythia-6 [99]), and performs a detector simulation on these events using Delphes-3 [91]. Subsequently CheckMATE can apply any of its pre-programmed and validated analyses to the generated signal events, and uses the resulting efficiencies along with published information, such as the 95% confidence limit on signal count, to produce results from which we can find the cross-section limit placed on our model by each analysis.
The signature of both processes that we consider: pp → h 1 h 1 j and pp → h 1 h 2 j, is that of a high-p T jet and a large missing transverse momentum, E miss T . In the case of pp → h 1 h 2 j, the h 2 will decay via a h 1 and a Z ( ) -boson. When ∆M is very small, the decay products of the Z will generally be too soft to be reconstructed in the detector and therefore in this case pp → h 1 h 2 j will give a mono-jet + E miss T signature. Using CheckMATE and HepMC files created with the i2HDM model implemented in CalcHEP, we calculated the limits given by all of the mono-jet + E miss T analyses currently implemented in CheckMATE [100,101,102,88] (3 ATLAS and 1 CMS analysis).
For both processes considered, we found that the lowest cross section limits for each benchmark point considered were provided by one of the ATLAS mono-jet + E miss T analysis [102], and these are the limits presented in this section. This analysis requires a leading jet with a p T > 120 GeV and |η| < 2.0, and the leading jet p T /E miss T > 0.5. Furthermore, to reduce multijet background where the large E miss T originates mainly from jet energy mismeasurement, there is a requirement on the azimuthal separation ∆φ(jet, p miss T ) > 1.0 between the direction of the missing transverse momentum and that of each jet. A number of different signal regions are considered with increasing E miss T thresholds from 150 GeV to 700 GeV. Full details are available in the ATLAS paper [102]. In order to project these limits for increased luminosity and to 13 TeV, we use Monte Carlo events to estimate the efficiencies for the signal and background at 13 TeV, which is a function of M h 1 and depends on the best analysis signal region for each mass. We make the assumption that the analysis cuts for 13 TeV data will be the same as for 8 TeV data, which does not take into account improvements in the signal to background ratio which would likely occur with new analysis cuts at 13 TeV. Therefore our projected limits will be slightly conservative.   The results for the process pp → h 1 h 1 j are shown for 8 TeV in Fig. 17 (a) with projections to 13 TeV and higher luminosities in Fig. 17 (b). The limits are denoted by the solid lines, whilst the cross sections for the i2HDM for different values of λ 345 are shown by the dashed lines. For M h 1 < M H /2, the maximum allowed value of λ 345 is given by the bound on the Higgs branching ratio to invisible in Eq. (24) (this constraint has not been applied on the (dashed blue) λ 345 = 0.1 curve), whilst when M h 1 > M H /2 the maximum allowed value is calculated using the constraints of Eq. (16). The cross section with this maximum value of λ 345 is denoted by the dashed green line. We see in Fig. 17 (a), that the 8 TeV LHC mono-jet + E miss T searches do not constrain the i2HDM via the pp → h 1 h 1 j process. However at 13 TeV, shown in Fig. 17 Fig. 17 We should note that a similar projection of CMS mono-jet limits [88] at 14 TeV has been studied previously [60], where the projected limits were slightly stronger than in Fig. 17 (b). Their projection was able to limit M h 1 for values of λ 345 as small as λ 345 = 0.01, while we require slightly larger values of λ 345 in order to limit M h 1 . We would like to note that in our paper the limits are based on the fast detector simulations rather than parton level ones used in [60] done for 14 TeV. Taking this into account we consider our results as more realistic projection of the future LHC data potential. In Fig. 18 we provide the limit on λ 345 versus M h 1 for different projected luminosities at the LHC@13TeV. This limit is derived from the analysis presented in Fig. 17 and could be more practical for comparison with limits on λ 345 from different experiments.
For pp → h 1 h 2 j, the results are shown in Fig. 19(a) for 8 TeV and in Fig. 19(b) for 13 TeV: we consider two scenarios with a small (∆M = 1 GeV in blue) and large (∆M = 100 GeV in yellow)  GeV σ(pp →h 1 h 2 j) (fb) at 13 TeV  In this case, it should be emphasised that as the couplings of the relevant diagrams (see Fig. 11) are fixed by the gauge couplings, this limit on M h 1 is independent of all parameters other than ∆M . At 13 TeV, and at higher luminosities, this lower limit on M h 1 in this degenerate mass scenario is improved slightly to 41 GeV, 43 GeV and 55 GeV for 20 fb −1 (solid red), 100 fb −1 (solid magenta) and 3000 fb −1 (solid black) of integrated luminosity respectively, as is shown in Fig. 19 (b). For ∆M = 100 GeV, the production cross section is much smaller and the model is not constrained via mono-jet and E miss T signatures from the pp → h 1 h 2 j process. However, in this region, other collider signatures, such as dilepton + E miss T from the decay h 2 → h 1 Z are available , and will provide stronger limits as studied for example in [71].

Constraints on the Parameter space due the extended limits
Taking into consideration these collider limits, and also adding the projections of the Direct Detection XENON1T experiment we are able to impose the complete set of constraints on the i2HDM parameter space. It is worth stressing that as before, we present the limits using the re-scaled DD cross sectionσ SI = R Ω × σ SI , where R Ω = Ω DM /Ω Planck DM , which allows us to take into account additional sources that could contribute to the amount of DM relic density. The results of the constraints are presented in Fig.20 as the color map of DM relic density in the (M h 1 , λ 345 ) plane together with the projected sensitivity of the LHC@13TeV with 3ab −1 of integrated luminosity using h 1 h 2 j and h 1 h 1 j channels, as well as a projection for the XENON1T experiment 95%CL exclusion regions which are indicated by black, dark grey and light grey colours respectively. In this figure we plot excluded points on the top of the allowed points demonstrating the coloured region of the parameter space which will be always allowed (AA). Fig.20a and Fig.20b present AA parameter space for the combined constraints (black on the top of dark grey and dark grey on the top of light grey) for large and zoomed (M h 1 , λ 345 ) regions respectively, while Fig.20c and Fig.20d present AA regions for separate XENON1T and h 1 h 1 +jet LHC13 constraints respectively. From Fig.20a-c one can see how constraints from the LHC and XENON1T are complementary to each other. One can see that XENON1T will exclude large M h 1 masses for large enough values of λ 345 while the LHC will probe the region of smaller values of λ 345 for M h 1 below the M H /2 threshold using the h 1 h 1 j channel, and will cover all values of λ 345 using the h 1 h 2 j channel for M h 1 below 55 GeV.
Besides the AA region it is informative to find and analyse the region with allowed points on the top of excluded points, therefore the black and grey colours present the region which will be always excluded (AE) or probed by the above experiments. Such region is presented in Fig.21a,b in exact analogy to Fig.20a Figure 21: Projection of the 5D random scan of the i2HDM into (M h 1 , λ 345 ) plane and the expected reach of the LHC@13TeV with 3ab −1 of integrated luminosity using h 1 h 2 j and h 1 h 1 j channels as well as for XENON1T experiment. Allowed points are on the top of the excluded ones, therefore the black and grey colours present the region which will be always excluded (AE) or probed by the above experiments.
When comparing Fig.20a,b and Fig.21a,b -i.e. plots with AA versus AE points -one observes a big difference between the order of the overlay of the excluded and allowed points. This is related to the fact that the Ωh 2 can substantially vary: even for fixed M h 1 and λ 345 values, a large M h + or M h 2 can provide respectively large quartic h 1 h 1 W + L W − L and h 1 h 1 Z L Z L which lead in their turn to an effective h 1 h 1 → W + W + /ZZ annihilation. This brings the relic density down and allows avoidance of XENON1T constraints (once we use DD rates re-scaled to relic density). In the (λ 345 , M h 1 ) plane, for example, these points overlap with the points where the quartic couplings mentioned above are small and the Ωh 2 (and respectively exclusion) is driven only by λ 345 . So the most complete picture comes from the combination of AA and AE plots: the most conservative allowed region comes from AA plots of Fig.20, while the most conservative exclusion region is presented by AE plots of  From Fig.20,21 one can see that imposition of the XENON1T constraint reduces substantially the parameter space, greatly expanding the previous limits imposed by LUX. This effect is not so evident in the other planes, presented in Fig.22 in analogy to Fig.21, because diagrams contributing to the spin-independent cross sectionσ SI for DD is given by the t-channel Higgs boson exchange and therefore is proportional to |λ 345 | 2 .
One should also stress again the importance of the pp → h 1 h 2 + j process, using which one can exclude M h 1 < 55 GeV when M h 2 − M h 1 is small for all values of λ 345 . This is shown clearly with the black dots in the Fig.(22)(b,d) where the (co)annihilation and respective mass degeneracy We have also found the projected limits from colliders of mono-jet signatures and the XENON1T DD experiment for the case where the i2HDM satisfy both the upper and lower PLANCK limits, Eq. (29). In this case, the scattering cross section σ SI is not re-scaled, because we are in the zone with the right amount of DM relic density. The results of these constraints are presented in Fig.(23) as a scatter plot where the red zones represent the right amount of DM relic density. In the first row we show the parameter space of the plane (M h 1 , λ 345 ) in the range ⊂[10 GeV-1000 GeV]). In the second row we present the planes (M h 1 , λ 345 ) and (M h 1 , M h 2 ) but in a low range of mass ⊂[50 GeV-80 GeV]).
As we can see, the incorporation of the DD constraint sets important restrictions on the parameter space. In Fig.(23)(a) however there is still a zone in the upper range of mass that is not ruled out. Also in the low mass range there is a zone which survives the restrictive constraint contained in the range between 55 GeV < M h 1 < 74 GeV, for small values of λ 345 . We zoom into the surviving low mass region, as we can see in Figs. (23)(b,c). Because of the improved limits of the DD experiment, the parameter λ 345 is very sensitive to scattering cross section, which set a limit of |λ 345 | < 0.01 for this range of mass.
The collider signatures of the pp → h 1 h 2 + j process, imposes an exclusion limit for M h 1 < 55 GeV (black dots) at the beginning of the region represented by the thin horizontal strip, for very small values of λ 345 , in Fig.(23)(b) as well as thin strip in the lower part of Fig.(23)(c) corresponding to the h 1 − h 2 (co)annihilation region. The pp → h 1 h 1 + j process imposes an extra constraint on the lower mass zone where the DM annihilates through Higgs boson exchange, and is visible in Fig.(23)(b) where in has a shape of two symmetric wings for negative and positive values of λ 345 . This excludes the M h 1 < 55 GeV region. At the same time, XENON1T will improve this constraint and exclude the M h 1 < 56.5 GeV region.

Concluding remarks
The i2HDM is a clear example of a minimal consistent DM model which is very well motivated by theoretical considerations. At the same time this model could provide mono-jet, mono-Z, mono-Higgs and VBF+E miss T signatures at the LHC complemented by signals in direct and indirect DM search experiments.
The model is implemented into the CalcHEP and micrOMEGAs packages, and is publicly available at the HEPMDB database together with the LanHEP model source, and is ready for a further exploration in the context of the LHC, relic density and DM direct detection.
In this paper we have performed detailed analysis of the constraints in the full 5D parameter space of the i2HDM from perturbativity, unitarity, electroweak precision data, Higgs data from LHC, DM relic density, direct/indirect DM detection and LHC mono-jet analysis as well as implications of experimental LHC studies on disappearing charged tracks relevant to high DM mass region. LHC mono-jet analysis for the i2HDM model has been performed at the fast detector simulation level and provides new results together with limits from disappearing charged tracks at the LHC. Our results on non-LHC constraints are summarised in Fig.6-9 presenting the effect of consequent application of constraints from: Cut-1) vacuum stability, perturbativity and unitarity; Cut-2) electroweak precision data, LEP constraint and Higgs LHC data; Cut-3) relic density constraints and Cut-4) constraints from LUX on DM from direct detection. In this paper we have explored for the first time the parameter space where DM from the i2HDM is underabundant implying an additional source of DM, as well as the parameter space where the model's DM contributes 100% of the DM budget of the Universe to the relic density.
Though in general the parameter space of the i2HDM is 5-dimensional, the parameter space relevant to the LHC mono-jet signature is only 1-2 dimensional, so the model can easily be explored at the LHC. There are two qualitatively different and complementary channels in monojet searches: pp → h 1 h 1 j and pp → h 1 h 2 j, with the second one being relevant to mono-jet signature when the mass gap between h 2 and h 1 is of the order of a few GeV. In case of h 1 − h 2 degeneracy, the rate for pp → h 1 h 1 j will be doubled as one can see from Eq. 28 and this can be easily taken into accounts for the estimation of constraints in the respective region of the parameter space. For a fixed M h 1 , the strength of the first process depends only on λ 345 because the Higgs boson is the only mediator, while the strength of the second process is fixed by the weak coupling since the Z-boson is the only mediator for this process. The last process is important to cover the h 1 h 2 co-annihilation region realised when 54 GeV < M h 1 < 74 GeV where the relic density agrees with the Planck data. The results about this process and about this region are new as far as we are aware. Therefore these two processes complement each other in covering the parameter space: for large values of λ 345 pp → h 1 h 1 j would be the dominant LHC signature, while for small or vanishing values of λ 345 the pp → h 1 h 2 j process will cover additional parameter space as demonstrated in Fig.20-23.
Talking about quantitative results, the LHC has a quite limited potential to probe M h 1 with the mono-jet signature. Even for a projected luminosity as high as 3ab −1 , we have found that the LHC could set a limit on M h 1 up to 83 GeV from the pp → h 1 h 1 j process with the maximal value allowed for λ 345 and M h 1 only up to 55 GeV from pp → h 1 h 2 j for any value of λ 345 , covering just a tip of h 1 h 2 co-annihilation region. Such a low sensitivity of the LHC is related to the similarity between the shapes of the E miss T distribution of the dominant Zj → ννj background and that of the signal which has the same Z-boson mediator, while the DM mass is not very different from M Z /2, which as shown in [15] is the reason for such a similarity in E miss T shape. At the same time, the potential of the LHC using a search for disappearing charged tracks is quite impressive in probing M h 1 masses up to about 500 GeV already at 8 TeV with 19.5 fb −1 luminosity as we have found in our study.
We have also explored the projected potential of XENON1T to probe the i2HDM parameter space and have found that it is quite impressive, confirming results of previous studies. In our study we have presented "absolutely allowed" and "absolutely excluded" points in different projections of the i2HDM 5D space demonstrating different features of the models and the potential of current and future experiments. In general, DM DD experiments and collider searches complement each other: the pp → h 1 h 1 j process is complementary in the region with a large λ 345 coupling where DM DD rates are low because of low relic density re-scaling, while the pp → h 1 h 2 j process is sensitive to the parameter space with low λ 345 where DM DD rates are low because of the low rate of DM scattering off the nuclei.