Observable Heavy Higgs Dark Matter

Dark Matter (DM), arising from an Inert Higgs Doublet, may either be light, below the $W$ mass, or heavy, above about 525 GeV. While the light region may soon be excluded, the heavy region is known to be very difficult to probe with either Direct Detection (DD) experiments or the Large Hadron Collider (LHC). We show that adding a second Inert Higgs Doublet helps to make the heavy DM region accessible to both DD and the LHC, by either increasing its couplings to the observed Higgs boson, or lowering its mass to $360 \gev \lesssim m_{DM}$, or both.


Introduction
The long-awaited Higgs boson, with a mass of m h ≈ 125 GeV, was famously discovered in 2012 by the ATLAS and CMS experiments at the CERN Large Hadron Collider (LHC) [1,2]. Although its properties are in accordance with the predictions of the Standard Model (SM), including Electro-Weak (EW) precision data, it remains an intruiging possibility that the observed Higgs boson, denoted here as h, may just be one member of an extended Higgs sector. A good motivation for the latter is the idea that it might provide a candidate for Cold Dark Matter (CDM).
Although the nature of Dark Matter (DM) is not yet known, according to the Standard Cosmological Lambda-CDM Model [3] it should be a particle which is stable on cosmological time scales, cold, i.e., non-relativistic at the onset of galaxy formation, non-baryonic, neutral and weakly interacting. Various candidates for such a state exist in the literature, the most well-studied being the Weakly Interacting Massive Particles (WIMPs) [4,5,6], with masses between a few GeV and a few TeV. Any such WIMP candidate must be cosmologically stable, usually due to the conservation of a discrete symmetry, and must freeze-out (i.e., drop out of thermal equilibrium) to yield the observed relic density [3] 1 : It is clear that the SM Higgs sector cannot provide a WIMP candidate, since its Higgs boson is unstable. However, it was suggested some time ago that the Higgs sector could be extended by the addition of an extra doublet, which may not develop a Vacuum Expectation Value (VEV), leaving a discrete Z 2 symmetry unbroken [7]. Independently, it was later shown that an extra scalar doublet with zero VEV, odd under a discrete Z 2 symmetry, could yield monojets at hadron colliders while being constrained by DM considerations (the first time to our knowledge that any connection with hadron colliders or DM was made) [8]. This possibility, which became known as the Inert Doublet Model (IDM), has been studied extensively for the last few years (see, e.g., [9,10,11]). Since the IDM involves 1 Inert Doublet plus 1 active Higgs Doublet, we shall also refer to it henceforth as the I(1+1)HDM.
In the IDM, aka the I(1+1)HDM, one extra spin-zero SU (2) L doublet with the same quantum numbers as the SM Higgs doublet is introduced. One of the possible vacuum states in this model involves the first doublet acquiring a VEV, henceforth called the active doublet, while the second doublet does not develop a VEV and is referred to as the inert doublet since it does not take part in EW Symmetry Breaking (EWSB). Since this doublet does not couple to fermions and it is by construction the only Z 2 -odd field in the model, it provides a stable DM candidate, namely the lightest state among scalar and pseudo-scalar Z 2 -odd particles.
The I(1+1)HDM remains a viable model for a scalar DM candidate, being in agreement with current experimental constraints. As of now, there are two regions of DM masses where one can expect viable solutions: a low DM mass region, 53 GeV m DM m W and a heavy DM mass region, m DM 525 GeV. The most recent experimental data, both from direct detection experiments and from the LHC, has reduced the viable parameter space in the low mass region [12,13]. In the heavy mass region, however, where the sensitivity of DM direct detection experiments decreases significantly with increasing DM mass, the DM candidate may escape possible detection in this model.
In a recent paper [14] we studied DM in a model with 2 inert Higgs plus 1 active Higgs doublet, which we referred to as the I(2+1)HDM. In particular we focused on the region of parameter space of the I(2+1)HDM where the DM candidate, the lightest inert scalar, is in the light mass region (m DM m W ). We found that the extended scalar sector can relax the exclusion limits from direct detection experiments, providing a viable DM candidate in a region of parameter space which would be excluded in the I(1+1)HDM. In this paper we study the heavy DM mass region of the I(2+1)HDM. We show that heavy Higgs DM in this model becomes more readily observable as a result of either lowering the DM mass to 360 GeV m DM , or increasing the DM-Higgs coupling, or both, while always maintaining the DM relic density within the required region.
The layout of the remainder of this paper is as follows. In section 2 we review the I(2+1)HDM and focus on a simplified version of the model based on a smaller number of parameters. In section 3 we calculate the relic density in the I(2+1)HDM, discussing the relevant DM annihilation scenarios, including the extended co-annhilating (pseudo-)scalar sector. Section 4 will be focused on new features of the I(2+1)HDM with respect to the I(1+1)HDM in the context of DM phenomenology, including enhanced DM-Higgs couplings and the new mass region 360 GeV m DM 525 GeV as well as on discussing the resulting improved prospects for direct detection. In section 5 we present heavy DM signals via Higgs mediation at the LHC in the I(2+1)HDM which look more promising than in the I(1+1)HDM. Finally in section 6 we conclude the paper.

