Two component dark matter with inert Higgs doublet: neutrino mass, high scale validity and collider searches

The idea of this work is to investigate the constraints on the dark matter (DM) allowed parameter space from high scale validity (absolute stability of Higgs vacuum and perturbativity) in presence of multi particle dark sector and heavy right handed neutrinos to address correct neutrino mass. We illustrate a simple two component DM model, consisting of one inert SU(2)L scalar doublet and a scalar singlet, both stabilised by additional Z2×Z2′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal{Z}}_2\times {\mathcal{Z}}_2^{\prime } $$\end{document} symmetry, which also aid to vacuum stability. We demonstrate DM-DM interaction helps achieving a large allowed parameter space for both the DM components by evading direct search bound. High scale validity puts further constraints on the model, for example, on the mass splitting between the charged and neutral component of inert doublet, which has important implication to its leptonic signature(s) at the Large Hadron Collider (LHC).


Introduction
Discovery of the 'Higgs' boson at Large Hadron Collider (LHC) in 2012 [1,2] strongly validates the Standard Model (SM) of particle physics as the fundamental governing theory of Strong, Weak and Electromagnetic interactions. However many unresolved issues persist. For example, the electroweak (EW) vacuum turns out to be metastable [3][4][5][6][7][8][9] with the present measured value of Higgs mass (m h ∼ 125.09 GeV [10]) and top quark mass (m t ∼ 173.1 GeV [10]). Large uncertainty in the measured value of m t can even make EW vacuum JHEP03(2020)090 neutrinos to yield correct light neutrino mass through additional Yukawa coupling with SM Higgs. This in turn may alter the Higgs vacuum stability condition at high scale (larger than the RH neutrino masses) with large neutrino Yukawa coupling [56][57][58][59][60][61][62][63][64][65][66][67][68]. Inclusion of an inert Higgs doublet along with right handed neutrinos can also generate the light neutrino mass radiatively [69]. In that case, the usual neutrino Yukawa coupling involving the SM Higgs doublet can be forbidden with the choice of appropriate discrete symmetry, while a new Yukawa interaction involving the inert Higgs doublet is present. However in such a scenario (with radiative light neutrino mass generation), the EW vacuum stability would not be much different from that of the SM as the running of the Higgs quartic coupling remains unaffected by the presence of this new Yukawa interaction.
Our model under scrutiny addresses some of the interesting features in view of above discussion. We address a multipartite dark sector consisting of a SU(2) L doublet scalar (the IDM) and a scalar singlet, both stabilized by additional Z 2 × Z 2 symmetry (for an earlier effort, see [41]) and provide a two component DM set up. The presence of DM-DM interactions enlarge the available parameter space significantly (by reviving the below 500 GeV mass region for both the DM candidates), while the inert DM can also produce leptonic collider signature at LHC. We augment the model with heavy RH neutrinos to address the light neutrino masses. Though the set-up with inert Higgs doublet and RH neutrinos opens up both the possibilities of generating light neutrinos masses, either through type-I seesaw or via radiative generation, we consider here the neutrino mass generation through type-I seesaw and accordingly choose the appropriate discrete charges of the fields so as to forbid the Yukawa interaction involving SM lepton doublet, RH neutrinos and inert Higgs doublet. The primary reason behind such a choice is to involve the study of the stability of the EW vacuum with large neutrino Yukawa couplings (involving SM Higgs doublet) as mentioned before. Then the parameter space obtained from DM phenomenology will be further constrained from vacuum stability criteria. In brief, we want to study a multipartite DM scenario which would be adversely affected by the nonzero light neutrino mass in terms of EW vacuum stability. Note that while the neutrino Yukawa coupling involving RH neutrinos and SM Higgs tends to destabilze the EW vacuum, the presence of the additional scalars in the set-up tends to stabilize it [56][57][58][59][60][61][62][63][64][65][66][67][68]. So, we take up an interesting exercise of validating the model from DM constraints, neutrino masses and high scale validity (absolute stability of the Higgs vacuum and perturbativity). This analysis provides some important conclusions, which are phenomenologically viable at LHC.
Let us finally discuss the plan of the paper. In section 2, we discuss the model construct in details. Section 3 presents possible theoretical and experimental constraints on the model parameters. Then in sections 4, 5 and 6 subsequently, we discuss DM phenomenology. Indirect search constraints are discussed in 7. In section 8, we investigate the high scale validity of the model. Section 9 summarises collider signature(s) in context of the proposed set up. Finally we conclude in section 10. Tree level unitarity condition is elaborated in appendix A. The high scale stability condition on the single component DM frameworks in presence of RH Neutrino is chalked out in appendix B for the sake of comparison. All the constraints together on the model parameter space along with different choices of RH neutrino mass and Yukawa couplings are listed in two tables in appendix C.

JHEP03(2020)090
BSM and SM Higgs Fields Table 1. Charge assignments of the BSM fields assumed in the model under G as well as that of SM Higgs. The U(1) Y hypercharge is chosen as Q = T 3 + Y /2.

The model
The model is intended to capture the phenomenology of two already established DM frameworks involving that of a singlet scalar and that of an inert scalar doublet together with right handed neutrinos to address neutrino mass under the same umbrella. Therefore, we extend SM by an inert doublet scalar (Φ) and a real scalar singlet (φ) and include three RH Majorona neutrinos N i (i = 1, 2, 3) in the set up. The lightest neutral scalar mode of the IDM and φ are the DM candidates provided an appropriate symmetry in addition to that of SM stabilizes both of them. This is minimally 1 possible by introducing an additional Z 2 ×Z 2 discrete symmetry under which all SM fields along with the right handed neutrinos transform trivially and the other additional fields transform non-trivially as tabulated in the table 1. We also note the charges of SM Higgs (H) explicitly in table 1, as it will be required to form the scalar potential of the model. Note here, that charges of the two DM candidates (Φ and φ) are complementary, i.e. odd under either Z 2 or Z 2 for their stability. We also point out that the U(1) Y hypercharge assignment of Φ is identical to SM doublet H. Therefore the only SU(2) L × U(1) Y invariant terms are H † H, Φ † Φ, H † Φ and its conjugate.
The scalar Lagrangian reads as: where and g 2 , g 1 denote SU(2) L and U(1) Y coupling respectively.
The most relevant renormalizable scalar potential in this case is given by,

2)
JHEP03(2020)090 where, The Lagrangian involving right handed neutrinos can be written as, whereH = iσ 2 H * . We have considered three generations of RH neutrinos with {i, j} = 1, 2, 3, which can acquire Majorana masses and can possess Yukawa interactions with SM lepton doublet l L . Note here, that the charge assignment of the N fields then aid us to obtain neutrino masses through standard Seesaw-I [54,55] mechanism (as detailed later), while it also prohibits the operator likel L ΦN due to Z 2 charge assignment, and hence discards the possibility of generating the light neutrino mass radiatively. The ingredients and interactions of the model set up is described in the cartoon as in figure 1.
After spontaneous symmetry breaking, SM Higgs doublet acquires non-zero vacuum expectation value (VEV) as H = (0 v+h √ 2 ) T with v=246 GeV. Also note that neither of the added scalars acquire VEV to preserve Z 2 ×Z 2 and act as DM components. After minimizing the potential V (H, Φ, φ) along different field directions, one can obtain the JHEP03(2020)090 following relations between the physical masses and the couplings involved: where λ L = 1 2 (λ 1 + λ 2 + λ 3 ) and m h , m H 0 , m A 0 are the mass eigenvalues of SM-like neutral scalar found at LHC (m h = 125.09 GeV), heavy or light additional neutral scalar and the CP-odd neutral scalar respectively. m H ± denotes the mass of charged scalar eigenstate(s). The mass for φ DM will be rescaled as m 2 The independent parameters of the model, those are used to evaluate the DM, neutrino mass constraints are as follows: (2.8)

