Dark Matter with Two Inert Doublets plus One Higgs Doublet

Following the discovery of a Higgs boson, there has been renewed interest in the general 2-Higgs-Doublet Model (2HDM). A model with One Inert Doublet plus One Higgs Doublet (I(1+1)HDM), where one of the scalar doublets is"inert"(since it has no vacuum expectation value and does not couple to fermions) has an advantage over the 2HDM since it provides a good Dark Matter (DM) candidate, namely the lightest inert scalar. Motivated by the existence of three fermion families, here we consider a model with two scalar doublets plus one Higgs doublet (I(2+1)HDM), where the two scalar doublets are inert. The I(2+1)HDM has a richer phenomenology than either the I(1+1)HDM or the 2HDM. We discuss the new regions of DM relic density in the I(2+1)HDM with simplified couplings and address the possibility of constraining the model using recent results from the Large Hadron Collider (LHC) and DM direct detection experiments.


Introduction
The ATLAS and CMS experiments at the Large Hadron Collider (LHC) have found evidence of a Higgs scalar with a mass of m h ≈ 125 GeV, which is in good agreement with earlier predictions performed using Electro-Weak (EW) precision data [1,2]. Further studies are needed in order to determine whether this particle belongs to the Standard Model (SM) or to one of its extensions. However, so far, there are no reports of detection of physics Beyond the SM (BSM), neither by discovery of new particles, nor by any significant deviation from the SM prediction of the Higgs signal strengths, and strong bounds are set for the most common BSM models.
On the other hand, new physics is expected for various theoretical and experimental reasons. One of the most important is the existence of Dark Matter (DM), stable on cosmological time scales, cold, i.e., non-relativistic at the onset of galaxy formation, non-baryonic, neutral and weakly interacting component of the Universe [3]. Strong premises for its existence come from the galactic, cluster and horizon scales, making the modified-gravity based explanations of the observed phenomena less likely. Various candidates for such a state exist in the literature, the most well-studied being the Weakly Interacting Massive Particles (WIMPs) [4,5,6].
WIMP's mass may change roughly between a few GeV and a few TeV, and the annihilation cross section is of approximately weak strength. The relic density of WIMPs is calculated with the assumption that they were in thermal equilibrium with the SM particles after inflation. Once the rate of reactions DM DM ↔ SM SM becomes smaller than the Hubble expansion rate of the Universe, the WIMPs freeze-out, i.e., drop 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 [3]: WIMPs are usually stable due to the conservation of a certain discrete symmetry. In case of the most-studied candidate in Supersymmetric (SUSY) models, neutralino (a Majorana fermion), it is the R-parity related to the imposed R-symmetry [7,8]. Bosonic candidates appear in the models with Universal Extra Dimensions (UED) and are made stable by the KK-parity, the remnant of momentum conservation in the extra dimension [9,10]. One could also consider the scalar candidates, stabilized, for example, by the conserved Z N discrete symmetry in the scalar potential, see, e.g., [11,12,13,14,15,16,17,18].
One of the simplest models that provide a scalar DM candidate is the model with One Inert Doublet plus One Higgs Doublet (I(1+1)HDM) 1 , proposed in 1976 [13], and which has been studied extensively for the last few years (see, e.g., [14,16,17]). In this model one SU (2) W doublet with the same quantum numbers as the SM Higgs doublet is introduced. One could simply extend the I(1+1)HDM by introducing an extra inert SU (2) W doublet with the same quantum numbers as the SM Higgs doublet, resulting in a 3-Higgs-Doublet Model (3HDM). As the next simplest example beyond 2HDMs, which has been extensively studied in the literature, 3HDMs are very well motivated. Furthermore, all possible finite symmetries in 3HDMs have recently been identified [25].
3HDMs may address the problem of the origin and nature of the three fermion families. Indeed it is possible that the symmetry of the three Higgs doublets could describe the symmetry of the three families of quarks and leptons. In a recent paper [26], we studied symmetric 3HDMs and derived the conditions under which the vacuum alignments (0, 0, v 3 ), (0, v 2 , v 3 ) and (v 1 , v 2 , v 3 ) are minima of the potential. Here we focus on the alignment (0, 0, v 3 ), which is of particular interest because of its I(1+1)HDM similarity and the absence of Flavour Changing Neutral Currents (FCNCs) 3 .
In this paper, we study a model with Two Inert Doublets plus One Higgs Doublet (I(2+1)HDM). The two inert doublets are Z 2 -odd and the active doublet is Z 2 -even which plays the role of the SM Higgs doublet. The I(2+1)HDM may be regarded as an extension to the I(1+1)HDM. In this scenario the Z 2 -odd particle content is doubled with respect to the I(1+1)HDM, and so new possibilities of DM (co)annihilation appear. One can have up to six (co)annihilating states, which introduce a very different behaviour with respect to models with fewer number of states in the inert sector.
The layout of the paper is as follows. In section 2 we construct the Z 2 -symmetric I(2+1)HDM potential and study the (0, 0, v 3 ) vacuum point. In section 3 theoretical and experimental constraints on the parameters of the model are derived and presented. In section 4 we list all phenomenologically viable DM (co)annihilation scenarios in our model. In section 5 we study in detail a simplified version of the I(2+1)HDM with reduced number of parameters and present relic density plots in different scenarios. We finally draw our conclusions in section 6. The Feynman rules are presented in the Appendix.