The scalar potential
It has been shown in [15] that an N-Higgs-Doublet Model potential symmetric under a group G of phase rotations can be divided into two parts; a phase invariant part, V 0 , and a collection of extra terms ensuring the symmetry group G, V G .
We construct our Z 2 -symmetric 3-Higgs Doublet Model potential generated by the group which is of the following form 2 : We shall not consider CP-violation in this paper, therefore, we require all parameters of the potential to be real.
The doublets are defined as where φ 1 and φ 2 are the two inert doublets and φ 3 is the one active doublet which plays the role of the SM Higgs doublet, with h being the SM-Higgs boson and G ± , G 0 are the would-be Goldstone bosons. We assign Z 2 charges to each doublet according to the Z 2 generator in Eq.(2): odd-Z 2 charge to the inert doublets, φ 1 and φ 2 , and even-Z 2 charge to the active doublet, φ 3 . It is clear that the symmetry of the potential is respected by the vacuum alignment (0, 0, v √ 2 ). The neutral fields from the inert doublets could then in principle be DM candidates. These neutral fields are stabilised from decaying into SM particles as a result of the conserved Z 2 symmetry of the potential after EWSB.
To make sure that the entire Lagrangian and not only the scalar potential is Z 2 symmetric, we assign an even Z 2 parity to all SM particles, identical to the Z 2 parity of the only doublet that couples to them, i.e., the active doublet φ 3 . With this parity assignment Flavour Changing Neutral Currents (FCNCs) are avoided as the extra doublets are forbidden to couple to fermions by Z 2 conservation.
The Yukawa Lagrangian of the model is identical to the SM Yukawa Lagrangian, with φ 3 playing the role of the SM Higgs doublet: 2 Note that adding extra Z 2 -respecting terms such as (φ † does not change the phenomenology of the model. The coefficients of these terms, therefore, have been set to zero for simplicity.

Mass eigenstates
The minimum of the potential sits at the point (0, 0, v The mass spectrum of the scalar particles are as follows.
• The fields from the active doublet The fields from the third doublet, G 0 , G ± , h, which play the role of the SM Higgs doublet fields have squared masses: • The CP-even neutral inert fields The pair of inert neutral scalar gauge eigenstates, H 0 1 , H 0 2 , which are rotated by • The charged inert fields The pair of inert charged gauge eigenstates, φ ± 1 , φ ± 2 , which are rotated by into the mass eigenstates, H ± 1 , H ± 2 , have squared masses: • The CP-odd neutral inert fields The pair of inert pseudo-scalar gauge eigenstates, A 0 1 , A 0 2 , which are rotated by R θa = cos θ a sin θ a − sin θ a cos θ a , with tan 2θ a = 2µ 2 12 into the mass eigenstates, A 1 , A 2 , have squared masses: We will refer to (H 1 , A 1 , H ± 1 ) as the fields from the first generation and to (H 2 , A 2 , H ± 2 ) as the fields from the second generation. Each of the four neutral particles could, in principle, be the DM candidate, provided it is lighter than the other neutral states. In what follows, without loss of generality, we assume the CP-even 3 neutral particle H 1 from the first generation to be lighter than all other inert particles, that is: (Note that this choice is arbitrary: if the CP-even particle from the second generation, H 2 , where to be assumed lighter than the other inert states, then H 2 will play the role of the DM candidate.) Assuming the CP-even neutral inert particles are lighter than the CP-odd and charged inert particles puts the following constraints on the parameters: In our DM analysis, we consider cases where the mass alignment is changed, but where H 1 is always the lightest inert state and hence is the DM particle. In the remainder of the paper the notations H 1 and DM particle will be used interchangeably.

Simplified couplings in the I(2+1)HDM
Due to the large number of free parameters in the I(2+1)HDM which makes it impractical to analyse the model in the general case, we focus on a simplified case where parameters related to the first inert doublet are k times the parameters related to the second doublet without introducing any new symmetry to the potential. The motivation for this simplified scenario is that in the k = 0 limit the model reduces to the well-known I(1+1)HDM. We assume no specific relation among the other parameters of the potential. It is important to note that the remaining quartic parameters do not influencethe discussed DM and LHC phenomenology of the model and thus their values have been fixed in agreement with the theoretical constraints discussed in the original paper [14] and compliant with the results on unitarity obtained in [16].
The k = 1 case In this paper we focus on the k = 1 case in the heavy DM mass region 4 . The mass spectrum in this case is simplified to: The quartic couplings in the potential can be written in terms of the masses of the physical particles as: , where g H 1 H 1 h = λ 23 + λ 23 + 2λ 2 is the Higgs-DM coupling. The Feynman rules for this model are presented in Appendix A.
3 Calculating the relic density in the I(2+1)HDM The relic density of the WIMP (identified in our model as the lightest inert scalar H 1 ) is calculated with the assumption that the WIMP was in thermal equilibrium with the SM particles after inflation. Once the rate of DM DM ↔ SM SM reactions becomes smaller than the Hubble expansion rate of the Universe, the WIMP freezes out, i.e., drops out of the thermal equilibrium. After freeze-out the co-moving WIMP density remains essentially constant with the current value estimated by the Planck experiment to be the one already given in Eq. (1).
As mentioned, in the I(2+1)HDM one of the neutral inert (pseudo-)scalar particles play the role of the DM. The relic density of a (pseudo-) scalar DM candidate, S, after freeze-out is given by the solution to the Boltzmann equation: where the thermally averaged effective (co)annihilation cross-section contains all relevant scattering processes of any S i S j pair into SM particles: where Therefore, only processes in which the mass splitting between a state S i and the lightest Z 2odd particle S (H 1 in our case) are comparable to the thermal bath temperature T provide a sizeable contribution to this sum.
In the I(2+1)HDM, the presence of additional inert particles has important consequences in the heavy mass regime. For lighter masses the most important channel for the annihilation of DM particles is the Higgs-mediated process Fig. 12a), as studied in [14]. However, coannhilation with H 2 , A 1 and A 2 may change the results significantly (see Fig. 12b).
For heavier masses the diagrams including one or two virtual gauge bosons, shown in Fig.  13 also contribute to the total annihilation cross-section. Finally, co-annihilation plays an important role in scenarios with multiple particles which are close in mass. This scenario is realised in the I(2+1)HDM for the heavy DM mass region. Particles up to 20% heavier than the DM candidate may influence the DM relic density. Therefore, the co-annihilation diagrams should be included in calculating the effective annihilation cross-section. These diagrams are presented in Figs. 14 and 15 -representing pure gauge channels and coannhilation channels involving the SM-like Higgs particle, respectively.

Relevant co-annihilation scenarios
In the I(2+1)HDM, the strength and importance of coannhilation processes depend on the mass splittings between the inert particles. We define δ A and δ C as the splitting between H 1 and the pseudoscalar and charged state from the first generation, respectively, δ A,C are related to the quartic couplings in the potential, which are constrained by the perturbativity (and unitarity) conditions, i.e., the λ i 's cannot be too large. As a result of this, all particles within one generation will have a similar mass. These masses, however, could have high values because of the (almost) unconstrained quadratic parameters µ 2 2 and µ 2 12 : It is important to stress that, even if bounds on λ i were relaxed leading to larger values of δ A,C , there exist very stringent limits from relic density analysis. Coannihilation must occur at least between H 1 , A 1 and H ± to achieve DM relic density in agreement with the current experimental measurements. This is a pattern followed by all general heavy scalar DM models. In the absence of these co-annihilation channels, the maximum relic density that can be achieved through H 1 H 1 → SM SM (even when H i H j → SM SM is allowed) is of order 10 −3 which is well below the observed value.
The other important mass splitting, ∆, is defined as the mass difference between H 1 and the other CP-even state H 2 ("splitting between doublets"): ∆ is related to the quadratic parameter µ 2 12 through Note that µ 2 12 is not limited by any theoretical constraints -similar to µ 2 2 -and therefore ∆ can in principle be very large. So, unless ∆ is forced to be small by limits put on µ 2 12 , one should also consider a case where the second doublet is decoupled from the first, leading to a scenario which was not listed in [14]. Therefore, in the very heavy mass region one can consider: • Case G: with small δ A , δ C , ∆, where all inert particles are close in mass and coannihilate with each other.
• Case H: with small δ A , δ C and large ∆, where the second generation is effectively decoupled from the first generation and does not influence relic density calculations. In this case, the relevant diagrams are the (co)annihilation channels between fields from the lighter generation only, Tab. 1 summarises the two possible scenarios for relic density studies.
Coannhilation between H 1 , H 2 is not efficient enough and DM density is below experimental bounds.

Large ∆ Case H is realised
There are no co-annihilatin channels open and gauge annihlation reduces DM relic density effectively below the experimental bounds.

The gauge limit
To illustrate the difference between cases G and H, let us first consider the gauge limit in both scenarios, which is the limit where all quartic couplings λ i are set to zero. Therefore, all scalar self-couplings, including the DM-Higgs coupling, are removed in this limit. As a result δ A,C = 0, leading to degenerate H 1 , A 1 , H ± 1 states. Note that this limit is excluded by results of direct detection experiments, nevertheless, it is an interesting limit to study as it represents the main difference between cases G and H. In this limit H 1 annihilates solely through the gauge annihilation channels presented in Fig. 14.
Non-zero λ i will lift this degeneracy and, at the same time, reduce the effective annihilation cross-section for a given mass. Therefore, the gauge limit corresponds to the minimum value of m H 1 , for which it is possible to obtain a proper relic density for any value of the Higgs-DM coupling.
These results are presented in Fig. 1 for the two scenarios: case G where all particles have degenerate mass (in the δ A,C , ∆ → 0 ⇔ λ i → 0, µ 2 12 → 0 limit) and case H with large ∆ where the second generation is decoupled from the first generation and all particles in the first generation are degenerate in mass (in the δ A,C → 0 ⇔ λ i → 0 limit).
It is clear that, for a given mass of m H 1 the destructive interference between an increased number of coannhilation diagrams in case G leads to a reduced cross-section, i.e. larger DM relic density with respect to case H. The important consequence of this interference is that for case G it is possible to obtain proper relic density for smaller masses of DM candidate in comparison to case H. Note also that case H behaves like the I(1+1)HDM (the Inert Doublet Model) in this limit. This similarity in behaviour will be repeated as we will show in the following sections.

The benchmark points
With non-zero scalar couplings and mass splitting more diagrams contribute to the coannihilation of DM -all diagrams shown in Figs. 14 and 15 contribute to the total annihilation cross-section. Here we present two benchmarks points, two sets of parameters, for which we have studied the DM relic density: • For case G with δ A = δ C = 1 GeV and ∆ = 1 GeV Here all inert particles have similar masses and therefore can co-annihilate with each other. The degeneracy between charged and "pseudo-scalar" particles is allowed and doesn't lead to any unacceptable results. The important degeneracy which must be avoided is H 1 -A 1 degeneracy leading to the scattering through the Z boson which is tightly constrained by direct detection experiments and puts a lower limit on δ A .
• For case H with δ A = δ C = 1 GeV and ∆ = 100 GeV Here the second generation of inert scalars is significantly heavier than the first one. Within each generation, however, particles are almost degenerate.
Note that there are certain differences between cases G and H. In case H, the heavier generation of inert particles is decoupled from the first generation particles and does not influence the relic density calculations. The model in this case resembles the I(1+1)HDM.
Furthermore, in case G the Higgs-DM couplings which result in a relic density in agreement with experiment are larger in comparison to case H for the same DM mass. This difference is explicit in Fig. 2.
Planck-3σ This leads to the fact that for case G we can obtain viable relic density values for m DM much smaller than in case H (or the I(1+1)HDM) in which the minimal value of m DM resulting in DM relic density in agreement with Planck limits is m H 1 ≈ 525 − 535 GeV. In case G (with ∆ = 1 GeV), however, the DM mass can be as low as ∼ 375 GeV. This result is shown in Fig. 3 which represents relic density plots for cases G (left) and H (right). Note that in case G the minimum m H 1 which touches the lowest acceptable relic density limit (the green solid line) is 375 GeV (for a given ∆ of 1 GeV), whereas in case H this minimum value is 525 GeV (the solid red line). Fig. 4 is meant to represent the same benchmark points as in Fig. 3, in the m H 1 -g H 1 H 1 h plane. The bands correspond to proper relic density in agreement with Planck measurements in case G (for an exemplary ∆ = 1 GeV, δ = 1 GeV) in red and case H (for an exemplary ∆ = 100 GeV, δ = 1 GeV) in red. Note that, for the a given DM mass (and same δ), the Higgs-DM coupling in case G is much larger than in case H.    Figure 4: Relic density bands in agreement with Planck measurements in case G (for an exemplary ∆ = 1 GeV, δ = 1 GeV) in red and case H (for an exemplary ∆ = 100 GeV, δ = 1 GeV) in red. Note that for the a given DM mass, the Higgs-DM coupling in case G is much larger than in case H.
Since δ A and δ C are related to the quartic parameters λ i 's, they are constrained from unitarity bounds and are required to be small. However, regardless of any limits on λ i from unitarity, relic density studies show that both δ A and δ C must be relatively small to allow for co-annihilation between particles which is crucial in the heavy DM mass region. The upper limit follows the following rough rule: co-annihilation effects take place when the mass difference between co-annihilating particles is of the order of 20% of their mass.
The lower bound on δ A comes from direct detection experiments where a degeneracy between H 1 and A 1 leads to the scattering through the Z boson which is tightly constrained. Further, the above limits in Eq. (23) are in agreement with LEP searches for exotic particles. Finally, below δ A,C ∼ 0.1 GeV there is no visible difference in the results.
As mentioned before, a large difference between δ A and δ C violates relic density limits. In the cases we study below, we have set δ A = δ C = δ for simplicity 5 . The mass splitting between the two generations, ∆, is proportional to and therefore unconstrained by unitarity. The maximum value for ∆ is proportional to the arbitrary maximum value chosen for µ 2 12 . In general, we allow for ∆ to vary in the following region 100 keV ∆ 200 GeV.
For large ∆ values, ∆ 20 − 50 GeV (depending on m H 1 , since the story-changing mass splitting is roughly 20%m H 1 ), co-annihilation effects between the two generations are not strong enough to compete with the standard (co)annhilation between H 1 , A 1 , H ± 1 , in which case the I(2+1)HDM acts just like the I(1+1)HDM. The second generation is effectively decoupled from the first generation and does not influence DM phenomenology. Therefore, scenario H is realised for: ∆ 20 GeV ⇒ scenario H The exact value of ∆, when above ∼ 20 − 50 GeV does not make any significant difference in the relic density calculations.
For small values of ∆, the co-annihilation effects between all particles are important. The smaller ∆ is, the more relevant particles from the second generation are for DM studies: for ∆ ≈ 1.5 GeV the relative contribution to relic density calculation coming from particles from the lighter generation to the heavier generation is 70%-30%. For ∆ ≈ 0.0001 GeV this relation is 50%-50%. The case G, is therefore realised when ∆ varies in the following window 0.0001 GeV ∆ 20 GeV ⇒ scenario G The closer ∆ gets to this upper limit, the weaker the coannnihilation effects and the more scenario H is realised. Finally, notice that, for our studies, the Higgs-DM coupling, g hH 1 H 1 , is kept within the |g hH 1 H 1 | < 1 range. Fig. 5 illustrates the effect of changing ∆ on the relic density. In all four plots m DM has been set to 400 GeV as the value of ∆ changes, 0.1 GeV in the top left plot, 1 GeV in the top right plot, 5 GeV in the bottom left plot and 10 GeV in the bottom right plot. In each plot different colours represent different δs. For small values of ∆ (0.1 GeV), the H 1 -H 2 co-annihilation leads to viable relic density values even for large δ (i.e. the H 1 , A 1 , H ± co-annihilation is absent). For large values of ∆ (10 GeV) H 1 -H 2 co-annihilation does not exist and even small values of δ cannot compensate this lack, thus, the relic density is below the acceptable limit.
In Fig. 6 the value of δ is set to 0.5 GeV for m DM = 400 GeV in the left plot and m DM = 550 GeV in the right plot. In each plot, different colours correspond to changing ∆s. 5 Cases with δ A = δ C do not lead to any new phenomenology and in fact the region of validity for g H1H1h decreases as we increase the difference between δ A and δ C .  Note that in the left plot only small values of ∆ lead to viable relic density values, which is where case G is realised. In the right plot, small values of ∆ correspond to case G and large values of ∆ correspond to case H, and they all lead to acceptable relic density values. Note that for ∆ 50 GeV all curves correspond to the same value.

Changes in m DM
Here we describe several sub-regimes where the DM mass can vary with important characteristics.
• In the region m DM 360 GeV neither scenario H nor G results in viable relic density values 6 . This lower limit can be reached in case G by very specific points in the parameter space: (a) when the mass

Direct detection
Direct detection experiments, which are mostly designed to hunt for the standard EW-scale WIMP, are the most sensitive to DM masses of the order of 100 GeV. For heavier DM masses the sensitivity of these experiments drops significantly. The most recent LUX results set the limit of the DM-nucleon scattering cross-section to be σ DM −N ≈ 10 −8 pb for DM masses ≈ 500 − 1000 GeV [17].
In the I(2+1)HDM, similar to the other scalar DM models, DM candidate can be detected through elastic scattering on nuclei through the exchange of a Higgs particle. Therefore, σ DM −N will depend on the value of DM-Higgs coupling, and the DM mass. These results are presented in Fig. 8 where the shaded region corresponds to the probed phase space of the I(2+1)HDM for various choices of ∆ and δ (as shown in Fig. 7), all of which have relic density in agreement with Planck measurements. Also, results for the benchmark points studied in section 3 are presented explicitly. We found that the current experimental limits do not constrain the heavy DM mass region, neither for case H nor for case G, even though the Higgs-DM coupling is larger in the latter case.
Recall that for certain choices of ∆ and δ one can obtain proper DM relic density for g H 1 H 1 h ≈ 0. This leads to a strongly suppressed scattering cross-section, which may not be detected as it lies within the coherent neutrino-nucleus scattering regime [18]. Fig. 8 also shows a limit from the future XENON1T experiment [19], with a proposed sensitivity of the order of 10 −9 − 10 −10 pb (dashed black line). We expect the next generation of DM detectors, such as XENON1T, to be able to test a large portion of the parameter space of the heavy DM in the I(2+1)HDM for m H 1 1 TeV. As a final note to this subsection, we should like to mention that a viable intermediate DM mass region, m W m DM 160 GeV, regarding relic density studies, has been found in the I(1+1)HDM. The correct relic density in this region is obtained due to cancellations between different diagrams contributing to DM annihilation into gauge bosons (W + W − and ZZ). In [20] it was shown that this scenario is realised if the inert particles, in particular the charged scalar, are heavy enough, ∼ 300 − 500 GeV. A relatively large DM-Higgs coupling is also required for the DM in this mass region to stay within relic density limits, however, this large Higgs-DM coupling is excluded by LUX. Similarly, we did not find any solutions in the medium mass region with viable relic density and in agreement with direct detection experiments in the I(2+1)HDM.

Heavy Higgs DM at the LHC in the I(2+1)HDM
A scalar DM candidate is a stable particle with limited interactions with all SM particles and therefore it cannot be directly detected at the LHC. However, its presence can influence the detectable properties of SM particles. One way to ascertain the influence of DM candidate on the properties of a Higgs particle is to look at the Higgs invisible decays, h → SS, where S is a scalar DM candidate with mass below m h /2. Invisible decays of the SM-like Higgs particle in the I(2+1)HDM were studied in [14,21], where limits for the mass of a light DM candidate combined with Planck limits for the relic density measurements provided constraints comparable or stronger than those from direct detection experiments.
Another strategy, useful for a heavy DM particle, is to look for a high p T monojet or a two jet/two lepton signal, accompanied by a large missing transverse energy E T . The monojet signature in the I(2+1)HDM, pp → H 1 H 1 + jet, corresponds to h coupling to an invisible pair of DM particles (yielding the large E T ) with produced in association with an energetic quark or gluon jet. The following processes are considered in our analysis.
1. gg → gH 1 H 1 (Fig. 16) via a triple gluon and a hgg effective vertex. Note, that the hgg effective vertex in the I(2+1)HDM is the same as in the SM, as the Higgs production here is not modified by presence of additional scalar states. This is the dominant contribution to the monojet process, as the gluon fusion is an enhanced production mechanism for the Higgs particle.
3. qg → qH 1 H 1 (Fig. 18), where q = u, d, c, s, b. The dominant contributions here come from gb → H 1 H 1 b with the Higgs boson radiated off of the b quark legs (Fig. 18a) and qg → qH 1 H 1 t-channel via a gqq tree-level vertex and the hgg effective coupling (Fig.  18b).
Note, that all above processes contain the h → H 1 H 1 vertex, therefore a strong dependency on the g H 1 H 1 h coupling is expected, i.e., a significant difference between scenarios G and H, as discussed in the previous section.
For the studies of pp → H 1 H 1 + 2jets we have considered the Vector Boson Fusion (VBF) process of the form q i q j → H 1 H 1 q k q l , with q = u, d where a pair of DM particles (with large E T ) is produced by the neutral (Fig. 19) or charged (Fig. 20) VBF processes, either directly or mediated by the Higgs particle or another neutral scalar.
We have also considered the Higgs-Strahlung (HS) processes of the form q iqj → V * H 1 H 1 where a pair of DM particles is radiated off the Z or W ± boson leg ( Fig. 21 and Fig. 22, respectively), either directly or mediated by the Higgs particle or another neutral scalar.
Notice, that in the dijet searches only one diagram out of each set depends on the g H 1 H 1 h , therefore we expect smaller differences between scenarios H and G than in the monojet searches. The strength of the other diagrams is set by the gauge interactions.
Also, it is important to stress that, given our initial choice of parameters, i.e. introducing the k = 1 relation between the doublets (see Eq.(12)), we have limited the number of possible diagrams, because vertices of the type ZH i A j and W ± H i H ∓ j (i = j) are absent when k = 1. Relaxing this initial assumption would in principle not only influence the evolution of DM relic density, but could also lead to a possibly stronger difference between scenarios G and H in the dijet analysis.
In the following subsections we present results for the monojet and dijet analysis, for the 14 TeV LHC. The following selections were used.
1. For the monojet searches, we require the following cuts on the transverse momentum of the jet, p j T , and the pseudo-rapidity of the jet, η j , p j T > 20 GeV and |η j | > 2.5 2. For the dijet searches, we require the following cuts on the invariant mass of the two jets, M (j, j), and the difference between the pseudo-rapidity of the forward and backward jet, M (j, j) > 700 GeV and |η j Calculations were done with the aid of LanHEP [22] and CalcHEP [23] packages.

Monojet results
In Fig. 9 results for monojet signals of scenarios G (δ A = δ C = 1 GeV, ∆ = 1 GeV) and H (δ A = δ C = 1 GeV, ∆ = 100 GeV) are shown. For comparison, we also present results for the I(1+1)HDM, with δ A = δ C = 1 from [12]. The DM-Higgs coupling (defined as λ 345 in [12]) is the same as g H 1 H 1 h in scenario H for equal DM masses, therefore the monojet diagrams in case H and in I(1+1)HDM are identical. Scenario G, which corresponds to much larger Higgs-DM couplings compared to that of scenario H or the I(1+1)HDM, results in a significantly larger cross-section in the monojet process. Also the special features of the model are more visible in this process. For masses up to m H 1 ≈ 450 GeV we observe a rise in the cross-section connected to an opening of the phase space combined with an increasing Higgs-DM coupling. After that peak, the cross-section decreases with increasing DM mass regardless of the rising of g hH 1 H 1 .
Notice that for the lower masses, the difference between scenario G and H is significant, as every diagram involved in the monojet process contains the hH 1 H 1 vertex, whose coupling differs significantly in cases G and H. Notice the low end of the allowed mass region in case G with a large cross-section for the mass region which is not even accessible by case H or the I(1+1)HDM. As the DM mass grows, results for both cases get closer together, stabilising for the very heavy mass region in the decoupling limit.  19 and 20 for diagrams involved in the neutral and charged VBF processes, respectively). We can still observe some differences in the lower mass range. The cross-section for case G is generally larger (as is the g hH 1 H 1 strength) than in case H. Also, charged channels have slightly larger cross-sections than the neutral ones, since the cross-section for producing the W ± boson is larger than the cross-section for producing the Z boson. For the heavier masses all results, for both scenarios G and H, as well as for charged and neutral channels, tend to converge.