Theoretical and experimental constraints
We would like to address possible theoretical and experimental constraints on model parameters here.
Stability. In oder to get the potential bounded from below, the quartic couplings of the potential V (H, Φ, φ) must have to satisfy following co-positivity conditions (CPC) as given by [70,71], where CPC{i} denotes i th copositivity condition and µ is the running scale. The above conditions show that the model offers to choose even negative λ 1,2,3,φh satisfying the above conditions. However, as demonstrated in eq. (2.8), we use λ L and physical masses to be the parameters. Therefore, if we choose a specific mass hierarchy as: m H ± ≥ m A 0 ≥ m H 0 with positive λ L , we are actually using λ 1 to be positive while λ 2,3 negative abiding by the above conditions (see eq. (2.7)). The conditions CPC{i} as in eq. (3.1), will be used later in demonstrating stability of the potential in section 8.
Perturbativity. In oder to maintain perturbativity, the quartic couplings of the scalar potential V (H, Φ, φ), gauge couplings (g i=1,2,3 ) and neutrino Yukawa coupling Y ν should obey: and  Tree level unitarity. Next we turn to the constraints imposed by tree level unitarity of the theory, coming from all possible 2 → 2 scattering amplitudes as detailed in appendix A follows as [72,73]: and |x 1,2,3 | < 16π , (3.5) where x 1,2,3 be the roots of the cubic equation as detailed in appendix A.
Electroweak precision parameters. There exists an additional SU(2) L doublet (Φ) in our model in addition to a gauge singlet scalar (φ). As the vev of Φ is zero, it does not alter the SM predictions of electroweak ρ parameter [74]. However IDM, being an SU(2) L doublet makes a decent contribution to S and T parameters [75,76] which we will identify as ∆S and ∆T . The experimental bound from the global electroweak fit results on ∆S and ∆T using ∆U = 0 are given by: (3.6) at 1σ level with correlation coefficient 0.91 [77]. We show the constraint from ∆S and ∆T on the model parameter space in figure 2 using the standard formula as presented in [75,76]. In left plot, we scan 1σ fluctuation on ∆S in m A 0 −m H 0 versus m H ± −m H 0 plane for two different values of m H 0 = {80, 500} GeV. We see that for smaller m H 0 , the constraint is larger. In right panel, we show 1σ and JHEP03(2020)090 2σ limits from ∆T in m A 0 − m H 0 versus m H ± − m H 0 plane for a range of IDM mass m H 0 = {80 − 500} GeV. We can clearly see, that ∆T constrains the mass splitting much more than ∆S.
Higgs invisible decay. Whenever the DM particles are lighter than half of the SM Higgs mass, the Higgs can decay to DM and therefore it will contribute to Higgs invisible decay. Therefore, in such circumstances, we have to employ the bound on the invisible decay width of the 125.09 GeV Higgs as [10]: where and Γ(h → SM) = 4.2 MeV [10]. In this analysis, we have mostly focused in the region where m φ , m H 0 > m h /2, actually larger than W mass, i.e. m φ , m H 0 ≥ m W , so that the above constraint is not applicable.
Collider search constraints. Experimental searches for additional charged scalars and pseudoscalar in LEP and LHC provide bound on IDM mass parameters and coupling coefficients of IDM with SM particles.
(i) Bounds from LEP. The observed decay widths of Z and W bosons from LEP data restrict the decay of gauge bosons to the additional scalars and therefore provide a bound on IDM mass parameters as In addition, neutralino searches at LEP-II, provides a lower limit on the pseudoscalar Higgs (m A 0 ) to 100 GeV when m H 0 < m A 0 [78]. The chargino search at LEP-II limits indicate a bound on the charged Higgs to m H ± > 70 GeV [79].
(ii) Bounds from LHC. Due to the presence of SM Higgs and IDM interaction, charged scalars H ± take part into the decay of SM Higgs to diphoton. Thus it contributes to Higgs to diphoton signal strength µ γγ which is defined as [80][81][82][83] Now when IDM particles are heavier than m h /2, one can further write The analytic expression of Γ(h → γγ) IDM can be obtained as [80][81][82][83]:  where A SM represents pure SM contribution (see [80][81][82][83]). And Therefore it turns out that IDM contribution to µ γγ is function of both mass of the charged Higgs (m H ± ) and the coefficient of the trilinear coupling hH + H − i.e. λ 1 .
The measured value of µ γγ are given by µ γγ = 1.17 ± 0.27 from ATLAS [84] and µ γγ = 1.14 +0.26 −0.23 from CMS [85]. In figure 3, we show the variation of µ γγ as function of m H ± for different values of λ 1 . We also present the experimental limits on µ γγ from ATLAS in figure 3. Excepting for the resonance at m h/2 , we see that our choice of λ 1 is consistent with experimental bound. The larger is m H ± GeV, µ γγ → 1 i.e. approches to SM value. Also, we see that λ 1 > 0 diminishes µ γγ , while λ 1 < 0 tends to enhance it. In this analysis, we mostly consider m H ± > m h/2 and positive λ 1 within correct experimental limit.
Relic density of DM. The PLANCK experiment [14] provides the observed amount of relic abundance Furthermore, strong constraints exist from direct DM search experiments. In our analysis we will consider the most recent bound on direct detection cross section provided by XENON 1T [26]. Relic density and direct search allowed parameter space of the model will be evaluated in details. neutrino masses m ν i ≤ 0.12 eV as provided by PLANCK data [14,86] is incorporated. The present values of neutrino mass hierarchies and mixing angle can be found in [87,88]. Due to the presence of RH neutrino in the set up, the constraint from lepton flavor violating decay (LFV) (dominantly from µ → eγ) will be applicable [89][90][91]   The allowed parameter space of φ is depicted in left hand side (l.h.s.) of figure 5 in m φ − λ φh plane by the red dots. Direct search allowed parameter space from XENON1T data [26,27] using spin-independent DM-nucleon scattering cross section is shown by the blue dots. We therefore see that the model can only survive either in the Higgs resonance region (∼ m h /2) or at a very heavy mass 900 GeV. The under abundant (shown in yellow) and over abundant regions are also indicated.
IDM (H 0 ) as a single component DM have annihilation and co-annihilation channels for freeze-out due to both gauge and Higgs portal interactions as shown by the Feynman graphs in figures 6 and 7. The parameters, which govern the IDM phenomenology are [94]: Here we see again that allowed region from relic density and direct search constraint for IDM lies either in the small mass region m H 0 m W ∼ 80 GeV or in the heavy mass region m H 0 550 GeV. It is a well known result coming essentially due to too much annihilation and coannihilation of IDM to SM through gauge interactions [94]. The disallowed region 80 GeV < m H 0 < 550 GeV is often called desert region, which is obviously under abundant. Another important point of single component IDM is that the direct search allowed parameter space beyond resonance (m H 0 550 GeV) have significant co-annihilation dependence with m H ± − m H 0 10 GeV and m A 0 − m H 0 10 GeV.