Constructing the I(2+1)HDM potential
In a general N-Higgs-Doublet Model (NHDM), the scalar potential which is symmetric under a group G of phase rotations can be written as where V 0 is invariant under any phase rotation and V G is a collection of extra terms ensuring the symmetry group G [28]. The most general phase invariant part of the I(2+1)HDM potential has the following form: Constructing the Z 2 -symmetric part of the potential depends on the generator of the group. The Z 2 generator which forbids FCNCs and is respected by the vacuum alignment (0, 0, v) has the following form The terms ensuring the Z 2 group generated by g are which need to be added to V 0 in Eq.
(3) to result in an I(2+1)HDM potential which is only Z 2 -symmetric. We shall not consider CP-violation in this paper, therefore, we require all parameters of the potential to be real.

Mass eigenstates
We define the doublets as with two inert doublets (φ 1 and φ 2 ) and one active doublet (φ 3 ) where the latter plays the role of the SM Higgs doublet, with h being the SM-Higgs boson. The CP-even/odd neutral Z 2 -odd fields from the inert doublets could in principle be DM candidates since only the fields from the active doublet couple to the fermions. To stabilise the DM candidate from decaying into SM particles, we make use 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 FCNCs are avoided as the extra doublets are forbidden to decay to fermions by Z 2 conservation. Note that the Yukawa Lagrangian in this model is identical to the SM one, with φ 3 playing the role of the SM Higgs doublet: The point (0, 0, v √ 2 ) becomes the minimum of the potential at Expanding the potential around this vacuum point results in the mass spectrum below, where the pairs of inert scalar/pseudo-scalar/charged base fields (H 0 1,2 , A 0 1,2 , H ± 1,2 ) are rotated by: 4 into the mass eigenstates identified in boldface fonts.
There are two generations of physical inert states; fields from the first generation, (H 1 , A 1 , H ± 1 ) are chosen to be lighter than the respective fields from the second generation, (H 2 , A 2 , H ± 2 ), with H 1 being the lightest of them all, i.e., a DM candidate: The mass spectrum has the schematic form shown in Fig.(2), provided the CP-even neutral inert particles are lighter than the CP-odd and charged inert particles, which puts the following constraints on the parameters: We also consider cases where the mass alignment is changed, but where H 1 is always the lightest inert state. In the remainder of the paper the notations H 1 and DM will be used interchangeably. Figure 2: Schematic mass-squared spectrum of the Z 2 symmetric I(2+1)HDM, where Σ = 4µ 4 12

Theoretical constraints
Theoretical requirements of positivity of mass eigenstates, bounded-ness of the potential and positive-definite-ness of the Hessian put the following constraints on the potential.
1. Positivity of the mass eigenstates Bounded-ness of the potential For the V 0 part of the potential to have a stable vacuum (bounded from below) the following conditions are required 5 • λ 11 , λ 22 , λ 33 > 0 (16) • λ 12 + λ 12 > −2 λ 11 λ 22 We also require the parameters of the V Z 2 part to be smaller than the parameters of the V 0 part: 3. Positive-definite-ness of the Hessian For the point (0, 0, v √ 2 ) to be a minimum of the potential, the second order derivative matrix must have positive definite determinant. Therefore, the following constraints are required:

Experimental constraints
Relevant constraints limit the parameters from different experiments.

Collider constraints
• LEP limits Measurements done at LEP limit the invisible decays of Z and W ± gauge bosons, require that [32,33] • Also, LEP provides a model-independent lower limit for the mass of the charged scalars: Searches for charginos and neutralinos at LEP have been translated into limits of region of masses in the I(1+1)HDM [33] where for m H < 80 GeV and m A < 100 GeV the following region is excluded We have taken this limit into account in our numerical studies for any pair of CP-even and CP-odd particles. • where h → inv. represents the SM-Higgs decay to invisible particles channel which in our case is the h → H 1 H 1 channel. Global fits on Higgs signal strengths require the invisible Br of a Higgs boson with SM couplings but additional invisible decay modes to be limited to [36] • Br(h → inv.) < 23% at 95% CL The invisible Higgs decay into two scalar particles, CP-even or CP-odd, denoted by S is given by where m S < m h /2 is the scalar particle mass and λ is its coupling to the SM Higgs boson. In first approximation (with only one scalar, say H 1 , having a mass below m h /2 ≈ 62.5 GeV) the invisible decay rate is such that The limit from Eq. (24) leads to strong constraints on the H 1 H 1 h coupling (roughly λ 0.02 for masses m H 1 m h /2). In general, this will lead to tension between the LHC limits, which favour smaller couplings, and the relic density limits (discussed below), which favour larger couplings needed for effective DM annihilation.
Note that the presence of additional charged scalar states, H ± 1,2 , may modify the Higgs diphoton decay channel and lead to deviation from the SM value defined as 6 : where Γ(h) SM and Γ(h) I(2+1)HDM are the total decay widths of the Higgs boson in the SM and the I(2+1)HDM, respectively, while Γ(h → γγ) SM and Γ(h → γγ) I(2+1)HDM are the respective partial decay widths for the process h → γγ.
However, if the future combined value with reduced uncertainties shows significant deviation from µ γγ = 1, it will provide strong constraints for multi-scalar models.