HS results
HS signatures depend on the W and Z decay patterns. While at the LHC, leptonic signatures are preferred, hadronic ones are also possible. The latter potentially interfere with the VBF topologies, but the effect is small so that we can safely ignore it here.
The results of the (on-shell) HS process cross-sections in terms of the DM mass are presented in Fig. 11. It is clear that the general picture is different from the VBF studies. The largest cross-section again appears in the lower mass region where only case G can be realised. Similarly to the VBF case, the charged channels have larger cross-section compared to the neutral channels since the W ± production in theses processes has a larger cross-section than the Z boson production. All cross-sections in neutral and charged processes in both cases G and H converge to a similar value as m DM approaches very large values. Note that the cross-section in both G and H are similar in the region of DM mass above 525 GeV, i.e., where both cases can be realised (unlike in the VBF processes). This similarity is the result of the fact that the difference in the Higgs-DM coupling does not translate into a difference in the cross-section between cases G and H. To explain this similarity let us focus, e.g., on the neutral VBF and HS processes ( Fig. 19 and Fig. 21) In the neutral VBF process, out of all the involved diagrams (Fig. 19a,b,c) there is only one diagram, Fig. 19a, which depends on the Higgs-DM coupling. The cross-section of this diagram (σ h ) for a given m DM relative to the cross-section of all three diagrams involved (σ tot ) is σ h /σ tot = 0.1445 for case G and σ h /σ tot = 0.2416 for case H. Recall that the main difference between cases G and H is that for a chosen m DM the Higgs-DM is larger in case G than in case H. We conclude that the Higgs-mediated diagram plays a much more important role in case H than it does in case G. As a result, even though the g hH 1 H 1 coupling is much smaller in case H, the total cross-section does not fall far below the total cross-section in case G, which is depicted in Fig. 10. Now, let us consider the HS neutral processes (Fig. 21a,b,c). Again, only one diagram, Fig. 21a, depends on the Higgs-DM coupling. Repeating the procedure above, we calculate the relative cross-section of this one diagram (σ h ) relative to the cross-section of all diagrams involved (σ tot ). For a given m DM , we obtain σ h /σ tot = 0.9736 in case G and σ h /σ tot = 0.9708 in case H. So, this diagram plays only a slight role in case H compared to case G. We therefore conclude that the difference in g hH 1 H 1 coupling does not play an important role in distinguishing case G and H in the HS processes. Thus, we do not expect to see a difference between cases G and H, which is depicted in Fig. 11, with the important exception that case G allows for the heavy Higgs DM mass to be below 525 GeV, making its cross-section significantly larger.