Coupled Boltzmann equations and direct search
In presence of two DM components (φ and H 0 ) DM-DM conversion plays a crucial role. The heavier DM can annihilate to the lighter component and thus contribute to the freezeout of heavier DM. The conversion processes are shown in figure 8, which shows that they are dictated by four point contact interactions as well as by Higgs portal coupling. It is clear that H ± , A 0 are not really DM, but belongs to dark sector, hence annihilation to them is broadly classified within DM-DM conversion. More importantly none of them contribute to direct search. The two component DM set up therefore requires following parameters for analysis: There are two self interacting quartic couplings present in the model; namely λ Φ and λ φ , which do not play an important role in DM analysis, but appear in vacuum stability constraint that we discuss later. When m φ > m H 0 , m H ± , m A 0 , then φ can annihilate to all possible IDM components. The dominant s-wave DM-DM conversion cross-sections (σv) of φ in non-relativistic ap-JHEP03(2020)090 proximation are given by: . On the other hand, when m φ < m H 0 , the conversion process will be as H 0 H 0 (A 0 ) → φφ or H + H − → φφ. The corresponding cross-sections can easily be gauged from eq. (5.2).
The evolution of DM number density for both components (φ and H 0 ) in early universe as a function of time is obtained by coupled Boltzmann equations (CBEQ) as described in eq. (5.3): We can clearly spot DM-DM conversion contributions in second and third lines of each equation, which actually make the two equations 'coupled'. The freeze-out of two component DM is therefore obtained by numerically solving the above CBEQ and yields relic density (for a detailed discussion see for example [39]). The total relic density (Ω DM ) will then have contributions from both DM components as: Now let us turn to direct search of two component DM set up. Both the DM candidates can be detected through the spin independent (SI) direct detection (DD) processes through t-channel Higgs mediation as depicted in figure 9. The SI DD cross section for H 0 (σ H 0 ) and for φ (σ φ ) at tree level turn out to be [20,39]:  [95,96] and m N = 0.939 GeV represents nucleon mass. Importantly, the effective direct search cross-section for each individual component is modified by the fraction with which it is present in the universe, 2 given by It is pertinent to note here that next to leading order (NLO) corrections to SI direct detection cross section can turn dominant for some region of parameter space in case of IDM [98,99]. The NLO corrections arise mainly from one loop triangular and box diagrams mediated by gauge bosons and also from Higgs vertex corrections. In [99], authors extend the analysis of [98] to include gluon contribution at two loop and also contribution of IDM self coupling. In our case, the self coupling of IDM can be considered to be sufficiently small (as it is not playing any significant role in DM phenomenology) and hence we can ignore the effect of it. As shown in [99], in presence of loop effects, the coupling (λ L ) that enters into σ eff H 0 gets replaced by λ eff L = λ L + δλ L where δλ L takes into account the loop contributions effectively. Following this, we rewrite σ eff H 0 as where κ is the enhancement factor to SI DD cross section of IDM candidate due to loop correction. It is shown in [98], κ can take different values depending on λ L . For smaller values of λ L , κ turns bigger; for example, with 0.001 λ L < 0.01 the correction turns out to be 100 κ > 10; for 0.01 λ L < 0.1 → 10 κ > 1 and κ 1 for λ L 0.1 for single component IDM relic density allowed points. Here we follow similar scaling and take κ ∼ 0.1 λ L to accommodate NLO corrections for IDM direct search in our analysis. To obtain relic density and tree level DD cross sections of both the DM candidates numerically, we have used MicrOmegas [100]. It is noteworthy, that version 4.3 of MicrOmegas is capable of handling two component DM and we have cross-checked the solution from the code to match very closely to the numerical solution of CBEQ in eq.

Role of DM-DM conversion in relic density
We first study the variation of relic density with respect to DM mass and other relevant parameters to extract the importance of DM-DM conversion in this two-component set up before elaborating on the relic density and direct search allowed parameter space of the model.
In figure We see that for m H 0 < m φ , relic density of H 0 changes significantly with the variation of λ c . We also see that λ c = 0 (in two component set-up) is way above the pure IDM case yielding a large relic density. Now, with slight increase in λ c = 10 −5 , relic density goes further up and then reduces significantly for larger λ c = 0.01, 0.2. The interesting point is that the case of λ c = 0.2 lies very close to the pure IDM case.