Dark Matter constraints
• Relic density constraints DM relic density Ω DM h 2 is constrained by the combined Planck and WMAP results to be [3] : which leads to the 3σ bound: If a DM candidate fulfils this requirement, then it constitutes 100% of the DM in the Universe. However, a subdominant DM candidate is allowed if its relic density is smaller than 0.1118. Regions of the parameter space corresponding to Ω DM h 2 larger than the Planck upper limit are excluded.
In our work we use the micrOMEGAs 3.5 package to compute the relic density [38]. All annihilation and coannihilation channels are taken into account, including final states with one or two virtual gauge bosons in all cases relevant for the chosen values of masses.
• Direct detection constraints Neutral and non-relativistic WIMPs are expected to interact mainly with the atomic nuclei, whose nuclear recoil energy is to be measured by the DM detector. The current strongest upper limit on the Spin Independent (SI) scattering cross section σ DM −N is provided by the LUX experiment [39]: Limits from XENON100 (2012) are slightly weaker [40], with the strongest exclusion limit • Indirect detection constraints The indirect evidence for DM can be provided by measurements of the excess in the cosmic ray fluxes coming from the annihilation of DM in the Milky Way halo. The strongest constraints for light DM 7 annihilating into bb or τ + τ − is provided by the measurements of the gamma-ray flux from Dwarf Spheroidal Galaxies by the Fermi-LAT satellite, ruling out the canonical cross-section [41,42]: For the heavier DM candidates PAMELA and Fermi-LAT experiments provide similar limits of in the bb, τ + τ − or W + W − channels [43]. HESS measurements of signals coming from the Galactic Centre set limits of σv ≈ 10 −25 − 10 −24 cm 3 /s for DM masses up to TeV scales [44].

DM (co)annihilation in the I(2+1)HDM
The relic density of the scalar DM candidate, S, after freeze-out is given by the solution of the Boltzmann equation: where the thermally averaged effective (co)annihilation cross-section contains all relevant annihilation processes of any S i S j pair into SM particles: where Therefore, only processes for which the mass splitting between a state S i and the lightest Z 2 -odd particle S are comparable to the thermal bath temperature T provide a sizeable contribution to this sum. The I(2+1)HDM studied here shares many features of a Higgs-portal DM model. In a large region of parameter space the most important channel for the DM annihilation is shown in Fig.(3). The efficiency of this annihilation channel depends both on the value of DM mass and its coupling to the Higgs particle. In general, if m DM < m h /2, then one needs a coupling that is relatively large to produce relic density in agreement with Eq. (1). In this case a small DM-Higgs coupling results in too large a relic density and leads to the overclosure of the Universe. Note that a relic density below the Planck value does not exclude the DM candidate in the model, but requires another component to the DM to complete the deficit in the relic density. Figure 4: Diagrams contributing to the total annihilation cross section when m DM > m W , where V is any of SM gauge bosons.
The diagrams shown in Fig.(4) also contribute to the total annihilation cross section, where V is any of SM gauge bosons. Contribution from these diagrams is suppressed when the DM mass is smaller than m W , however, as studies have shown, diagrams with off-shell gauge bosons may be very important for m DM < m W in models such as the I(1+1)HDM. In our analysis the diagrams shown in Fig.(5) are also included. Coannihilation effects play an important role in scenarios with multiple particles that are close in mass. Particles up to 20% heavier than the DM candidate may influence the DM relic density. Therefore, the coannihilation diagrams should be included in calculating the effective annihilation cross section. The coannihilation channels shown in Fig.(6) appear in our studies.