Conclusion
We have calculated the relic density for heavy Higgs DM in the I(2+1)HDM and shown that the prospects for its discovery at both DD experiments and the LHC are significantly enhanced as compared to the I(1+1)HDM, where the heavy Higgs DM particle must have a mass above about 525 GeV and is weakly coupled to the observed Higgs boson. Adding a second inert Higgs doublet helps to make the heavy Higgs DM region accessible to both DD and the LHC, by either increasing its couplings to the observed Higgs or lowering its mass to 360 GeV m DM , or both. In particular we have presented LHC signatures of the I(2+1)HDM in the monojet, VBF (dijet) and HS processes and shown that the prospects for heavy Higgs DM discovery are significantly brighter for all channels.
In DD experiments, although the standard values of annihilation cross-section for the heavy Higgs DM masses in the I(2+1)HDM are well below current experimental exclusion limits for DM decaying into pairs of gauge bosons or fermions, the prospects for a future DD discovery remain open due to the complementary nature of collider vs cosmological limits and the fact that the DD cross-sections are higher than in the I(1+1)HDM.
Turning to indirect detection signatures of the I(2+1)HDM, there is the possibility of internal bremsstrahlung in the processes of H 1 H 1 → W + W − γ, generated through the exchange of any of the two charged scalars H ± 1,2 . It was shown that one can expect such signatures in the I(1+1)HDM [24], mediated by a charged scalar in the t-channel, which would correspond to scenario H considered in the I(2+1)HDM. In principle, the signal could even be stronger for scenario G, as the scalar couplings are enhanced.
Finally we comment on the observed h → γγ channel where, in the I(1+1)HDM, only in the heavy DM mass region are both proper DM relic density and (minimal) enhancement in the h → γγ channel realised. By contrast, in the I(2+1)HDM there exist two charged scalars, H ± 1 and H ± 2 , which contribute to the h → γγ loop which may enhance the rate for a wide range of parameters.
In conclusion, adding a second inert Higgs doublet significantly improves the prospects for observability of heavy Higgs dark matter in future experiments both underground and at the CERN LHC. nanced in part through the NExT Institute and from the STFC Consolidated ST/ J000396/1. He also acknowledges the H2020-MSCA-RICE-2014 grant no. 645722 (NonMinimalHiggs). VK's research is financially supported by the Academy of Finland project "The Higgs Boson and the Cosmos" and project 267842.

C Feynman diagrams for the LHC analysis
Here we present the DM (co)annihilation diagrams which play a role in our LHC studies.