JHEP03(2020)090
It is therefore evident that the presence of the second DM component φ plays an important role in Ω H 0 h 2 through the coupling λ c . For m H 0 < m φ , relic density of the heavier component φ can easily inherit the annihilation to other DM components {X, Y } = {H 0 , A 0 , H ± } in addition to SM as [39,102]: where in the last step, we have used g * = 106.7 and x f = 20 [102]. However it is difficult to envisage the relic density for the lighter DM component due to such DM-DM conversion. This can be understood from the CBEQ for the two component DM as in eq. (5.3). Let us define the following notations first: as earlier. The CBEQ with above notation, turns out to be (assuming m φ > m H 0 ): In eq. (5.9), we neglected the equilibrium number densities as they are tiny near freeze-out, where the dynamics is under study. With λ c = 0, and λ φh = 10 −5 , annihilation cross-sections for φ, ( σv φ φ→SM SM and σv φ φ→X Y ) are very small. Hence φ freezes out early and the number density of φ turns out to be large since Now, it is easy to appreciate that with λ c = 0, and λ φh = 10 −5 , annihilation of φ to SM is larger than conversion to other DM (H 0 ), i.e. σv φφ→ SM SM σv φφ→XY . However due to large φ abundance 3 (n φ ∼ 0.1 cm −3 ), F φΦ becomes comparable with F H 0 . As these two terms (F φΦ and F H 0 ) appear in the evolution of n H 0 (eq. (5.9)) with opposite sign, it is quite evident that effective annihilation cross-section for H 0 becomes small and hence n H 0 after freeze out turns out to be much larger than the pure IDM case. Next let us consider non zero but small λ c (= 10 −5 ). Then σv φφ→ X Y increases compared to the earlier case of λ c = 0. However due to smallness of the coupling λ c , this does not make any significant change in the number density of φ and n φ ∼ 0.1 cm −3 remains the same (as φ φ → SM SM still dominantly contributes to the total annihilation cross section of φ). Therefore, F φΦ increases and reduces the separation with F H 0 . Hence the effective annihilation cross-section for H 0 turns out to be even smaller than λ c = 0 case. Therefore, for λ c = 10 −5 relic density increases further than that of λ c = 0. For larger value of λ c = 0.01, contribution from DM-DM conversion, σv φφ→ X Y significantly rises and therefore JHEP03(2020)090 the number density of n φ drops to n φ ∼ 10 −5 cm −3 , and therefore F φΦ becomes much smaller than F H 0 . This increases the effective annihilation for H 0 and reduces relic density. This trend continues for higher values of λ c and eventually leads to a vanishingly small F φΦ to closely mimic the case of single component IDM. The case of m H 0 > m φ can also be understood from the CBEQ in this limit: For m H 0 > m φ , relic density of H 0 can be written simply as The annihilation to SM ( σv H 0 X→ SM SM ) due to gauge coupling is much larger than the conversion cross-section σv H 0 X→φφ for all the choices of λ c , so we do not find any distinction between all those cases.
• Top right of figure 10: In the top right panel, the same parameter space is used to show the variation of Ω φ h 2 with respect to m H 0 . The dynamics is much simpler With larger λ c , the conversion to other DM (F φΦ ) becomes larger and relic density drops accordingly. For m φ < m H 0 , we see from eq. (5.10), that there is a competition between F φ and F Φφ . With increasing m H 0 , σv H 0 X→φφ decreases, therefore F Φφ decreases and eventually it becomes vanishingly small for m H 0 400 GeV. The equation for n φ then becomes equivalent to the single component case of φ where m H 0 is no more relevant for Ω φ h 2 .
• The bottom panel figures of figure 10 essentially indicate that with larger λ φh , annihilation of φ to SM becomes large, resulting a smaller n φ after freeze out. Therefore F φΦ in eq. (5.9) turns insignificant. On the other hand, F Φφ also becomes smaller than F φ in eq. (5.10). Together, relic density of the lighter DM component is not affected by the presence of a heavy DM component.
Before we move on, let us summarise the outcome of figure 10. We see here that relic density of lighter DM component is affected by the heavier one, when the annihilation cross-section to SM is tiny.
In figure 11. we plot the variation of IDM relic density as function of DM mass m H 0 for different choices of ∆m = m A 0 − m H 0 : {1 − 50} GeV shown by different coloured lines. The other parameters are kept fixed and mentioned explicitly in the figure insets. Let us first focus on the left panel plot. We notice that for m H 0 < m W , with larger ∆m = m A 0 − m H 0 , relic density becomes larger. Keeping in mind that in this case (m H 0 < m W ), co-annihilations (see figure 7) dominate the total σv , this can be explained from the usual convention of co-annihilation cross-section being reduced by Boltzmann suppression (∼ JHEP03(2020)090 e −∆m/T ) with larger splitting (∆m). However, for m H 0 > m W , the phenomena seems to be reversed, i.e., with larger mass splitting, relic density decreases. Note that for m H 0 > m W , DM annihilations (in particular, annihilations into gauge bosons) become dominant. It turns out that the contribution of σv H 0 H 0 →W + W − or ZZ increases with larger ∆m (for the relevant cross-section, see [41]). This happens due to presence of the destructive interference between the s and t channel contributions. With larger ∆m, the destructive interference reduces and hence the annihilation cross section increases which in turn lowers the relic density of H 0 in this region. The same feature prevails in the right panel figure. Here, for a larger λ c , the relic density of H 0 goes further down due to annihilation to φ beyond m H 0 > m φ .
6 Relic density and direct search allowed parameter space One of the important motivations of this analysis is to study interacting multi-component DM phenomenology for DM mass lying between 80 ≤ m DM ≤ 500 GeV for both the components in view of relic density (0.1166 ≤ Ω DM h 2 (≡ Ω H 0 h 2 + Ω φ h 2 ) ≤ 0.1206 [14]) and direct detection (XENON 1T [26]) bounds. We would like to recall that for the individual scenarios, none of the DM components satisfy relic and direct detection bounds simultaneously within this mass range (see section 4). Now in order to find a consistent parameter space in the set up, we perform a numerical scan of the relevant parameters within the specified ranges as mentioned below.
We also note here that as λ L and λ φh enters into direct search cross-sections for H 0 and φ DMs respectively, we keep those couplings in a moderate range, while λ c governs DM-DM interactions, but do not directly enter into direct search bounds, therefore we choose a JHEP03(2020)090 larger range (with natural values) for scanning λ c . λ φ = 0.001 and λ Φ = 0.001 are kept fixed throughout the analysis since those are not relevant for DM phenomenology. There are two possible mass hierarchies for the two-component DM set up relevant for phenomenological analysis: (i) m φ > m H 0 and (ii) m φ ≤ m H 0 , which we address separately below.
Primarily, in such a scenario, the main physics arises due to annihilation of φ to H 0 , on top of their individual annihilation to SM to govern the freeze-out.
In figure 12, we show relic density allowed parameter space for the model in terms of different relevant parameters. In top left panel of figure 12, we have shown the relic density satisfied (considering contribution from both φ and H 0 components) points in m φ − λ φh plane for different ranges of λ c , as depicted in the figure with different colour codes. As seen from the plot, for a fixed m φ , there is a maximum λ φh . All possible values less than the maximum λ φh is also allowed subject to different choices of λ c . The larger is λ c , the smaller is the required λ φh thanks to the conversion of φ → H 0 to yield relic density. It is also noted that for large λ c 0.5, as the DM-DM conversion is very high, the DM mass m φ has to lie in the high mass region ( 400 GeV) to tame the annihilation cross-section within JHEP03(2020)090 acceptable range. To summarise, this figure shows that due to the presence of second DM component, much larger parameter space (actually the over-abundant regions of the single component framework, compare figure 5) is allowed. Top right panel figure shows the relic density allowed parameter space in m H 0 − λ L plane again for different ranges of λ c as in the left plot. It naturally depicts that λ L is insensitive to m H 0 as the annihilation and co-annihilation cross-section of H 0 is mainly dictated by gauge interaction. However, we see a mild dependence on λ c , such that when λ c 0.5, the DM mass m H 0 has to be heavy 400 GeV. This is because with large λ c , DM-DM conversion is large; to achieve relic density of correct order, m φ requires to be large and the conversion can only be tamed down by phase space suppression, i.e. by choosing H 0 mass as close as possible to φ mass (m H 0 ∼ m φ ). Bottom left figure correlates the DM masses to obtain correct density within m φ > m H 0 . We see that for small λ c , particularly with higher m φ , large m H 0 values are disfavoured in order to keep the DM-DM conversion in the right order. While for large λ c , mass degeneracy is required (m H 0 ∼ m φ ) to tame the DM-DM conversion. Bottom right figure shows the relative contribution of relic density of the two DM components. First of all, this shows that φ contributes with larger share of relic density, while the relic density of H 0 can at most be 40% of the total. This is natural from the perspective of IDM to be the lighter DM component, as we know it is always suppressed below 500 GeV (with a maximum contribution 40%) and the conversion channels do not alter IDM's annihilation for m φ > m H 0 . Although for 0.01 < λ c 0.1, individual relic abundance of IDM could be large for very small (and fine tuned) λ φh ∼ 10 −5 , as shown in upper left of figure 10; with the increase of λ φh , this enhancement vanishes. Since we scan the parameter space where λ φh 0.001, such effect can not be found in figure 12. For small λ c , contribution from H 0 is smaller, as relic density contribution from φ gets larger due to small DM-DM conversion. With larger λ c 0.1, the DM-DM conversion for φ becomes larger and therefore the relic density of φ can easily span between 60 − 100%. With very high λ c ∼ 1, DM-DM conversion becomes too large, therefore to keep relic density in the correct ballpark, the DM mass (m φ ) has to be heavy and almost degenerate with the heavier DM (m H 0 ∼ m φ ). λ φh in such cases, requires to be very small, which are validated by some dark blue points with Ω φ h 2 /Ω DM h 2 ∼ 60%.
Next in figure 13, we show the SI direct detection cross section for H 0 at relic density satisfied points in σ eff H 0 − m H 0 plane. Direct search cross-section at tree level is shown by light black points and NLO effects as prescribed in eq. (5.6), is shown by orange points. For comparison, we also point out the most recent XENON 1T experimental bound in dotted red line. Upper left panel of figure 13 shows that for 0.001 λ L < 0.01, the one loop correction to σ eff H 0 is really significant considering κ ∼ 0.1 λ L and rules out some of the DD bound satisfied points at tree level. Similar kind of feature can be found for 0.01 λ L < 0. similar to figure 12. We have already mentioned that spin independent (SI) DM-nucleon cross-section depends on square of Higgs portal couplings of the respective DM candidates, λ φh for φ and λ L for H 0 scaled by a pre-factor of the relative number density Ω i h 2 /Ω DM h 2 (i = φ, H 0 ). Since in this two component scenario, the dominant contribution is coming from φ DM, the pre-factor Ω φ h 2 /Ω DM h 2 ∼ 1. On the other hand, for H 0 the pre-factor is small, Ω H 0 h 2 /Ω DM h 2 < 1, and will help H 0 reducing the effective direct search crosssection. Therefore, portal coupling λ φh is tightly constrained from XENON 1T bound to  figure). In absence of a lighter mode, the relic density of φ is essentially governed by its annihilation to SM and due to small conversion cross-section the production of φ is also not large enough to change the conclusion. However, the situation changes significantly with larger λ c (cyan and blue points), where we see again that the overabundant region of the single component scenario is getting allowed by relic density. In order to understand this let us remind ourself of the CBEQ for m φ < m H 0 as depicted in eq. (5.10). In particular, the number density of φ is dictated bẏ n φ +3Hn φ −F φ +F Φφ . With larger λ c and larger conversion, F Φφ increases to reduce the effective F φ that determines the relic of φ. Therefore, to keep the relic density of φ to cor-JHEP03(2020)090 rect order, F φ has to increase. Now recall, F φ = σv φφ→SM SM (n φ ) 2 ∼ 1/ σv φφ→SM SM as n φ ∼ 1/ σv φφ→SM SM . Therefore, to increase F φ , one has to reduce the annihilation cross-section σv φφ→SM SM . This is possible by reducing λ φh as we see here in the plot. Next let us discuss the top right figure. This figure in m H 0 − λ L plane essentially depicts that with larger λ c , larger m H 0 is favoured to tame the DM conversion as well as annihilation cross-section to keep the relic within limit. The dependence however is not that much significant due to the presence of large number of co-annihilation channels which remain unaffected by λ c . In the bottom left panel, mass correlation has been plotted and carries no information. Lastly, bottom right figure shows the relative relic density contributions of these two components. It is well understood that an additional channel for annihilation of H 0 only reduces the possibility of bringing Ω H 0 h 2 in the correct ballpark due to already existing gauge mediated annihilation and co-annihilation channels. Therefore, for small λ c , it is still possible to get a contribution from Ω H 0 h 2 ∼ 40%, but that becomes harder with large λ c , where the relic density contribution of H 0 is further limited to Ω H 0 h 2 ∼ 20%.
Direct search constraints from XENON 1T 2018 on the relic density allowed points are shown in figure 16. To emphasise again, the demand of these plots are simultaneous satisfaction of XENON1T limit for both DM components. The main outcome of this plot is to see the absence of small λ c points (red dots) upto λ c ∼ 0.30. This is simply due to JHEP03(2020)090 the fact that, with small λ c , the required λ φh is high enough for φ DM to be discarded by XENON1T data. The other important feature is that with larger λ c , larger DM masses are favoured. Lastly a very important conclusion comes from the bottom panel figure in the mass correlation plot. This shows, as only high λ c region is allowed, the conversion of H 0 to φ still needs to be restricted and therefore the mass difference between the two DM components (m H 0 − m φ ) has to be very small. These features are all distinct from that of m φ > m H 0 region. So far our discussion has been focused on DM mass region m W ± ≤ m H 0 , m φ ≤ 500 GeV. But if m H 0 < m W ± or m φ < m W ± , while other DM mass is heavier than m W ± , the only region available for lighter DM with mass < m W ± are the Higgs and Z resonance regions: It is important to remind that the resonance regions are already available in absence of second DM component and therefore brings no new phenomenological outcome.