Coannihilation scenarios
We introduce the following parameters for the mass splitting between the DM candidate H 1 and other inert particles: • between H 1 and the other CP-even state H 2 • between H 1 and the pseudoscalar/charged state from the (same) lighter generation • between H 1 and pseudoscalar/charged state from the (other) heavier generation The following scenarios are relevant for DM relic density studies: The DM particle is significantly lighter than other inert particles. This case is similar to the standard Higgs-portal approach, especially if m H 1 is very light. One should, however, take into account the importance of annihilation channels with virtual gauge bosons even in the m DM < m W region.
There is a small difference between the masses of H 1 and H 2 , thus the coannihilation between those two particles occurs through H 1 H 2 → h →f f while there is no coannihilation between scalar and pseudo-scalar particles.
There is a small difference between the masses of H 1 and A 1 while the CP-even particles mass splitting is large. The situation in this case is similar to the I(1+1)HDM with small CP-even/CP-odd mass splitting.
There is a small difference between the masses of H 1 , A 1 , H 2 and possibly A 2 (depending on the size of δ A and ∆), which can lead to coannihilation between all neutral inert particles.
The above scenarios are the only phenomenologically relevant cases in the m DM < m h /2 region. This is due to the LEP limits on the charged scalar excluding m ± H i < 70 − 90 GeV, which means that it cannot be close in mass with the DM candidate in this mass region.
However, for m DM > m W , coannihilation with H ± i is also allowed leading to the two following cases, which are analogous scenarios to C and D with δ A ( ) ↔ δ C ( ): These two scenarios can be realized in the m DM > m W region which is dominated by very effective annihilation into gauge bosons. Usually one needs strong cancellation effects to obtain proper relic density.
A final possible coannihilation scenario is the following.
All inert particles have similar masses, which resembles the situation in the heavy mass region in the I(1+1)HDM.
In principle, a scalar DM candidate with acceptable relic density abundance can have a mass range from a few GeV to a few TeV. A light DM candidate usually requires a relatively large DM-Higgs coupling. As the mass grows, the annihilation into Higgs bosons becomes more effective, especially around the Higgs resonance (m DM ≈ m h /2) and a smaller DM-Higgs coupling is required. For m DM > m W annihilation into gauge bosons is very effective which usually leads to relic density values well below the observed limit. The cancellation effects between diagrams in Fig.(4) and coannihilation effects can help restore the relic density value.
In section 5 we will numerically investigate all possible scenarios for DM (co)annihilation in the I(2+1)HDM in different regions of DM mass. Here, we briefly comment on the gross feature of each benchmark configuration.
For a start, note that, in the I(2+1)HDM, DM-nucleon scattering is through the exchange of the Higgs particle in the t-channel. This cross section depends on the Higgs-DM coupling, g H 1 H 1 h , and, as discussed, in many cases in the Higgs-portal models the coupling needed for the proper relic density is too large to reconcile with results from direct detection experiments. As it will be shown in section 5, it is often the case in scenario A.
Furthermore, note that the coannihilation effects may help in satisfying these constraints. In scenarios B-G an acceptable value of relic density is obtained through coannihilation processes, while the scattering still depends only on the t-channel Higgs-exchange. DM-Higgs coupling in such cases may be smaller than in case A, and so the scattering cross section may lie below the current experimental limits.
Finally, note that the direct detection experiments set strong limits on the scattering through Z-exchange in regions where scalar and pseudoscalar states are (nearly) degenerate. A non-zero mass splitting between the lightest scalar and the lightest pseudoscalar particle needs to be larger than the kinetic energy of DM in our galactic halo, so that scattering through Z-exchange is forbidden kinematically. This sets the lower limit on δ A to be a few 100 keV [45]. In principle, this limit has important consequences for NHDMs, where degenerated states appear in a natural way. In the I(2+1)HDM though, the degeneracy between H i and A i is lifted by non-zero λ 2 and λ 3 .

Simplified couplings in the I(2+1)HDM
We study the simplified case of the I(2+1)HDM where all the parameters related to the first inert doublet are k times the parameters related to the second doublet: However, we assume no specific relation among the other parameters of the potential. The masses of the CP-even neutral states in this case are simplified to with the mixing angle between the CP-even states given by and m 2 A 2 , m 2 A 1 have similar values with Λ φ 2 and θ h replaced by Λ φ 2 , Λ φ 2 and θ c , θ a , respectively.
Note that the positivity of mass eigenstates puts the following limits on the acceptable values of k: We will study several cases in the simplified I(2+1)DM which are listed in Tab. (1). First, cases of k = 0 (section 5.2) and k = 1 with µ 2 12 = 0 (section 5.3) are discussed briefly for completeness, but they do not provide any solution to the problems of the Higgs-portal DM scenario. We then study the case of k = 1 with µ 2 12 = 0 in detail since it represents all features of this model clearly. Our numerical DM studies are done mostly for the selected benchmark points which exhibit typical behaviour of a particular scenario. In each case we present the resulting relic density plots and in section 5.5 the different cases are compared to each other. The other cases from Tab.(1) are discussed briefly since they do not present any new features of the model. k mixing parameters section , g H1H1h , µ 2 12 , k 5.7 Table 1: The cases with different values of k and mixing between the inert doublets studied here alongside the paper section where they are dealt with.