Indirect searches
It is possible to probe DM signals in the context of different indirect detection experiments as well. The experiments look for astrophysical sources of SM particles produced through JHEP03(2020)090 DM annihilations or via DM decays. Amongst these final states, photon and neutrinos, being neutral and stable can reach indirect detection probes without getting affected much by intermediate regions. For WIMP type DM, emitted photons lie in the gamma ray regime that can be measured at space-based telescopes like the Fermi-LAT or ground-based telescopes like MAGIC [104]. The photon flux in a specific energy range is written as where r(b, l, x) is the distance of the DM halos from the galactic center. Galactic coordinates are defined by b, l and ρ(r) represents the DM density profile. Eq. (7.1) is further rewritten as where J = dxρ 2 (r(b, l, x)) is conventionally known as J-factor and LOS means line of sight. Now from the observed Gamma ray flux produced due to DM annihilations, one can constrain the DM annihilation into different charged final states like µ + µ − , τ + τ − , W + , W − and b + b − . It is pertinent to mention here that along with the continuous gamma ray spectrum, gamma ray lines originated from IDM annihilations are also of special interest [105]. It has been shown in [106] that such a feature, characteristic of inert doublet dark matter (IDM) model, can indeed be important for DM mass above TeV where Sommerfeld enhancement becomes relevant. Since we mostly focus on the intermediate mass range of IDM in this study (∼ 80 -500 GeV), aforementioned constraint is extremely subdued.
In the multicomponent dark matter set up as we have here, it is expected that the γ ray flux will be the sum of the contributions from both φ and H 0 . Then we can write the total γ ray flux can be written as: where Φ φ F and Φ H 0 F are the individual contribution to Φ T F from φ and H 0 . Following eq. (7.2), it would be natural to write where X, X denote the annihilation products and m R indicates some effective mass scale. One can further write the effective annihilation cross-section in eq. (7.4) considering m R to be the reduced mass of the two component system as It turns out that for real singlet scalar dark matter indirect searches does not pose strong bound [18,107]; while for single component IDM candidate, the indirect search rules out the region of m H 0 < 400 GeV due to large thermal average cross section of H 0 H 0 → W + W − channel [106][107][108][109]. Hence in our proposed two component DM set up JHEP03(2020)090  we shall mostly focus on the σv DM,DM→W + W − and compare it with the indirect search experimental bound from Magcic+Fermi Lat bound [104]. In left panel of figure 17, we show the contribution of φ DM in σv eff W + W − for all relic and SI DD cross section satisfied points and compare it with the limit coming from Fermi Lat experiment [104]. Similar comparison is drawn for IDM (H 0 ) in right panel of figure 17. Finally in figure 18, we show the order of magnitude of total σv eff W + W − for all relic density and SI DD (at NLO) satisfied points and find it to lie well below the upper limit imposed by Fermi Lat experiment. The relative abundance of individual relic densities and ratio of reduced mass to the DM mass helps to evade the indirect search bound on σv DM,DM→W + W − successfully. We have also confirmed that the bounds on σv eff for bb, µ + µ − , τ + τ − and other annihilation final products from Magic+Fermi Lat experiment [104] are easily satisfied in our set up.

JHEP03(2020)090 8 Electroweak vacuum stability and high scale perturbativity in presence of RH neutrino and DM
One of the motivations of this work is to show that the presence of right handed neutrinos to yield correct neutrino masses in presence of multipartite DM. Although the neutrino sector considered here seems decoupled from the dark sector, is not completely true. The effect of the RH neutrinos alter the allowed DM parameter space when the model is validated at high scale. As already stated before, we employ type-I seesaw mechanism to generate the light neutrino mass, for which three RH neutrinos are included in the set up.

Neutrino Yukawa couplings
We first describe the strategy in order to study their impact on RG evolution. For simplicity, the RH neutrino mass matrix M N is considered to be diagonal with degenerate entries, i.e. M N i=1,2,3 = M R . It is to be noted that in the RG equation, Tr[Y † ν Y ν ] enters. In order to extract the information on Y ν , we use the type-I seesaw formula for neutrino 2M R . Then, naively one would expect that large Yukawa couplings are possible with even heavier RH neutrino masses. For example with M R ∼ 10 14 GeV, Y ν comes out to be 0.3 in order to obtain m ν 0.05 eV. However, contrary to this, it can be shown that even with smaller M R one can achieve large values of Tr[Y † ν Y ν ], but effectively keeping Y T ν Y ν small, using a special flavor structure of Y ν through Casas-Ibarra parametrization [58]. Note that our aim is to study the maximum effect coming from the right handed neutrino Yukawa, i.e. with large Tr[Y † ν Y ν ], on EW vacuum stability. For this purpose, we use the parametrisation by [110] and write Y ν as where m d ν is the diagonal light neutrino mass matrix and U PMNS is the unitary matrix diagonalizing the neutrino mass matrix m ν such that m ν = U * PMNS m d ν U † PMNS . Here R represents a complex orthogonal matrix which can be written as R = Oexp(iA) with O as real orthogonal matrix and A as real antisymmetric matrices. Using above parametrisation, then one gets, Note that the real antisymmetric matrix A however does not appear in the seesaw expression for neutrino mass as m ν = Y T ν Yν v 2 2M R . Therefore with any suitable choice of A, it would actually be possible to have relatively large Yukawa even with light M R .

β functions and RG running
To study the high scale validity of this multi-component DM model with neutrino mass (including perturbativity and vacuum stability), we need to consider the RG running of the associated couplings. The framework contains few additional mass scales: one extra scalar singlet, one inert doublet and three RH neutrinos with mass M N i=1,2,3 (= M R ). Hence

JHEP03(2020)090
in the study of RG running, different couplings will enter into different mass scales. In DM phenomenology, we have considered the physical masses of DM sector particles within ∼ 500 GeV. Then for simplicity, we can safely ignore the small differences between the dark sector particle masses and identify the single additional mass scale as m DM . On the other hand, RH neutrinos are considered to be heavier than the scalars present in the model (M R > m DM ). Hence, for running energy scale µ > m DM , we need to consider the couplings associated to DM sector in addition to SM while for µ > M N i=1,2,3 , the neutrino couplings will additionally enter into the scenario. Below we provide the one loop β functions for the additional for the couplings involved in our model, when µ > M R .
β functions of gauge couplings [111] β g 1 = 21 5 g 3 1 , β functions of Yukawa couplings [111,112]: β functions of quartic scalar couplings [111] β λ H = 27 200 g 4 1 + 9 20 g 2 1 g 2 2 + 9 8 g 4 2 − The above β functions are generated using the model implementation in the code SARAH [113]. The running of λ H as in eq. (8.5) is influenced adversely mostly by top Yukawa JHEP03(2020)090 coupling y t ∼ O(1) and right handed neutrino Yukawa coupling as Tr[Y † ν Y ν ]. On the other hand, multiparticle scalar DM couplings present in the model help in compensating the strong negative effect from y t and Tr[Y † ν Y ν ]. To evaluate Tr[Y † ν Y ν ], we employ eq. (8.2). As an example, let us consider magnitude of all the entries of A to be equal, say a with all diagonal entries to be zero. Then using the best fit values of neutrino mixing angles and cosnidering the mass of lightest neutrino zero, one can find forM R = 1 TeV, Tr[Y † ν Y ν ] can be as large as 1 with a = 8.1 [110,114]. Equipped with all these β functions and strategy to estimate the relevant couplings present in them, let us now analyse the SM Higgs vacuum stability at high energy scale. Below we provide the stability and metastability criteria.
Stability criteria. The stability of Higgs vacuum can be ensured by the condition λ H > 0 at any scale. However, we have multiple scalars (SM Higgs doublet, one inert doublet and one gauge real singlet) in the model. Therefore we should ensure the boundedness or stability of the entire scalar potential in any field direction. This can be guaranteed by using the co-positivity criterion in eq. (3.1). Note that once λ H becomes negative the other copositivity conditions no longer remain valid.
Metastability criteria. On the other hand, when λ H becomes negative, there may exist another deeper minimum other than the EW one. Then the estimate of the tunnelling probability P T of the EW vacuum to the second minimum is essential to confirm the metastability of the Higgs vacuum. Obviously, the Universe can be in a metastable state only, provided the decay time of the EW vacuum is longer than the age of the Universe. The tunnelling probability is given by [4,5], where T U is the age of the Universe, µ B is the scale at which tunnelling probability is maximized, determined from β λ H = 0. Therefore, solving above equation, we see that metastable Universe requires [4,5]: As noted in [4], for µ B > M P (Planck Scale, M P = 1.22×10 19 GeV), one can safely consider λ H (µ B ) = λ H (M P ). One should also note, that even with meatstability of Higgs vacuum, the instability in other field direction may also occur. The conditions to avoid the possible instability along various field directions for λ H < 0 are listed below [115]: (i) λ Φ > 0 to avoid the unboundedness of the potential along H 0 , A 0 and H ± directions, (ii) λ 1 > 0 to ensure the stability along a direction between H ± and h, (iii) λ L > 0 to ensure the stability along a direction between H 0 and h, (iv) λ S > 0, to avoid unboundedness of potential along a direction between A 0 and h, (v) λ φ > 0, otherwise the potential will be unbounded along φ direction, (vi) λ φh > 0, to ensure the stability along a direction between φ and h.

JHEP03(2020)090
BPs  Now to begin the vacuum stability analysis in the present multi-component DM scenario, we first choose two benchmark values of DM masses (φ and H 0 ) which satisfy both the relic density and direct detection bounds. These along with corresponding values of other relevant parameters are denoted by the set of two benchmark points, BP1 and BP2, as tabulated in table 2. These parameters are mentioned in the caption of table 2 for both benchmark points. We also show the value of electroweak parameters and µ γγ for these two benchmark points in table 3. We run the two loops RG equations for all the SM couplings and one loop RG equations for the other relevant couplings in the model from µ = m t to M P energy scale. We use the inital boundary values of all the SM couplings as provided in [4]. The boundary values have been evaluated at µ = m t in [4] by taking various threshold corrections and mismatch between top pole mass and M S renormalised couplings into account.
In figure 19, we show the running of λ H for BP1 as a function of energy scale µ. The left panel shows running of λ H for different values of RH neutrino mass M R , considering top quark mass m t = 173.1 GeV, Higgs mass m h = 125.09 GeV and Tr[Y † ν Y ν ] = 0.5. As it is visible that for low value of M R ∼ 10 4 GeV, λ H enters into unstable region at very early stage of its evolution (blue line in the figure). This has happened as the scalar couplings in β λ H (eq. (8.5)) are not sufficiently large to counter the strong negative impact of Tr[Y † ν Y ν ] which brings down the λ H shraply towards unstable region. On contrary, for large value of M R ∼ 10 14 GeV, the effect of Y ν in β λ H starts at a comparatively larger energy scale than the earlier case. As a consequence, although λ H becomes negative at some high energy scale, it stays in metastable region till M P energy scale (violate line). Green region here describes stable, white region describes metastable (see eq. (8.8)) and the red region indicates unstable solution for the potential. For comparison, we also display the evolution of λ H (black dotted curve) in absence of RH neutrinos in the theory with the inclusion of scalar couplings (related to DM) only. This clearly shows that in absence of RH neutrinos, EW vacuum could be absolutely stable till M P energy scale. In right panel of figure 19, we study the evolution of λ H considering 2σ uncertainty of measured top mass m t , keeping JHEP03(2020)090     One important distinction between the left and right panel figures, corresponding to two different benchmark points, is that BP2 has significantly larger parameter space available from high scale stability. This is because of the larger values of λ 1,2,3 parameters used in BP2 compared to BP1 (see table 4 for details). Therefore, it can be concluded, that larger is the mass splitting in IDM sector, the more favourable it is from the stability point of view. However there is an important catch to this statement, which we will illustrate next.
In figure 24 we plot the running of all the individual couplings present in the set up. We see that (fortunately) for the two benchmark points used in the analysis, we are still okay with the perturbative limit at the high scale, i.e. all the couplings obey the perturbative limit, |λ i | < 4π, |Tr[Y † ν Y ν ]| < 4π at Planck scale. However, for BP1, with M R = 10 8 GeV, and Tr[Y † ν Y ν ] = 0.5 as shown in the left panel yields unstable solution to EW vacuum with λ H running negative. On the other hand, BP2 with same choices of right handed neutrino mass and Yukawa yields a stable vacuum. Therefore, once again we find that larger splitting in IDM sector is more favourable for EW vacuum stability, as we have also inferred before from figure 23. But, it turns out that as larger splitting in the IDM sector also uses JHEP03(2020)090 larger values of λ 1,2,3 , it is possible, that many of those points become non-perturbative i.e. |λ i | > 4π when run upto Planck scale. We will show in the next section, that this very phenomena disallows many of the relic density and direct search allowed parameter space of the model. Another point to end this section is to note that our available parameter space from DM constraints heavily depend on large DM-DM conversion, which naturally comes from large choices of the conversion coupling λ c . It is natural, that some of those cases will also be discarded by the perturbative limit |λ c | < 4π, when we evaluate the validity of the model at high scale.