The k = 0 case
With k = 0 the model is reduced to the two doublet case, the I(1+1)HDM, which we briefly review here. The Z 2 -symmetric I(1+1)HDM potential of two doublets, one active and one inert, is commonly written as: with all real parameters and g Z 2 = diag(1, −1). This symmetry is respected by the vacuum alignment (v, 0) and the neutral fields from the inert doublet, φ 2 , are viable DM candidates.
In this model there are three distinctive regions of m H where one can expect to obtain proper relic density [16,17,32,46,47,48,49,50,51,52,53].  The I(1+1)HDM parameter space is strongly constrained by the recent LHC and direct detection results which lead to exclusion of the low DM mass region for m DM 55 GeV [54,55,56]. This is due to the incompatibility between the relic density limits, which require the Higgs-DM coupling to be relatively big to ensure the efficient-enough DM annihilation, and the LHC µ γγ and Br(h → inv.) constraints, which prefer much smaller values of such a coupling. The medium mass region for m DM < m W , is instead in agreement with the current experimental results, however, an enhancement in the diphoton Higgs decay channel is disfavored [54]. Finally, LHC and direct detection limits do not provide significant constraints on the heavy mass region, however, if a significant enhancement in µ γγ is confirmed, then it is not possible to reconcile it with relic density limits.

The k =1 case with vanishing mixing
In this case the two inert doublets are completely degenerate: Furthermore, the µ 2 12 = 0 leaves only 4 base parameters, and the 5 parameters describing the self-interactions of inert particles, i.e., λ 2 , λ 11 , λ 22 , λ 12 , λ 12 , which are not relevant for the relic density analysis. Note that the potential becomes Z 2 × Z 2 symmetric after imposing these equalities, as The inert particle mass spectrum in this case has the following form: The parameters of the model in terms of the physical parameters are: The relevant Feynman rules are: In first approximation, i.e., by neglecting the self-interactions of the inert doublets (for example, H 1 H 1 ↔ H 2 H 2 ), one can treat this case as a doubled I(1+1)HDM, with two DM candidates, H 1 and H 2 . They have degenerated masses and identical interactions, as noted above, and they contribute equally to the DM relic density Ω DM h 2 : where Ω H h 2 is the relic density of a single DM candidate from the I(1+1)HDM.
To fulfil the Planck limit, Ω H h 2 needs to lie between 0.0559 and 0.064, meaning that the DM annihilation should be more effective than in the I(1+1)HDM. This requires even bigger values of the DM-Higgs coupling, which, for m DM < m h /2, would lead to even larger values of Br(h → inv.), making it even harder to satisfy relic density and LHC bounds simultaneously, unless it is in the Higgs-resonance region.

The k = 1 case with non-vanishing mixing
Similar to the previous case, the two inert doublets are perfect copies of each other, with equal parameters. The non-zero CP-even mixing angle from Eq. (43) is The CP-even mass spectrum is therefore of the following form: Assuming µ 2 12 > 0 and θ h belonging to the 1 st quadrant, one has which indeed makes H 1 the lightest among the inert particles and therefore our DM candidate. Further assuming that all θ i 's are in the 1 st quadrant, the mass spectrum has the following form: The base parameters can then be expressed in terms of Finally, the following equations relate different parameters: , is the coefficient of the H 1 H 1 h term in the potential. It is interesting to note that the equations for λ 23 and λ 2 are identical to the corresponding relation for λ 4 and λ 5 in the I(1+1)HDM case and that the phenomenology of the model depends on masses of inert particles and one coupling only.
The relevant Feynman rules are: In the following subsections, we study the k = 1 with µ 2 12 = 0 case in detail for the scenarios proposed in section 4. In the second part of this section, scenarios with DM mass from the m W > m H 1 > m h /2 range are discussed where the effects of DM annihilation into gauge bosons play an important role and lead to a rather different phenomenology.