Allowed parameter space of the model from high scale validity
Finally, we come to the point where we can assimilate all the constraints together, from DM constraints to high scale validity. In order to do that, we choose relic density and direct search at NLO allowed parameter space of the model as discussed in section 6 and impose the high scale validity of the model till some energy scale µ by demanding: • Stability of the scalar potential (eq. (3.1)) determined by satisfying copositivity conditions for any energy scale µ.
• Non-violation of perturbativity and unitarity of all the relevant couplings present in the model as defined in eqs. (3.4) and (3.5).
Note that the high scale validity of the models does not depend on the structure of mass hierarchy of the DM candidates (i.e. on the sign of mass difference m H 0 − m φ ). To study the EW vacuum stability we demand the positivity of λ H at each energy scale from EW to high scale µ. As evident from β λ H in eq. (8.5), a particular combination of the scalar couplings in the form of 2λ 2 1 + λ 2 2 + λ 2 3 + 2λ 1 λ 2 + 1 2 λ 2 φh determines the positivity of λ H during its running. However the factor λ φh 0.1 is strongly constrained from the SI DD cross section bound for m φ < 500 GeV. Hence, we can assume safely that the factor 2λ 2 1 + λ 2 2 + λ 2 3 + 2λ 1 λ 2 without λ φh effectively determines the stability of Higgs vacuum in relic density direct search and indirect search allowed parameter space. It turns out that when λ H > 0. all other copositivity conditions for all relic and DD cross section satisfied points in our model stays positive from µ = m t to µ = M P .
In figure 25, we constrain relic density direct search and indirect search allowed points of the model in 2λ 2 1 + λ 2 2 + λ 2 3 + 2λ 1 λ 2 − λ 1 plane to additionally satisfy perturbativity and vacuum stability conditions following eqs. (3.1)-(3.5). Orange points satisfy relic density, DD bounds and indirect search bound, while the green, blue and red points on top of that satisfy perturbativity and vacuum stability conditions upto high energy scales µ = 10 10 GeV, 10 16 GeV and 10 19 GeV respectively. This very figure essentially shows that all those points with either small values of λ 1 (i.e. small m H ± − m H 0 )are discarded due to stability of EW vacuum, while those with large λ 1 (i.e. large m H ± − m H 0 )are discarded by perturbative limits of the coupling at high scale.
In figure 26 we study the correlation between the individual scalar couplings to satisfy the DM constraints, perturbativity limits and vacuum stability criteria. In left panel of figure 26, we show the DM relic density and DD cross section satisfied points in λ 2 − λ 3 JHEP03(2020)090  plane for different values of λ 1 . In right panel we first identify the relevant parameter space in the same plane which satisfy the DM constraints, the perturbativity bound and vacuum stability criteria till EW energy scale. Then we further impose perturbativity bound and vacuum stability conditions considering the high scale as µ = 10 10 GeV, 10 16 GeV and 10 19 GeV in addition to the DM constraints. It is seen that the lower portion of the available parameter space gets discarded by vacuum stability or high value of λ c while perturbativity bounds constrain the higher values of the couplings. In this plot also we kept M R = 10 8 GeV, Tr[Y † ν Y ν ] = 0.5 with m t = 173.1 GeV. We must also note that with larger M R and smaller Yukawa Tr[Y † ν Y ν ], we could obtain a larger available parameter space from high scale validity.
In figure 27, The orange points are corresponding to relic and direct search allowed parameter space at EW scale. We see that larger mass difference between inert higgs components, ∆m − ∆M which are related with quartic couplings, λ 1,2,3 (see eq. (2.7)), are discarded from perturbativity conditions mentioned in eq. (3.4). While the small mass differences between inert componets are also excluded from stability criteria of Higgs potential.
Till now, while discussing the effect of stability and high scale validity of the proposed set up, we have considered fixed right handed neutrino mass, M R = 10 8 GeV and corresponding Yukawa coupling Tr[Y † ν Y ν ] = 0.5. For sake of completeness we extend our study for few different values of RH neutrino mass and Tr[Y † ν Y ν ] and find out the allowed ranges of the relevant parameters considering both the hierarchies m φ > m H 0 and m φ ≤ m H 0 separately. We note the corresponding results in table 9 and table 10   LHC in presence of a second scalar singlet DM component. It is also worth mentioning here that the real scalar singlet, which interacts with SM only through Higgs portal coupling, does not have any promising collider signature excepting mono-X signature arising out of initial state radiation (ISR), where X stands for W, Z, jet. Such signals are heavily submerged in SM background due to weak production cross-section of DM in relic density and direct search allowed parameter space [51]. The charge components H + , H − of inert DM can be produced at LHC via Drell-Yan Z and γ mediation as well as through Higgs mediation. Further decay of H ± to DM (H 0 ) and leptonic final states through on/off shell W ± yields hadronically quiet opposite sign di-lepton plus missing energy (OSDL+ / E T ), 4 as shown in left panel of figure 28. In this study, we focus on this particular signal of inert dark matter as detailed below: Around m H ± ∼ m h /2, there is a sharpe fall of production cross-section. This is because, for m H ± ≤ m h /2, there is a significant contribution arising from Higgs production and its subsequent decay to the charged scalar components, which otherwise turns into an off-shell propagator to yield a subdued contribution to Drell-Yan production. Following [116], a conservative bound on the charge scalars is applied here as m H ± ≥ 70 GeV, as indicated in the r.h.s. plot of figure 28. BPs Table 4. DM masses, quartic couplings, relic densities and spin independent effective DM-neucleon cross-section of selected benchmark points for collider study. All benchmark points chosen here have m H ± − m H 0 < m W ± for off-shell production of W ± . The maximum scale (µ) of Higgs vacuum satability and peturbativity in presence of right handed neutrinos are also noted. All masses and scales are in GeV. BPs Table 5. DM masses, quartic couplings, relic densities and spin independent effective DM-neucleon cross-section of selected benchmark points for collider study. All benchmark points chosen here have m H ± − m H 0 > m W ± for on-shell production of W ± . The maximum scale (µ) of Higgs vacuum satability and peturbativity in presence of right handed neutrinos are also noted. All masses and scales are in GeV.

JHEP03(2020)090
We next choose a set of benchmark points (BPs) allowed from DM relic, direct search constraints as well as from Higgs invisible decay constraints for performing collider simulation, shown in

JHEP03(2020)090
The simulation technique adopted here is as follows. We first implemented the model in FeynRule [117] to generate UFO file which is required to feed into event generator Madgraph [118]. Then these events are passed to Pythia [119] for hadronization. All parton level leading order (LO) signal events and SM background events 5 are generated in Madgraph at √ s = 14 TeV using cteq6l1 [120] parton distribution. Leptons ( = e, µ) isolation, jet and unclustered event formation to mimic to the actual collider environment are performed as follows: (i) Lepton isolation. The minimum transverse momentum required to identify a lepton ( = e, µ) has been kept as p T > 20 GeV and we also require the lepton to be produced in the central region of detector followed by pseudorapidity selection as |η| < 2.5. Two leptons are separated from each other with minimum distance ∆R ≥ 0.2 in η − φ plane. To separate leptons from jets we further imposed ∆R ≥ 0.4.
(ii) Jet formation. For jet formation, we used cone algorithm PYCELL in built in Pythia.
All partons within a cone of ∆R ≤ 0.4 around a jet initiator with p T > 20 GeV is identified to form a jet. It is important to identify jets in our case because we require the final state signal to be hadronically quiet i.e. to have zero jets.
(iii) Unclustered Objects. All final state objects with 0.5 < p T < 20 GeV and 2.5 < |η| < 5 are considered as unclustered objects. Those objects neither form jets nor identified as isolated leptons and they only contribute to missing energy.
The main idea is to see if the signal events rise over SM background. For that there are three key kinematic variables where the signal and background show different sensitivity. They are: • Missing Energy ( / E T ). The most important signature of DM being produced at collider. This is defined by a vector sum of transverse momentum of all the missing particles (those are not registered in the detector); this in turn can be estimated form the momentum imbalance in the transverse direction associated to the visible particles. Thus missing energy (MET) is defined as: where the sum runs over all visible objects that include the leptons, jets and the unclustered components.
• Transverse Mass (H T ). Transverse mass of an event is identified with the scalar sum of the transverse momentum of objects reconstructed in a collider event, namely lepton and jets as defined above. • Invariant mass (m ). Invariant mass of opposite sign dilepton hints to the parent particle, from which the leptons have been produced and thus helps segregating signal from background. This is defined as: The distribution of missing energy ( / E T ), invariant mass of opposite sign dilepton (m + − ) and transverse mass (H T ) for the BPs along with dominant SM background events are shown in figure 29 top left, top right and bottom panel respectively. All BPs depicted in figure 29 correspond to m H ± − m H 0 ≡ ∆M < m W ± where opposite sign di-lepton are produced from off-shell W ± mediator. All the distributions are normalised to one event.
Missing energy (as well as H T ) distributions of BPs (BPC1-BPC3) show that the peak of distribution for the signal is on the left of SM background. This is because the benchmark points are characterised by small ∆M , where the charged scalars and inert DM have small mass splitting. Therefore, such situations are visibly segregated from SM background by MET and H T distribution. Clearly, when ∆M becomes m W ± (for example, BPC3) the distribution closely mimic SM background. Therefore the signal events for this class of benchmark points can survive for a suitable upper / E T and H T cut while reducing SM backgrounds. It is important to take a note that OSDL events coming from ZZ background JHEP03(2020)090 naturally peaks at m Z in m distribution. Therefore, we use invariant mass cut in the Z mass window to get rid of this background. The situation is reversed for larger splitting between the charged scalar component with DM, i.e ∆M ≡ m H ± − m H 0 > m W ± corresponding to BPs (BPD1-BPD2) as in table 5. The distributions of / E T , m and H T therefore become flatter and peak of the distribution shifts to higher value as shown in figure 30. In such cases the signal events for large ∆M can be separated from SM background at a suitable lower end cut of / E T and H T .
Therefore the selection cuts used in this analysis are summarised as follows: • Invariant mass (m ) cuts: m < (m Z − 15) GeV and m > (m Z + 15) GeV.
BPs  Table 7. Signal cross-section for BPD1-BPD2 after the selection cuts are employed.
We next turn to signal and background events that survive after the selection cuts are employed. The signal events are listed in table 6 for BPC1-BPC3. Similarly, signal events for BPD1-BPD2 are listed in table 7 using a lower cut on / E T and H T . The signal cross-section and event numbers for this class of points are much smaller due to the fact that the charged scalar masses are on higher side as for all of the cases ∆m > m W independent of DM masses. However, the SM background events get even more suppressed with the set of / E T , H T cuts. The SM background cross-section and event numbers after cut flow is mentioned in table 8. Therefore, in spite of smaller signal events in this region of parameter space with BPD benchmark points, the discovery potential of the signal requires similar luminosity to that of BPC cases.
Finally we present the discovery reach of the signal events in terms of significance σ = S √ S+B , where S denotes signal events and B denotes SM background events in terms of luminosity. This is shown in figure 31. This shows that the benchmark points that characterise the two component DM framework, can yield a visible signature at high luminosity with L ∼ 500 fb −1 depending on the charged scalar mass and its splitting with DM for the case m H ± − m H 0 > m W ± .

Summary and conclusions
We have studied a two component scalar DM model in presence of right handed neutrinos that address neutrino mass generation through type I seesaw. The DM components are SM Bkg. (i) a singlet scalar and (ii) an inert scalar doublet, both studied extensively as single component DM in literature. We show that the presence of second component enlarges the available parameter space significantly considering relic density, direct search and indirect search constraints. In particular, the inert scalar DM will now be allowed in the so called 'desert region': {m W − 550} GeV. Also for singlet scalar, we can now revive it below TeV, which is otherwise discarded (except Higgs resonance) from direct search in single component framework. The results obtained for DM analysis crucially depends on DM-DM conversion, which have been demonstrated in details.
We also study the high scale perturbativity and vacuum stability of the Higgs potential by analysing two loop RGE β functions. This in turn puts further constraints on the available DM parameter space of the model. One of the important conclusions obtained are that the mass splitting of the charged scalar component to the corresponding DM JHEP03(2020)090 The presence of RHNs in the model not only helps us addressing the neutrino masses but also controls the high scale validity of the model parameters, for example, low ∆M regions. This is how the neutrino and dark sector constraints affect each other.
Inert Higgs having charged scalars have collider detectability. We point out that the collider search prospect of the charged components are not only limited to low DM masses (< m W ), but is open to a larger mass range in presence of the second DM component, even after taking the high scale validity constraints. We exemplified this at LHC for hadronically quiet dilepton channel with missing energy, where ∆M , turns out to be a crucial kinematic parameter, constrained from DM, high scale validity and neutrino sector. At LHC, due to its tt background, high ∆M regions can be segregated from SM background more efficiently at the cost of small production cross section. Low ∆M regions are more affected by SM background, although the signal cross-section is high. Therefore the model can be probed within High Luminosity reach of LHC. On the other hand, e + e − annihilation have better possibility to explore low ∆M regions absent tt background. Here, we would like to comment that there are several studies that have been done in this direction, but the high scale validity constraint may alter the conclusion significantly as we demonstrate.
Finally, we would like to mention that the analysis performed here, although focus on a specific model set up, but there are some generic conclusions that can be borrowed. For example, if the two DM components have sufficient interaction in between, the available parameter space will be enlarged significantly from both relic density and direct search.