Open invisible channels (m DM < m h /2)
• Case A with ∆ = 50 GeV We choose the following values of masses as an example for our numerical studies: Mass splittings δ A , δ C are of the order of 50 GeV, with H 1 being much lighter than other inert particles. The resulting relic density is plotted versus the DM-Higgs coupling in Fig.( annihilate through gauge bosons. However, the existence of these channels is massand case-dependent (see also Fig.(12)).
Note that, if the mass splittings δ A,C , δ A,C are large enough to forbid any coannihilation effects, the lightest Z 2 -odd particle is effectively decoupled from the Z 2 -odd sector and the exact values of masses do not influence the DM phenomenology. In that sense DM studies, i.e., relic density measurements and direct detection experiments, do not put any additional constraints on heavier Z 2 -odd particles and input from particle physics is needed. In this scenario, with ∆ relatively small, one expects the H 2 H 1 → h coannihilation effects to show in the relic density plots. Note that the result is sensitive to the value of ∆. Here we show the relic density plot for the ∆ = 8 GeV case (referred to as case B 8 ) and comment on the ∆ = 1 GeV case (referred to as case B 1 ).
One should note first that for the discussed k = 1 case there is no tree-level H 2 H 1 h coupling and so H 1 H 2 → h → ff diagrams do not exist. Therefore, in comparison with case A the annihilation of DM is not affected, which is shown in the the first considered setup, ∆ = 8 GeV, presented in Fig.(9a). Note that Fig.(9a) is almost identical to case A in Fig.(8a). The reason being that the ∆ mass splitting is large enough and so the H 2 H 2 → h → ff diagrams do not interfere with the thermal evolution of DM relic density.
Case B with ∆ = 1 GeV (case B 1 ) With smaller mass splitting, ∆ ≈ 1 GeV, the relic density evolution could be affected. In this case the second CP-even particle H 2 acts almost like the second DM candidate discussed in the k = 1 no-mixing scenario in section 5.3. The obtained relic density is larger than in case A and larger couplings are needed to fulfil the Planck bounds. This is even harder to reconcile with limits from Br(h → inv.), especially since now there are two invisible channels open. Therefore, we conclude that scenario B cannot provide any solution to the problems of Higgs-portal DM with m DM < m h /2. • Case C with δ A = 8 GeV and ∆ = 50 GeV In this case H 1 and A 1 are very close in mass, while other inert particles are heavy in comparison: As a result, there is coannihilation between H 1 and A 1 whose effects are visible in Fig.(10a) compared to the previous cases.
For this setup coannihilation effects lead to an enhanced cross section. Coannihilation becomes so effective that, even for small couplings, relic density does not reach the observed value. In fact, coannihilation processes are so strong that for every value of DM mass (below m h /2) relic density is below current Planck/WMAP limits. This situation does not result in the exclusion of this parameter space but rather corresponds to a subdominant DM candidate.
Note that the δ A chosen here is the boundary value of mass splitting between scalar and pseudoscalar in the I(1+1)HDM, obtained by translation of the null-searches for charginos and neutralinos at LEP-II [33]. Increasing this value to δ A = 10 GeV and thus reducing the strength of coannihilation -while still keeping it possible -allows for Ω DM h 2 within Planck limits for relatively small values of DM-Higgs coupling.
(a) • Case D with δ A = 7 GeV and ∆ = 1 GeV In this case the masses of all neutral inert particles H 2,1 , A 2,1 are relatively close. Two important coannihilation effects taking place here are the following. Firstly, the H 1 A 1 → Z coannihilation, which leads to a decrease in the relic density (similar to most SUSY models). As discussed in case C above, usually it leads to Ω DM h 2 below the observed value. Secondly, the H 2 H 2 → h (co)annihilation effect, which leads to an increase in the relic density (similar to the case in UED theories) by affecting the DM production rate, as presented in case B 1 with small ∆. These two effects combined will allow for small g H 1 H 1 h values and result in sufficient relic density, which is desirable since smaller g H 1 H 1 h leads to less stringent bounds from direct detection experiments.
The benchmark points studied above represent a typical behaviour in this region of DM mass. Clearly, as we have shown, the mass splitting between a DM candidate and other Z 2 -odd particles influences the freeze-out mechanism and the final value of DM relic density. Fig.(12) shows the allowed g H 1 H 1 h coupling in different mass regions, where the grey area inside the red (case A) and blue (case D) curves are excluded by relic density data. The white region outside the curves represents smaller relic density abundance than the observed (a) value. It is easy to see that, apart from the Higgs resonance region, couplings that lead to the proper value of relic density are much smaller in case D than they are in case A for the same values of masses. Direct detection limits for m H 1 < m h /2 Fig.(13) compares the direct detection limits in cases A and D, where the points above the LUX limit (black line) are excluded. The vertical line at 62.5 GeV represents the Higgs resonance mass region. For masses below m h /2 the direct detection limits constrain case A much more severely than they limit case D. Masses below 53 GeV in case A, corresponding to large g H 1 H 1 h (see Fig.(12)), are completely excluded in case A and only points around the Higgs-resonance region (denoted by the black vertical line) are allowed. In case D, however, almost the whole mass region (in the m H 1 < m h /2 range) surviving the relic density constraints agrees with the direct detection limits. Around the Higgs-resonance region very small g H 1 H 1 h ≈ 10 −3 couplings lead to the Ω DM h 2 in agreement with Planck. This corresponds to a very small DM-nucleon cross section, σ DM −N ≈ 10 −50 cm 2 and makes the direct detection experiments impractical due to coherent neutrino scattering.
Similar to other Higgs-portal DM models, in the small g H 1 H 1 h region, loop effect contributions are in principle important and one has to move beyond tree-level calculations for more detailed descriptions, both in the relic density estimates and the scattering cross section analysis.

Higgs invisible decays
Constraints arising from limits on Higgs invisible decays can easily be estimated with making some assumptions; Firstly, the Higgs-decay channels in the I(2+1)HDM are similar to the Case D Excluded Figure 12: Relic density constraints on the mass of the DM candidate and its coupling to SM Higgs boson, with the white and gray regions representing too little and too much relic abundance, respectively. The red band is where sufficient amount of relic abundance is produced in case A (and identically in case B 8 ), and the blue region shows the accepted window in case D. In case C the relic density is always below the observed value. ones in the SM (in particular, contributions to the h → γγ are negligible). Secondly, the total decay width in the I(2+1)HDM is changed with respect to the SM only through the invisible decays. Under these assumptions the Higgs invisible Br is given by the ratio of the decay widths: The sum runs over particles S i of masses smaller than m h /2, meaning that in the I(2+1)HDM there can be up to four particles contributing to Br(h → inv.) 8 . If only one particle (H 1 ) is lighter than m h /2 (Case A, Fig.(14a)  It is interesting to compare these limits, firstly, with the regions of proper relic density and, secondly, with exclusion limits obtained from the direct detection experiments. From  Fig.(14a) it is clear that, for case A, it is not possible to fulfil Planck measurements and LHC measurements for masses m H 1 53 GeV, as the region in agreement with LHC corresponds to having too much DM in the Universe. These limits are comparable with those provided by the LUX experiment and delineating a region which fulfils both requirements is related to entering the Higgs resonance region, where a very small coupling still results in the good relic density, without violating LHC or direct detection bounds. The situation is different in case D, see Fig.(14b), which in general is not constrained by direct detection limits. Here we can see that the LHC results provide severe constraints; while it is possible to fulfil Br(h → inv) < 37% for masses below m h /2 (such limits again are comparable to those provided by LUX), using the global fit value for Br(h → inv) excludes DM candidates with masses m H 1 53 GeV, just like in case A.
To summarize, cases A and D depend on different mechanisms to obtain proper relic density and therefore are differently constrained by direct detection experiments. However, LHC limits constrain them equally, leaving only m H 1 53 GeV.

Closed invisible channels (m W > m DM > m h /2)
In this mass region the DM phenomenology is heavily influenced by the annihilation into gauge bosons, which leads to a different behaviour compare to the m H 1 < m h /2 region .
The relic density values are dominated by three couplings, g DM V V , g hV V and g H 1 H 1 h . The first two couplings are set by gauge interactions, therefore, the behaviour of the relic density plots are ruled by the value of g H 1 H 1 h . Since this type of annihilation is given by the strength of gauge couplings and therefore is usually very effective, proper relic density is obtained due to the cancellation effects between H 1 H 1 → V V and H 1 H 1 → h → V V channels. For g H 1 H 1 h > 0 the annihilation cross section is large leading to small relic density values. As g H 1 H 1 h goes towards more negative values the annihilation cross section reduces, leading to larger values in relic density.
Below we present in detail the numerical results obtained for case A (∆, δ i ≈ 50 GeV) and case D (∆ = 1 GeV, δ A = 7 GeV) studied in the previous section.
• Case A Fig.(15 • Case D Results for case D are presented in Fig.(16). The existence of coannihilation channels drives the relic density to smaller values in comparison to case A in Fig.(15). However, the difference is not nearly as pronounced as it was in the m H 1 < m h /2 region.   (a) The white region outside the curves represents a relic density abundance smaller than the observed value. It is clear that the two scenarios correspond to very similar values of coupling and the mass splitting does not play as important a role here as previously. Therefore, the direct detection exclusions will be similar in both cases (as shown in Fig.(13)), with case D being slightly less constrained than case A. Let us finally comment briefly on the other two scenarios discussed in the previous section, namely cases B and C. As before, for case B 8 with ∆ = 8 GeV, we reproduce results from  Figure 17: Relic density constraints on the mass of the DM candidate and its coupling to SM Higgs boson, with the white and gray regions representing too little and too much relic abundance respectively. Note that the regions are overlapping.
case A. If ∆ = 1 GeV, the couplings which correspond to the proper relic density will be 30 − 50% larger than for the same mass in case A. This will lead to a larger DM-nucleon scattering cross section and stronger exclusion from direct detection experiments. Case C again corresponds to a subdominant DM candidate with relic density below the observed value. Detailed plots showing differences between discussed cases are presented in section 5.5.

Heavy DM mass (m DM m W )
As the mass splitting between inert particles is given by the quartic couplings λ 2 and λ 3 , unitarity bounds for these parameters will lead to an almost degenerated mass spectrum in the heavy mass regime. Therefore, the only scenario leading to acceptable relic density values is case F, discussed in section 4.1, where all inert particles have similar masses. Proper relic density here is obtained through cancellations among diagrams (see Figs.(4),(6), (7)) and a relatively large value of g H 1 H 1 h is also needed. Similar to other multi-scalar models, this mass region escapes both the current direct detection limits and LHC constraints. However, interesting indirect detection signatures (connected to the possibility of internal bremsstrahlung H 1 H 1 → W + W − γ processes generated through the exchange of any of the two charged scalars H ± 1,2 in the t-channel) could arise here, which will require further studies.

Summary of k = 1 results for fixed DM mass
Below we present the detailed comparison between different scenarios studied in section 5.4 for k = 1 and fixed values of DM mass. Figs. (18) and (19) show the relic density of the DM candidate vs. DM-Higgs coupling for DM masses less than and greater than half the Higgs mass, respectively.
In all plots case A (green) and case B 8 with ∆ = 8 GeV (black) are indistinguishable; coannihilation effects with H 2 play no role here. To show the relevance of the value of ∆ we also plot case B 1 with ∆ = 1 GeV (red). It is clear that in case B 1 the coupling which gives Ω DM h 2 in agreement with Planck measurements is equal or larger than the one from case A. Therefore exclusion limits in case B 1 are stronger and this scenario does not provide any solution to the problems of Higgs-portal DM models.
Case C (blue), which present an equivalent scenario to that of coannhilation in the I(1+1)HDM, is always below the Planck limit. Case D (purple) generally corresponds to couplings smaller than in cases A and B.
As a function of the DM mass, we observe the following: • • With m H 1 growing towards m W , the resulting relic density decreases.

Conclusions
We have studied a model with Two Inert Doublets plus One Higgs Doublet (I(2+1)HDM). The two inert doublets are Z 2 -odd and the active doublet is Z 2 -even which plays the role of the SM Higgs doublet. The I(2+1)HDM may be regarded as an extension to the I(1+1)HDM (also known as the Inert Doublet Model (IDM)). The I(2+1)HDM contains 4 neutral and 4 charged inert particles, which is double the particle content of the I(1+1)HDM. The lightest particle amongst the inert ones, stabilised by an exact Z 2 symmetry, provides a viable DM candidate. We have then studied the DM phenomenology in different parameter scenarios of this model. Similar to I(1+1)HDM, the I(2+1)HDM contains all features of a Higgs-portal scalar DM model, with certain modifications. In such models, the Higgs-DM coupling, g H 1 H 1 h , governs the DM annihilation rate σv , the DM-nucleon scattering cross section σ DM −N and the Higgs invisible decays. This is the scenario studied in case A, which is the closest to the pure Higgs-portal DM model. Imposing current Planck, LUX and LHC constraints we found that, as in the I(1+1)HDM, it is not possible to have a DM candidate with mass below 53 GeV.
The tension in simultaneously fulfilling current experimental constraints on the processes described above can be lifted by introducing coannihilation processes offered by the rich inert scalar sector. In such case the relation between Ω DM h 2 and σ DM −N is destroyed, which is the scenario studied in case D. These coannihilation processes allow for sufficient relic density values in different DM mass regions compared to the I(1+1)HDM, where, for m DM < m W , there is only a possibility of destructive coannihilation with a pseudoscalar inert particle. Constructive coannhilation with the remaining charged scalar is not possible in this region, as sufficiently light charged scalars are excluded by LEP measurements. We found that in the I(2+1)HDM coannihilation between all neutral particles, denoted as case D, can lead to a proper relic density for masses below m W . Furthermore, this region, especially for masses m H 1 50 GeV, is in agreement with the current direct detection limits provided by the LUX experiment.
This new region of masses can survive also the current direct limits on Higgs invisible decays provided by the ATLAS and CMS collaborations at the LHC, i.e. 37%. However, we found that it is not possible to be in agreement with the global fit value of Br(h → inv) < 20%, unless the DM mass is above 53 GeV. For this model then, the LHC limits are stronger than those provided by direct detection experiments.
In heavier DM mass regions (m W > m DM > m h /2), where annihilation through gauge bosons contributions is dominant and reduces the relic density values, the g H 1 H 1 h coupling is driven towards negative values. The coannihilation effects are not as visible in the m DM > m h /2 region since the contribution from gauge boson annihilation plays an important role in this region. This contribution is more pronounced closer to m W threshold, where case D corresponds to σ DM −N of an order of magnitude smaller than case A.
All studied cases, as the general Higgs-portal models, are in agreement with the current experimental constraints in the Higgs resonance region (∼ m h /2). Here, very small g H 1 H 1 h couplings are required to fulfil relic density constraints. Such small g H 1 H 1 h values lead to very small DM-nucleon cross section, which makes the direct detection experiments impractical in this region, especially since it reaches the boundary of coherent neutrino scattering which has then to be taken into account. Also, in such small g H 1 H 1 h region, loop-effect contributions are important and one has to move beyond tree-level calculations for more detailed descriptions, which is the case in all Higgs portal models.
It is also interesting to note that the strong limits on the scattering through Z-exchange by direct detection experiments do not apply to this model since a natural mass splitting between the scalar and pseudoscalar particles is provided by non-zero λ 2 and λ 3 parameters in the potential.
We would also like to comment here on the special case C studied in section (5.4), which always results in relic density values below the observed amount and therefore cannot account for the whole DM in the Universe. It could however provide a subdominant DM candidate in a multi-component DM model. One of the examples of such model is the I(4+2)HDM [57] in which the scalar sector consists of two copies of the scalar sector of our I(2+1)HDM. Case C could then become a viable scenario with sufficient relic density values.

Note Added
As this paper was being finalised a related paper appeared which also analyses DM in a model with 'Two Inert Doublets plus One Higgs Doublet' [58].