Minimal Two-component Scalar Doublet Dark Matter with Radiative Neutrino Mass

TL;DR

This paper introduces a minimal extension of the Standard Model with two inert scalar doublets and a right-handed neutrino, enabling two-component dark matter and radiative neutrino mass generation, expanding viable dark matter mass range.

Contribution

It presents a novel two-component scalar doublet dark matter model with radiative neutrino masses, unifying dark matter and neutrino mass mechanisms in a minimal framework.

Findings

Two-component dark matter model allows for a wider mass range.

Inter-conversion between dark matter candidates enhances relic abundance.

Model successfully explains neutrino masses radiatively.

Abstract

We propose a minimal extension of the Standard Model to accommodate two-component dark matter (DM) and light neutrino mass. The symmetry of the Standard Model is enhanced by an unbroken $\mathbb{Z}_2 \times \mathbb{Z}'_2$ such that being odd under each $\mathbb{Z}_2$, there exists one right handed neutrino and one inert scalar doublet. Therefore, each of the $\mathbb{Z}_2$ sectors contribute to ($i$) light neutrino masses radiatively similar to the scotogenic models while ($ii$) the two neutral CP even scalars present in two additional inert doublets play the role of dark matters. Focusing on the intermediate range of inert scalar doublet DM scenario: $M_W \leq M_{\rm DM} \lesssim 500 \; {\rm GeV}$, where one scalar doublet DM can not satisfy correct relic, we show that this entire range becomes allowed within this two-component scalar doublet DM, thanks to the inter-conversion between the two DM candidates in presence of neutrino Yukawa couplings with dark sector.

Tables

Table 1. Table 1: New particle content of the model and their charge assignments.

Field	$SU{\left(3\right)}_c\times SU{\left(2\right)}_L\times U{\left(1\right)}_Y$	${\mathbb{Z}}_2$	$\mathbb{Z}_2^{}{}_{}{}^{\prime }$
${\eta }_1$	(1, 2,$\frac{1}{2}$)	-	+
${\eta }_2$	(1, 2,$\frac{1}{2}$)	+	-
$N_1$	(1, 1,$0$)	-	+
$N_2$	(1, 1,$0$)	+	-

Table 2. Table 2: Global fit 3 ​ σ 3 𝜎 3\sigma values of neutrino oscillation parameters Esteban et al. ( 2019 ) .

Parameters	Normal Hierarchy (NH)	Inverted Hierarchy (IH)
$\frac{\Delta m_{21}^2}{{10}^{-5}{\text{eV}}^2}$	$6.79-8.01$	$6.79-8.01$
$\frac{\left|\Delta m_{31}^2\right|}{{10}^{-3}{\text{eV}}^2}$	$2.427-2.625$	$2.412-2.611$
${\sin }^2{\theta }_{12}$	$0.275-0.350$	$0.275-0.350$
${\sin }^2{\theta }_{23}$	$0.418-0.627$	$0.423-0.629$
${\sin }^2{\theta }_{13}$	$0.02045-0.02439$	$0.02068-0.02463$
$\delta {(}^{\circ })$	$125-392$	$196-360$

Table 3. Table 3: The couplings λ L 1 ​ and ​ λ L 2 subscript 𝜆 subscript 𝐿 1 and subscript 𝜆 subscript L 2 \lambda_{L_{1}}~{}\rm{and}~{}\lambda_{L_{2}} are considered to be 0.01 with λ 12 = 0 subscript 𝜆 12 0 \lambda_{12}=0 .

BP	Masses of (in GeV)				${\Omega }_1{h}^2$	${\Omega }_2{h}^2$	${\sigma }_{1,eff}$[pb]	${\sigma }_{2,eff}$[pb]	${\sigma }_{pp\to H_1^{+}{H}_1^{-}}$[fb]	${\sigma }_{pp\to H_2^{+}{H}_2^{-}}$[fb]
	$H_1$	$H_1^{+}$	$H_2$	$H_2^{+}$
Our model	$300$	$301$	$350$	$351$	$0.0497$	$0.070$	$1.6\times {10}^{-11}$	$1.7\times {10}^{-11}$	3.86	2.09

Equations

$$-{\mathcal{L}}^{\mathrm{new}}=Y_{\alpha 1}{\bar{L}}_{\alpha }{\tilde{\eta }}_1{N}_1+Y_{\alpha 2}{\bar{L}}_{\alpha }{\tilde{\eta }}_2{N}_2+\frac{1}{2}{M}_1\overline{N_1^c}{N}_1+\frac{1}{2}{M}_2\overline{N_2^c}{N}_2+h.c,$$

$$V(H,{\eta }_1,{\eta }_2)=V_H+V_{{\eta }_1}+V_{{\eta }_2}+V_{\mathrm{int}},$$

$$V_H$$

$$V_{{\eta }_1}$$

$$V_{{\eta }_2}$$

$$V_{\mathrm{int}}$$

$${\eta }_1$$

$$H$$

$$m_h^2$$

$$m_{H_1}^2$$

$$m_{A_1}^2$$

$$m_{{\eta }_1^{+}}^2$$

$$m_{H_2}^2$$

$$m_{A_2}^2$$

$$m_{{\eta }_2^{+}}^2$$

$$(m_{H_1},m_{A_1},m_{{\eta }_1^{+}},m_{H_2},m_{A_2},m_{{\eta }_2^{+}},{\lambda }_{L1},{\lambda }_{L2},{\lambda }_{{\eta }_1},{\lambda }_{{\eta }_2},{\lambda }_3^{\prime },{\lambda }_4^{\prime },{\lambda }_5^{\prime }).$$

$${\mu }_1^2$$

$${\lambda }_3$$

$${\lambda }_4$$

$${\lambda }_5$$

$${\mu }_2^2$$

$$\widetilde{{\lambda }_3}$$

$$\widetilde{{\lambda }_4}$$

$$\widetilde{{\lambda }_5}$$

$${\Omega }_{\mathrm{DM}}{h}^2=\left[1.07\times 10^9{\mathrm{GeV}}^{-1}\right]\frac{x_f{g}_{\ast }^{1/2}}{g_{\ast s}{M}_{\mathrm{Pl}}⟨\sigma v{⟩}_f},$$

$$⟨\sigma v⟩=\frac{1}{8m_{\mathrm{DM}}^{4}T{K}_2^2\left(\frac{m_{\mathrm{DM}}}{T}\right)}\underset{4m_{\mathrm{DM}}^2}{\overset{\infty }{\int }}\sigma (s-4m_{\mathrm{DM}}^2)\sqrt{s}{K}_1\left(\frac{\sqrt{s}}{T}\right)ds,$$

$${\Omega }_{\mathrm{DM}}{h}^2=0.120\pm 0.001={\Omega }_1{h}^2+{\Omega }_2{h}^2,$$

$$\frac{d{y}_{H_1}}{dx}$$

$$\frac{d{y}_{H_2}}{dx}$$

$$\displaystyle Y_{i}^{EQ}(x)=0.145\frac{g}{g_{*}}x^{3/2}\bigg{(}\frac{m_{H_{i}}}{\mu}\bigg{)}^{3/2}e^{-x\big{(}\frac{m_{H_{i}}}{\mu}\big{)}}.$$

$${\Omega }_i{h}^2$$

$$\frac{d{y}_{H_1}}{dx}$$

$$\frac{d{y}_{H_2}}{dx}$$

$${\sigma }_i^{\text{SI}}=\frac{{\lambda }_L^2{f}_n^2}{4\pi }\frac{{\mu }_{i,n}^2{m}_n^2}{m_h^4{m}_i^2},$$

$${\sigma }_{\mathrm{i},\mathrm{eff}}^{\text{SI}}=\frac{{\Omega }_{H_i}}{{\Omega }_{\mathrm{total}}}{\sigma }_{H_i}^{\text{SI}}.$$

$$\frac{d\Phi }{dE}(△\Omega )=\frac{1}{8\pi }⟨\sigma v⟩\frac{J(△\Omega )}{M_{\text{DM}}^2}\frac{dN}{dE},$$

$$J(△\Omega )={\int }_{△\Omega }d{\Omega }^{\prime }{\int }_{\mathrm{LOS}}dl{\rho }^2(l,{\Omega }^{\prime }).$$

$$(m^{\nu }{)}_{\alpha \beta }$$

$$=$$

Full text

I Introduction

There have been irrefutable amount of evidences in favour of the existence of non-luminous, non-baryonic form of matter in the universe, popularly known as dark matter (DM). The presence of this form of matter has also been supported by astrophysical observations like the ones related to galaxy clusters by Fritz Zwicky Zwicky (1933) back in 1933, observations of galaxy rotation curves in 1970’s Rubin and Ford (1970), the more recent observation of the bullet cluster by Chandra observatory Clowe et al. (2006) along with several galaxy survey experiments which map the distribution of such matter based on their gravitational lensing effects. There is equally robust evidence from cosmology as well, suggesting that around $26\%$ of the present universe’s energy density is in the form of dark matter. In terms of density parameter $\Omega_{\rm DM}$ and $h=\text{Hubble Parameter}/(100\;\text{km}~{}\text{s}^{-1}\text{Mpc}^{-1})$, the present DM abundance is conventionally reported as Aghanim et al. (2018): $\Omega_{\text{DM}}h^{2}=0.120\pm 0.001$ at 68% CL.

In spite of these astrophysical and cosmology based evidences, there have been no detection of particle DM at any experiments. The direct detection experiments like LUX Akerib et al. (2017), PandaX-II Tan et al. (2016); Cui et al. (2017) and Xenon1T Aprile et al. (2017, 2018) have continued to produce null results so far. As the Standard Model (SM) of particle physics can not accommodate such a form of matter, several beyond the Standard Model (BSM) proposals have been put forward Taoso et al. (2008) out of which the weakly interacting massive particle (WIMP) paradigm is the most widely studied one. While such interactions are capable of explaining correct relic density of DM, the same interactions can also give rise to production of DM particles at the large hadron collider (LHC) Kahlhoefer (2017). However, nothing is found so far in these searches also, putting strict bounds on DM coupling to the SM particles, particularly quarks. Another detection prospect lies in the indirect detection frontier where searches are going on to find excess of antimatter, gamma rays or neutrinos, originating perhaps from dark matter annihilations (for stable DM) or decay (for long lived DM). While no convincing DM signal has been observed yet, there are tight constraints on DM annihilations into SM particles Ahnen et al. (2016), specially the charged ones which can finally lead to excess of gamma rays for WIMP type DM.

Though the null results mentioned above have not ruled out all the parameter space for a single particle DM models yet, it may be suggestive of presence of a much richer DM sector. The idea may be natural given the fact that the visible sector is made up of several generations for single type of particles. There have been several proposals for multi-component WIMP dark matter during last few years, some of which can be found in Cao et al. (2007); Zurek (2009); Chialva et al. (2013); Heeck and Zhang (2013); Biswas et al. (2013); Bhattacharya et al. (2013); Bian et al. (2014, 2015); Esch et al. (2014); Karam and Tamvakis (2015, 2016); DiFranzo and Mohlabeng (2017); Bhattacharya et al. (2017a); Dutta Banik et al. (2017); Klasen et al. (2017); Ghosh et al. (2018); Ahmed et al. (2017); Bhattacharya et al. (2017b); Ahmed et al. (2017); Borah and Dasgupta (2017); Bhattacharya et al. (2018, 2019); Aoki and Toma (2018); Dutta Banik et al. (2018); Barman et al. (2018a); Yaser Ayazi and Mohamadnejad (2019); Poulin and Godfrey (2018); Chakraborti and Poulose (2018); Chakraborti et al. (2018); Bernal et al. (2018); Elahi and Khatibi (2019); Borah et al. (2019). Such multi-component DM scenarios, even if both the DM candidates are of the same type, can have very interesting signatures at direct as well as indirect detection experiments, as studied in Profumo et al. (2009); Fukuoka et al. (2011); Cirelli and Cline (2010); Aoki et al. (2012, 2013); Modak et al. (2015); Geng et al. (2014a); Gu (2013); Aoki et al. (2014); Geng et al. (2014b, 2015); Biswas et al. (2015); Borah et al. (2015, 2016, 2017); Herrero-Garcia et al. (2017, 2018). Since direct and indirect detection (considering annihilations only, for stable DM) rates of DM are directly proportional to the DM density and DM density squared respectively and thereby producing tight constraints on single DM models, multi-component DM models can remain safe from being ruled out by such direct and indirect search constraints if the relative densities of different DM components are within appropriate limits. On top of that, such multi-component DM often comes with additional features like giving rise to interesting indirect detection signatures like monochromatic X-ray or gamma ray lines, as explored in several works, see for example Biswas et al. (2013); Borah and Dasgupta (2017); Borah et al. (2015, 2016, 2017); Dutta Banik et al. (2017) and references therein.

Multi-component DM scenarios may also be connected with other sector of particle physics. One such immediate possibility evolves through a probable connection with neutrino physics, particularly with the origin of neutrino mass and mixing. Results of several experiments in last two decades like T2K Abe et al. (2011), Double Chooz Abe et al. (2012), Daya Bay An et al. (2012), RENO Ahn et al. (2012) and MINOS Adamson et al. (2013) have confirmed the existence of non-zero but tiny neutrino mass and large (compared to quark mixing) leptonic mixing Fukuda et al. (2001); Ahmad et al. (2002a, b); Abe et al. (2008, 2011, 2012); An et al. (2012); Ahn et al. (2012); Adamson et al. (2013); Patrignani et al. (2016). Similar to the case of DM, these experimental observations provide clear indication for BSM physics as neutrino mass can not be explained within SM framework. Several BSM models attempt to explain tiny neutrino mass by incorporating additional fields. Apart from the conventional type I seesaw Minkowski (1977); Gell-Mann et al. (1979); Mohapatra and Senjanovic (1980); Schechter and Valle (1980), there exist other variants of seesaw mechanisms also, namely, type II seesaw Mohapatra and Senjanovic (1981); Lazarides et al. (1981); Wetterich (1981); Schechter and Valle (1982); Brahmachari and Mohapatra (1998), type III seesaw Foot et al. (1989) and so on.

It is particularly interesting to think of possible connection between the origin of neutrino mass and dark matterBhattacharya et al. (2016); Bhattacharya et al. (2017c), and perhaps such a connection is most straightforward in scotogenic scenarios, originally proposed by Ma Ma (2006a). In scotogenic type of model, the $\mathbb{Z}_{2}$ odd particles take part in radiative generation of light neutrino masses, while the lightest $\mathbb{Z}_{2}$ odd particle plays the role of DM. The salient feature of this framework is the common origin of light neutrino mass and DM where one can constrain the model from observables in both the sectors, and hence enhancing the predictive power of the model. While there are several possible implementation of this scotogenic framework (for a review, please refer to Cai et al. (2017)) in single component DM scenarios, there have not been much studies on the role of multi-component DM on the origin of neutrino mass. For a recent work on the role of two component DM on radiative origin of neutrino mass, one may refer to Aoki et al. (2017) 111In another recent work Bernal et al. (2018), two component fermion DM was proposed as a new anomaly free gauged $B-L$ model.. We adopt a similar setup here while sticking to the most minimal scenario in order to accommodate two component scalar doublet DM with correct total relic abundance, satisfying the neutrino oscillation data and other relevant constraints from dark matter direct, indirect detections and lepton flavour violation (LFV).

We consider a $\mathbb{Z}_{2}\times\mathbb{Z}^{\prime}_{2}$ extension of the Standard Model such that each $\mathbb{Z}_{2}$ sector consists of a singlet neutral fermion and a scalar field doublet under $SU(2)_{L}$. While the minimal scotogenic model consists of three singlet neutral fermion and one scalar doublet, our scenario contains two singlet fermions and two scalar doublets, keeping the particle content as minimal as the original model. Although there exists the possibility of singlet fermion DM also in either or both the $\mathbb{Z}_{2}$ sectors, we stick to scalar doublet DM. The reason behind this choice is three fold: (a) the gauge interactions of DM due to which the correct thermal abundance can be obtained easily without requiring large dimensionless couplings to enhance annihilation cross section222Fermion singlet DM in such models typically require large Yukawa couplings to satisfy correct relic and often run into the problems of vacuum stability Lindner et al. (2016)., (b) the gauge interactions of scalar doublet DM enhance its detection prospects at direct, indirect search experiments, (c) there exists an intermediate region for single component scalar doublet DM where relic abundance criteria can not be satisfied which makes it worth studying if two component scalar doublet DM can fill the void. Single component inert doublet DM (IDM) and its extensions have been studied by several authors in varieties of contexts Borah and Cline (2012); Dasgupta and Borah (2014); Cirelli et al. (2006); Barbieri et al. (2006); Ma (2006b); Lopez Honorez et al. (2007); Hambye et al. (2009); Dolle and Su (2009); Lopez Honorez and Yaguna (2010, 2011); Gustafsson et al. (2012); Goudelis et al. (2013); Arhrib et al. (2014); Plascencia (2015); D??az et al. (2016); Ahriche et al. (2018); Borah and Gupta (2017); Borah et al. (2018a, b); Barman et al. (2018b). We find that the total relic abundance of two component scalar doublet DM can be satisfied in the intermediate region while being consistent with neutrino oscillation data. The model also predicts the lightest neutrino mass to be zero. While the DM candidates satisfy the constraints from direct and indirect detection experiments, there lies the tantalising possibility to probe these scenarios at future searches in these frontiers and also in rare decay experiments like $\mu\rightarrow e\gamma$.

This paper is organised as follows. In section II, we discuss the model and the particle spectrum. In section III we discuss the details of two component dark matter pointing out the different annihilation channels contributing to the individual and total relic abundance, constraints from direct, indirect search followed by discussions on the constraints from neutrino oscillation data in section IV. We then briefly comment on lepton flavour violation in section V. We discuss our results in section VI and finally conclude in section VII.

II The Model

We have extended the particle content of the Standard Model by introducing two $SU(2)_{L}$ scalar doublets $\eta_{1}$ and $\eta_{2}$ and two right handed (RH) neutrinos, $N_{1,2}$. Furthermore, we include additional discrete symmetries, $\mathbb{Z}_{2}\times\mathbb{Z}^{\prime}_{2}$ under which all SM fields transform trivially. The charge assignments of these additional fields under SM gauge symmetry as well as additional global discrete symmetries are indicated in Table 1. The two neutral CP even scalars out of these extra doublets form the multi-component dark matter framework while the presence of RH neutrinos are instrumental in realising the light neutrino mass similar to that of scotogenic model Ma (2006a). The additional discrete symmetries not only explain the stability of individual DM components, but also prevent Dirac neutrino mass at tree level by forbidding the couplings involving lepton doublets ($L$), singlet neutral fermions and the SM Higgs $H$. The framework therefore serves as the minimal set-up in getting multi (two) component DM model which can accommodate light neutrino mass.

The Yukawa Lagrangian involving interactions between the new fields of our model and SM fields can be written as

[TABLE]

where $\alpha,\beta=e,\mu,\tau$ stand for different generations of SM leptons. Note that as each RH neutrinos are odd under two different $\mathbb{Z}_{2}$ sectors, the corresponding RH neutrino mass matrix remains diagonal.

The most general renormalisable scalar potential of our model, $V(H,\eta_{1},\eta_{2})$, consistent with $SU(2)_{L}\times U(1)_{Y}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{\prime}$ consists of (i) $V_{H}$: sole contribution of the SM Higgs , (ii) $V_{\eta_{1}}$: $\eta_{1}$ contribution, (iii) $V_{\eta_{2}}$: $\eta_{2}$ contribution and (iv) $V_{\rm{int}}$: interactions among $H,\eta_{1},\eta_{2}$. This can be written as follows.

[TABLE]

where

[TABLE]

and

[TABLE]

In order to keep $\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{\prime}$ unbroken so as to guarantee the stability of DM components, the neutral components of $\eta_{1}$ and $\eta_{2}$ are chosen not to acquire any non-zero vacuum expectation value (VEV) and hence they can be identified as two inert Higgs doublets (IHD). These IHDs can be parametrised as

[TABLE]

On the other hand, the neutral component of the SM Higgs field acquires a non-zero VEV (denoted by $v$) which is responsible for electroweak symmetry breaking (EWSB). We parametrise the Higgs field $H$ as

[TABLE]

As follows from the scalar field Lagrangian, after EWSB, the physical scalars have the following masses

[TABLE]

The scalar potential should be bounded from below in order to make the electroweak vacuum stable. This poses some constraints on the scalar couplings of the model. In addition to this, all the relevant couplings should also maintain the perturbativity. These bounds together with the unitarity limit have been studied for the three Higgs doublets and can be found in Chakrabarty (2016), from which we obtain the limits in our case.

We are effectively left with two-component inert DM scenario where the presence of light neutrinos are also taken cared of. As specified before, we consider the CP even neutral components $H_{1},H_{2}$ as two DM candidates without any loss of generality. Analogous to the single IHD scenario, we define $\lambda_{L1}\equiv\frac{\lambda_{3}+\lambda_{4}+\lambda_{5}}{2}$ and $\lambda_{L2}\equiv\frac{\tilde{\lambda_{3}}+\tilde{\lambda_{4}}+\tilde{\lambda_{5}}}{2}$, which denote the individual Higgs portal couplings of two DM candidates respectively. For our analysis purpose, we first implement this model in LanHEP Semenov (2016) choosing the independent parameters in the scalar sector as

[TABLE]

For simplicity, we will consider the couplings of quartic interaction between two IHDs ($\eta_{1}$ and $\eta_{2}$) to be the same and identify it by $\lambda_{12}$, $i.e.\lambda_{3}^{\prime},\lambda_{4}^{\prime},\lambda_{5}^{\prime}=\lambda_{12}$. We express other couplings of the scalar sector in terms of these independent parameters as follows,

[TABLE]

Note that the parameters $\lambda_{\eta_{1}}$ and $\lambda_{\eta_{2}}$ do not take part directly in the DM phenomenology.

III Dark Matter Phenomenology

The set-up contains two dark matter components: $H_{1}$ and $H_{2}$. In order to find the final relic density, their annihilations and co-annihilations are to be considered. In addition, the role of neutrino Yukawa couplings are also important. Below we provide a systematic approach to calculate the relic density in our scenario. Constraints from dark matter search are also mentioned.

III.1 Relic Density

For a single component WIMP type DM, the DM candidate with mass $m_{DM}$ is part of the thermal plasma in the early universe which eventually freezes out when the rate of annihilations fall below the rate of expansion of the universe. The final abundance can then be obtained by solving the Boltzmann equation for the DM number density. In fact, for DM annihilations dominated by $s$-wave processes only, the relic abundance can be approximated as  Kolb and Turner (1990)

[TABLE]

where $g_{*}$ and $g_{*s}$ are the effective relativistic degrees of freedom that contribute to the energy density and entropy density, respectively. $x_{f}$ is to be determined from the parameter $x=m_{DM}/T$ evaluated at the freeze-out temperature $T_{f}$ and is given by $x_{f}\equiv[\frac{m_{\rm DM}}{T}]_{T=T_{f}}\ =\ \ln\left(0.038\frac{g}{g_{*}^{1/2}}M_{\text{Pl}}m_{\rm DM}\langle\sigma v\rangle_{f}\right)$, with $g$ being the number of internal degrees of freedom of the DM. The thermally averaged annihilation cross section, given by Gondolo and Gelmini (1991)

[TABLE]

is evaluated at $T_{f}$ and denoted by $\langle\sigma v\rangle_{f}$. The freeze-out temperature $T_{f}$ is derived from the equality condition of DM interaction rate $\Gamma=n_{\rm DM}\langle\sigma v\rangle$ with the rate of expansion of the universe $H(T)\simeq\sqrt{\frac{\pi^{2}g_{*}}{90}}\frac{T^{2}}{M_{\rm Pl}}$. In the above expression of Eq.(10), $K_{i}(x)$’s are the modified Bessel functions of order $i$. As is well known, if the mass splitting within the DM multiplet is relatively small, there can be additional contributions from coannihilation channels Griest and Seckel (1991), whose importance in IDM has already been discussed in several earlier works.

In the present model, we have two dark matter candidates $H_{1}$ and $H_{2}$. Since both the candidates contribute to the dark matter relic density obtained from PlanckAghanim et al. (2018) experiment, one must satisfy the following relation:

[TABLE]

where $h$ denotes the reduced Hubble parameter and the relic density of the $H_{1}$ and $H_{2}$ is given by $\Omega_{1}h^{2}$ and $\Omega_{2}h^{2}$ respectively. Since there exists annihilation channels through which $H_{1}$ can go into $H_{2}$ (or vice versa, depending upon which one is heavier), the Boltzmann equations for the two DM candidates are coupled, in general. In such a case, there is no approximate formula which we can use for individual DM abundance like we had in case of single component DM discussed above.

Before going to the details of coupled Boltzmann equations, we first identify and categorise different annihilation channels of $H_{1},H_{2}$ as shown in Fig. 1, co-annihilations channels involving individual DM candidates in Fig. 2, along with conversion between DM candidates $H_{1},H_{2}$ (depending upon which one is heavier) in Fig. 3. Feynman diagrams including contribution from neutrino Yukawa interactions are included in Fig. 4.

First we write down the coupled Boltzmann equations without considering the effect of neutrino Yukawa coupling ($i.e.$ excluding Fig. 4). Due to the involvement of two components of DM in our case, we modify the definition of parameter $x$ from before as $x=\mu/T$, where $\mu$ is the reduced mass defined through:$~{}\mu=\frac{m_{H_{1}}m_{H_{2}}}{m_{H_{1}}+m_{H_{2}}}$. Therefore the coupled Boltzmann equations, written in terms of the common variable $x=\mu/T$ and the comoving number density $Y_{\rm DM}=n_{\rm DM}/s$ ($s$ being the entropy density), are obtained as333We adopt the notation from a recent article on two component DM Bhattacharya et al. (2019),

[TABLE]

Here $y_{i}$ ($i=H_{1,2}$) is related to $Y_{i}$ by $y_{i}=0.264M_{Pl}\sqrt{g_{*}}\mu Y_{i}$ and similarly for equilibrium density, $y_{i}^{EQ}=0.264M_{Pl}\sqrt{g_{*}}\mu Y_{i}^{EQ}$, where the equilibrium distributions are now recasted in terms of $\mu$ having the form

[TABLE]

Here $M_{\rm Pl}=1.22\times 10^{19}~{}{\rm GeV}$ and $g_{*}=106.7$ and $X$ represents SM particles. One should note that the contribution to the Boltzmann equations coming from the DM-DM conversion (corresponding to Fig. 3) will depend on the mass hierarchy of DM particles. This is described by the use of $\Theta$ function in the above equations. These coupled equations can be solved numerically to find the asymptotic abundance of the DM particles, $y_{i}\left(\frac{\mu}{m_{H_{i}}}x_{\infty}\right)$, which can be further used to calculate the relic:

[TABLE]

where $x_{\infty}$ indicates a very large value of $x$ after decoupling.

In presence of neutrino Yukawa couplings, the Boltzmann equations get modified and are given by

[TABLE]

The total relic density of DM follows from the combined contribution of both the components and is given by $\Omega_{\rm{total}}=\Omega_{H_{1}}h^{2}+\Omega_{H_{2}}h^{2}$.

III.2 Direct Detection

As mentioned earlier, DM parameter space can be constrained significantly by the null results at different direct detection experiments such as LUX Akerib et al. (2017), PandaX-II Tan et al. (2016); Cui et al. (2017) and Xenon1T Aprile et al. (2017, 2018). There are two ways scalar DM can scatter off nuclei at tree level in our model. One is elastic scattering mediated by SM Higgs boson while the other one is the inelastic one mediated by electroweak gauge bosons. The latter can be kinematically forbidden by considering large mass splitting between IHD components, typically larger than the average kinetic energy of DM particle. The spin independent elastic scattering cross section mediated by SM Higgs (shown in Fig. 5) is given as Barbieri et al. (2006)

[TABLE]

where $\mu_{i,n}=m_{n}m_{i}/(m_{n}+m_{i})$ is the ${\rm DM}$-nucleon reduced mass and $\lambda_{L_{i}}$ is the quartic coupling involved in ${\rm DM}$-Higgs interaction. The index $i$ stands for DM candidate in our scenario: $H_{1},~{}H_{2}$. A recent estimate of the Higgs-nucleon coupling $f$ gives $f=0.32$ Giedt et al. (2009) although the full range of allowed values is $f=0.26-0.63$ Mambrini (2011). In this two-component DM framework, the spin-independent cross section relevant for each of the candidate can be expressed as

[TABLE]

Latest results from Xenon-1T experiment provides a strong constraint on single component IDM as it restricts $\lambda_{L}$ coupling significantly. However due to the presence of two DM components here in our set-up, such tight constraints can be evaded by suitable adjustment of relative DM abundance. We will discuss the status of our model at direct detection frontier in subsequent sections.

III.3 Indirect Detection

As mentioned earlier, WIMP DM candidates have good prospects at indirect detection experiments looking for excess of gamma rays. DM particles can annihilate and produce SM particles, out of which photons (and also neutrinos), being electromagnetically neutral, have better chances of reaching the detector from source without getting deflected. Following the notations of Ahnen et al. (2016), the observed differential gamma ray flux produced due to the DM annihilation can be computed as

[TABLE]

where $\triangle\Omega$ is the solid angle corresponding to the observed region of the sky, $\langle\sigma v\rangle$ is the thermally averaged DM annihilation cross section, $dN/dE$ is the average gamma ray spectrum per annihilation process and the astrophysical $J$ factor is given by

[TABLE]

In the above expression, $\rho$ is the DM density and LOS corresponds to line of sight. Therefore, measuring the gamma ray flux and using the standard astrophysical inputs, one can constrain the DM annihilation first into different charged final states like $\mu^{+}\mu^{-},\tau^{+}\tau^{-},W^{+}W^{-},b\bar{b}$ which in turn produces the gamma rays. As discussed in case of direct detection, here also our set-up carries some flexibility as far as indirect detection constraints are concerned. We incorporate the global analysis of the Fermi-LAT and MAGIC observations of dSphs Ahnen et al. (2016) for this purpose. The bounds quoted in Ahnen et al. (2016) consider $100\%$ annihilation of DM into particular final states as well as assume a single DM component which fills the entire $26\%$ of the universe. Since our construction involves deviation from these consideration, we can make the bounds weaker by playing with the branching fraction to a particular final states and simultaneously changing the relative fractional abundance. This is because the DM annihilation rates are directly proportional to number density squared of DM in local neighbourhood.

IV Neutrino Mass

Similar to the minimal scotogenic model, here also light neutrino masses are generated at one loop level, as shown in Fig. 6. However, since each $\mathbb{Z}_{2}$ sector contains only one singlet fermion and one IHD, combination of one loop contribution from both $\mathbb{Z}_{2}\times\mathbb{Z^{\prime}}_{2}$ sectors leads to the light neutrino mass matrix elements as

[TABLE]

where $m_{H_{i}}$ and $m_{A_{i}}$ are provided in Eq. (7). Expressions of $\Lambda_{ii}$ is then found to be

[TABLE]

The light neutrino mass matrix $m_{\nu}$ can be diagonalised through

[TABLE]

where $m^{\rm diag}_{\nu}={\text{diag}}(m_{{\nu}_{1}},m_{{\nu}_{2}},m_{{\nu}_{3}})$. Considering the charged lepton mass matrix as a diagonal one, the light neutrino diagonalising matrix $U_{\text{PMNS}}$ coincides with the lepton mixing matrix. The $U_{\text{PMNS}}$ matrix can be parametrised as

[TABLE]

where $c_{ij}=\cos{\theta_{ij}},\;s_{ij}=\sin{\theta_{ij}}$ and $\delta$ is the leptonic Dirac CP phase. The diagonal phase matrix $U_{\text{P}}=\text{diag}(1,e^{i\alpha},e^{i\beta})$ contains the Majorana CP phases $\alpha,\beta$ which remain undetermined at neutrino oscillation experiments. We summarise the $3\sigma$ global fit values in Table 2 from the recent global fit Esteban et al. (2019), which we use in our subsequent analysis.

Note that the two $\mathbb{Z}_{2}$ sectors in our model can generate at most two light neutrino masses444With the involvement of a third RH neutrino to start with and making it very heavy, it effectively leads to the two RH neutrino scenario we consider here. As a consequence of this limit, it is shown Antusch et al. (2012) that one of the light neutrino remains massless., which needs to be consistent with light neutrino data mentioned above. This leaves us with two possibilities (a) $m_{1}=0,m_{2}<m_{3}$, (b) $m_{1}<m_{2},m_{3}=0$ corresponding to normal and inverted hierarchies respectively. Since the inputs from neutrino data are only in terms of the mass squared differences and mixing angles (phases are considered to be zero here), it would be useful for our purpose to express the Yukawa couplings in terms of light neutrino parameters. This is possible through the Casas-Ibarra (CI) parametrisation Casas and Ibarra (2001) extended to radiative seesaw model Toma and Vicente (2014) which allows us to write the Yukawa couplings as

[TABLE]

where $\Lambda$ is the 2$\times$2 diagonal matrix with eigenvalues defined in Eq.(22) and $R$ is a complex orthogonal matrix Ibarra and Ross (2003, 2004); Antusch et al. (2012) having the form,

[TABLE]

Using Eq.(25), the elements of $3\times 2$ Yukawa matrix can be obtained with specific choices of the complex angle $z$. The same calculation can be repeated for inverted hierarchy as well. In the subsequent sections, we discuss how the constraints from neutrino sector can play a non-trivial role in the dark matter parameter space from relic abundance criteria.

V Lepton Flavour Violation

Since the charged lepton flavour violating (LFV) decays are very much suppressed in the SM, any observation of such effects will be a clear signature of beyond the SM physics. In our model, due to the coupling of each $\mathbb{Z}_{2}$ sector particles ($N_{i}$ and $\eta_{i}$) to the SM leptons, one may expect some contribution to such LFV effects at one loop level. The same fields that take part in the one-loop generation of light neutrino mass as shown in Fig. 6 also mediate LFV processes like $\mu\rightarrow e\gamma$. The neutral scalar in the internal lines of Fig. 6 will be replaced by their charged counterparts (which emit a photon) whereas the external fermion legs can be replaced by $\mu,e$ respectively, generating the one-loop contribution to $\mu\rightarrow e\gamma$.

As the couplings and masses involved in this process are the same as the ones that generate light neutrino masses and play a role in DM relic abundance, we can no longer choose them arbitrarily. It should be noted that each $\mathbb{Z}_{2}$ sector contributes separately to this process and hence we have to add the respective contributions at amplitude level. Adopting the notations from Baek et al. (2014), we can write

[TABLE]

Here $\alpha_{\rm em}=e^{2}/4\pi$ is the electromagnetic fine structure constant, $G_{F}$ is the Fermi constant, $F(x)$ = $(1-6x+3x^{2}+2x^{3}-6x^{2}\log{x})/(6(1-x)^{4})$ is the loop factor and $x=\frac{M_{i}^{2}}{m_{\eta_{i}^{+}}^{2}}$. The MEG experiment provides the most stringent upper limit on the branching ratio: ${\rm Br}(\mu\rightarrow e\gamma)<4.2\times 10^{-13}$ Baldini et al. (2016).

VI Results and Discussion

As mentioned earlier, we write the model in LanHEP Semenov (2016) and then extract the model files to use in micrOMEGAs 4.1 Belanger et al. (2015) for two component DM framework. To have a better understanding of the role of neutrino Yukawa coupling on the phenomenology of DM, we first discuss the results (i) without any involvement of neutrino Yukawa interactions and then (ii) incorporate the Yukawa interactions. Finally constraints from light neutrino data and LFV decays are studied in order to constrain the parameter space.

VI.1 Results in the absence of neutrino Yukawa coupling

As mentioned earlier in the introduction, one of the key motivations for studying the two component scalar doublet dark matter over other existing possibilities is related to the existence of an intermediate region of a single scalar doublet dark matter mass $80\;\text{GeV}\leq M_{\rm DM}\leq 500\;\text{GeV}$ where correct relic abundance can not be satisfied irrespective of the choices of other parameters. While in the single component DM framework, this intermediate region can be revived by additional production mechanisms Molinaro et al. (2014); Borah and Gupta (2017), there still exists severe bounds from indirect detection experiments which effectively rule out single scalar doublet DM mass below around 400 GeV Borah and Gupta (2017). Reviving this intermediate region by incorporating an additional DM component in the form of a fermion and vector boson have been studied recently by the authors of Chakraborti and Poulose (2018) and Chakraborti et al. (2018) respectively. In this work, we consider the other remaining possibility of introducing a second scalar doublet DM component which carries non-trivial implications for light neutrino masses (which most of the multi-component DM models do not address) also. Since intermediate mass region is special in a sense that one component scalar doublet DM is not enough to produce correct relic density, we stick to that mass range (roughly) for our numerical analysis.

The left panel of Fig. 7 shows the variation of individual contribution to the relics $\Omega_{1(2)}h^{2}$ against the mass of dark matter $m_{H_{1(2)}}$ such that $\Omega_{\rm total}h^{2}$ satisfies the Planck 2018 limit. The same variation for different values of the conversion coupling $\lambda_{12}$ is shown with three different choices of $\lambda_{12}$ = 0 (green for $H_{1}$ and orange for $H_{2}$ contribution), 0.1 (black for $H_{1}$ and purple for $H_{2}$ contribution), 0.5 (brown for $H_{1}$ and red for $H_{2}$ contribution). Furthermore, we consider $m_{H_{2}}>m_{H_{1}}$ within our range of interest of DM mass. The other relevant parameters are kept fixed at their respective benchmark values: mass splitting within $\eta_{1}(\eta_{2}),~{}i.e.~{}\Delta M_{i=1,2}=m_{A_{i}}-m_{H_{i}}=1$ GeV and $\lambda_{L_{1,2}}=0.01$. This mass splitting plays role in individual coannihilation diagrams, similar to single component IDM. Such choice of $\Delta M_{i}$ is motivated from the fact that the maximum contribution from a single DM component toward relic density can obtained for such small mass difference. $\lambda_{L_{1,2}}$ values are chosen for representative purpose.

It is intriguing to note also from Fig. 7 (both from left and right panel) that a change of pattern of the patches of different colours happened around $m_{H_{1}}=m_{H_{2}}\simeq 380$ GeV and that small region (where all these points meet) is almost independent of the change of $\lambda_{12}$ value. At the meeting point, where both the DM masses are around 380 GeV, almost $50\%$ individual contribution follows from both the DM candidates. Since $H_{1,2}$ are degenerate at this point, no number changing process can take place at this point irrespective of the fact whether $\lambda_{12}$ is zero or non-zero. Hence it explains why all the patches meet over this region. Note that while obtaining these plots, we keep the IHD masses within the stipulated intermediate range: 80-500 GeV. Finally what we observe from the plots is that the mass range 230-500 GeV becomes allowed range while incorporating a second inert doublet.

Fig. 7 is clearly an indicative of the importance of the parameter $\lambda_{12}$, which plays a non-trivial role in conversion of one DM candidate into the other. In order to understand this, let us begin with $\lambda_{12}~{}=~{}0$ case, as shown with green and orange patches. When $\Omega_{1}h^{2}$ is small, $\Omega_{2}h^{2}$ should be large such that $\Omega_{\rm total}h^{2}$ can satisfy the Planck 2018 limit. This indicates that a point (say $m_{H_{1}}$=320 GeV with $\Omega_{1}h^{2}$ = 0.0402 ) in the lower side of green patch ($i.e.$ with low $m_{H_{1}}$) is actually correlated to a single point (say $m_{H_{2}}$= 445 GeV with $\Omega_{2}h^{2}$ = 0.0789) near the higher side orange patch (with relatively large $m_{H_{2}}$). Effect of the conversion coupling $\lambda_{12}$ would be clear with the plot in right panel of Fig. 7, where the total relic satisfied points are placed in $m_{H_{1}},m_{H_{2}}$ plane. Given a fixed $m_{H_{1}}$, required masses of $m_{H_{2}}$ are shown (hence the code of $H_{2}$ are given as used in left panel) with $m_{H_{2}}>m_{H_{1}}$. Note that for $m_{H_{1}}=320$ GeV, the relic constraint would be satisfied by $m_{H_{2}}$ = 491 GeV with $\lambda_{12}=0.1$. Hence it is clear that due to conversion $H_{2}H_{2}\rightarrow H_{1}H_{1}$, somewhat higher value (as compared to that of $\lambda_{12}$ = 0) of $m_{H_{2}}$ becomes allowed. One should keep in mind that here the combination of couplings $\lambda^{\prime}_{3}$, $\lambda^{\prime}_{4}$, $\lambda^{\prime}_{5}$ appears in the conversion processes via $\lambda_{12}$ and hence even for the smaller values of individual couplings, large conversion effects can be seen.

Also note that once we are away from this specific point, a typical pattern (parabolic behaviour) of the plots sets in. In order to have a better understanding of it, we plot the variation of individual and total relic abundance with DM masses for fixed values of conversion parameter $\lambda_{12}=0.5$ and $m_{H_{2}}$ = 450 GeV, while keeping other parameters to their benchmark values. The results are shown in Fig. 8. For $m_{H_{1}}$ below 450 GeV, conversion like $H_{2}H_{2}\rightarrow H_{1}H_{1}$ is effective. As the mass of $H_{1}$ crosses the threshold of $H_{2}$ mass (here 450 GeV), the conversion ($H_{1}H_{1}\rightarrow H_{2}H_{2}$) becomes efficient thereby reducing the abundance of $H_{1}$ while increasing the latter’s. The total relic curve (green line) is almost symmetric around this cross-over point which explains the typical nature of the left panel plot of Fig. 7 in the high mass regime with large conversion coupling.

We also check the effects of considering different mass splittings $\Delta M_{1}\neq\Delta M_{2}$. As we increase $\Delta M_{2}$ to 10 GeV while keeping $\Delta M_{1}$ at 1 GeV, the relic abundance of $H_{2}$ decreases, a feature observed in single component IDM also in the high mass regime. The results are shown in Fig. 9 (considering $m_{H_{2}}>m_{H_{1}}$) which shows that the masses of $H_{2}$ are shifted to higher mass regions in order to satisfy the total DM relic abundance. This is expected as we know based on our knowledge of single component inert doublet DM analysis that making the mass splitting related to one IHD more, the corresponding relic density would be less. This clarifies why there is a separation between the green and orange patches (with $\lambda_{12}=0$). Similar situation prevails when nonzero $\lambda_{12}$ value is switched on as well, which can be seen from black (brown) and purple (red) patches. Also notice that the existence of the symmetric point similar to the case with equal $\Delta M_{i}$ (about which the relic contour lines gain a typical shape with non-zero $\lambda_{12}$) is lost here. This is related to the fact that now two components of DM have different type of co-annihilations and hence contribute to the total relic differently. This is also reflected from a mixed up distribution pattern of correct relic satisfied points put in $m_{H_{1}},m_{H_{2}}$ plane in the right panel of Fig. 9. One more point to notice is that there is a shift toward the higher mass range of DM as compared to the case displayed in Fig. 7 with $m_{H_{2}}>m_{H_{1}}$.

As mentioned earlier, scalar doublet DM has good prospects for detection at direct search experiments. This is due to the presence of gauge as well as Higgs portal interactions with nucleons. While the gauge portal interactions at tree level can be forbidden kinematically, the Higgs portal interactions can still give rise to sizeable spin independent elastic scattering of DM off nucleons. For representative purpose of our model, we fix the Higgs portal coupling $\lambda_{L_{1,2}}$ at 0.01 for both the DM candidates and provide their direct detection cross section for different conversion parameters $\lambda_{12}$ in the left panel of Fig. 10. Note that we include those points only which satisfy the total DM relic abundance. We have taken into account the relative abundance of individual DM candidates while calculating the direct detection cross section. The colour code used here is consistent to the discussion above related to Fig. 7. As can be seen from this plot, the scattering cross sections in our model remains well below the latest direct search bounds. However, with increase in the value of Higgs portal coupling, these points will be closer to the experimental upper bound and hence providing the model a good prospect of being detected at ongoing and near future experiments like LZ Akerib et al. (2015), XENONnT Aprile et al. (2016), DARWIN Aalbers et al. (2016) and PandaX-30T Liu et al. (2017).

One should note that it will be difficult to distinguish the direct detection signals coming from both the DM candidates in reality. As pointed out in Herrero-Garcia et al. (2017) it is possible to distinguish the direct detection signals for a two-component scenario from that of the single-component dark matter provided the two dark matter masses differ significantly (by an order of magnitude or more). In our present set-up, both the dark matter masses are with in the intermediate mass range ($80-500$) GeV and hence the possiblity is no more useful here.

We also check the prospects for indirect detection of DM in our model by specifically focussing on $W^{+}W^{-}$ final states from DM annihilations. This is due to the chosen mass range of DM where annihilation to this final states is the most dominant one. We incorporate the relative abundance of the two DM candidates while calculating their annihilation rates. We include the factor coming from the branching fraction to this final state from DM annihilations. Our findings are displayed in the right panel of Fig. 10. It can be seen that for some part of the parameter space, the indirect detection cross section lies close to the experimental upper bound. This plot also points out the increase in allowed parameter space for two component inert doublet DM compared to one component inert DM where DM mass below 400 GeV was ruled out by Fermi-LAT constraints Borah and Gupta (2017).

VI.2 Results in the presence of neutrino Yukawa couplings

After discussing the new results for two component inert DM with salient features arising due to gauge and Higgs portal interactions along with four point interactions, we now move on to discuss the role of neutrino Yukawa interactions with two-component DM candidates. This has non-trivial connections to the light neutrino sector as the same Yukawa couplings also play a role in generating light neutrino masses. The new annihilation channels that will come into play now are the ones shown in Fig. 4. We now have an additional free parameter which is important for dark matter phenomenology defined as $\Delta M_{NiHi}~{}=~{}M_{i}-m_{H_{i}}$. For analysis purpose, we fix it to $\Delta M_{NiHi}~{}=~{}10~{}{\rm GeV}$ value. Such small mass splitting between neutral singlet fermion and DM candidates enhances the co-annihilation cross section effectively having significant effect on the relic abundance as we will see soon.

In obtaining the relic density and direct/indirect detection cross sections in our model, we vary the IHD masses within the so called intermediate range of DM mass. Since we have chosen $\Delta M_{NiHi}~{}=~{}10~{}{\rm GeV}$, this can automatically fix the two RH neutrino masses for chosen values of $m_{H_{1},H_{2}}$. Hence with the benchmark values (as specified in last subsection) of $\Delta M_{1,2}$, we would get the elements of $\Lambda$ matrix. Then using Eq. (25), one can determine the neutrino Yukawa couplings for some choices of $z$ involved in $R$ matrix. We have made the simplest choice by setting $z$ = 0 in order to determine $R$ matrix. The light neutrino masses are evaluated from the best fit values of solar and atmosphere mass-square splittings as one of the light neutrino remains massless in the scenario.

We first reproduce the parameter space which gives rise to the total DM relic abundance in agreement with Planck 2018 limit and also allowed by direct detection bounds for different benchmark values of conversion coupling $\lambda_{12}$. The individual contributions of $H_{1}$ and $H_{2}$ to the total relic are displayed in Fig. 11(a). The same results, in terms of the relic contours in $m_{H_{1}}$-$m_{H_{2}}$ plane (while maintaining $m_{H_{2}}>m_{H_{1}}$), are shown in Fig. 11(b) which if compared with the one in the absence of Yukawa interactions (see Fig. 7(b)) reveals interesting shift in parameter space. For example with $\Delta M_{1,2}$ = 1 GeV, the symmetric point (near $m_{H_{1}}\sim 380$ GeV) about which patterns of relic contours corresponding to different $\lambda_{12}$ values evolve in Fig. 7 now shifted to 330 GeV. As this point corresponds to $m_{H_{1}}=m_{H_{2}}$ having no contribution followed from co-annihilations involving DM components, new co-annihilations channels involving neutrino Yukawa couplings (see Fig. 4) would anyway be present. This causes the observed shift.

It turns out that the required Yukawa couplings are of order $10^{-4}-10^{-5}$ for both normal and inverted hierarchies. This estimate is obtained from $m_{H_{1}}-m_{H_{2}}$ contours (allowed by both the relic and direct detection constratints)555We have checked that with the inverted hierarchy, the results in terms of the values of parameters involved, do not change significantly. along with the expression in Eq.(25) with the use of Eq.(22). The other interesting observation is: incorporating Yukawa interactions allow smaller values of DM masses that satisfy correct total DM abundance and hence the entire intermediate region for single inert doublet DM becomes allowed with two-component inert DM model involving neutrinos. As stated above, the presence of the Yukawa couplings would introduce additional co-annihilation channels (see Fig. 4). While Fig. 4(a) and (b) contribute to both the DMs in a similar way, the DM-DM conversions (via Fig. 4(c) and (d)) can effectively alter the two DM components differently (depending on their mass hierarchy) as seen in the Boltzmann equation via Eq.(15a) and Eq.(15b). So with $m_{H_{2}}>m_{H_{1}}$, in general the processes $H_{2}N_{1}\rightarrow H_{1}N_{2}$ (with non-zero $\lambda_{12}$) and $H_{2}N_{2}\rightarrow H_{1}N_{1}$ would enhance the number density of $H_{1}$ component 666With the choice of $\Delta M_{N_{i}H_{i}}=10$ GeV only $H_{2}N_{2}\rightarrow H_{1}N_{1}$ would contribute. and hence its contribution towards relic ($\Omega_{1}h^{2}$) becomes more while compared with a similar $m_{H_{1}}$ having no Yukawa interactions. This can be seen by comparing Fig.7(a) and Fig. 11(a). For $m_{H_{1}}=250$ GeV, $\Omega_{1}h^{2}$ was 0.026 as seen from Fig. 7(a). Now for same $m_{H_{1}}$, but with Yukawa interaction, $\Omega_{1}h^{2}$ becomes 0.033. Accordingly $m_{H_{2}}$ is shifted toward a lower mass value to satisfy the total relic to be consistent with observation. In this particular example, $m_{H_{2}}$ comes down to 380 GeV (Fig. 11(a)) from its earlier value of 492 GeV (Fig. 7(a)). Although the Yukawa interactions are not supposed to change the direct and indirect detection rates at tree level, due to the change in the prefactor $\frac{\Omega_{H_{i}}}{\Omega_{total}}$ of $\sigma^{SI}_{H_{i}}$, some changes in $\sigma^{SI}_{i,eff}$ are expected. This is shown in Fig. 12.

We finally show the predictions for LFV decay, $\mu\rightarrow e\gamma$, in Fig. 13 by choosing similar benchmark values of mass splittings as in the dark matter analysis (used in Fig. 11) against variation with $m_{H_{1}}$. In obtaining this plot, we have used that pair of $m_{H_{1}},m_{H_{2}}$ values which satisfy relic density constraint with their total contribution. The lower set of plots is for $\Delta M_{i}=1$ GeV while upper set is with $\Delta M_{i}=1$ MeV. Note that we have carried out the DM phenomenology with $\Delta M_{i}=1$ GeV only. If we work with further small values of mass-splitting, the DM phenomenology would not change, particularly when we concentrate our interest within the intermediate range of DM mass: 80 GeV - 500 GeV. However while making both of these splittings smaller, it will increase branching ratio related to the prediction of LFV decays as seen from the Fig. 13. This is particularly due to the increase in Yukawa couplings for smaller values of mass splitting, in accordance with the CI parametrisation. Hence with not-so-small values of the mass splittings, we can conclude that the range of the parameter space allowed by the DM relic density, direct and indirect search results and neutrino data is mostly unaffected by the LFV constraints.

Finally it is perhaps pertinent to discuss the collider propects of the present scenario. Due to the involvement of the associated charged components of the two IHDs, the typical opposite sign di-lepton plus missing energy signature can be quite interesting at LHC Bhattacharya et al. (2019). On top of that, in our present model, both the DM masses are within the intermediate range ($<500$ GeV). Hence it would be more promising in our present set-up in terms of detection at LHC compared to the IHD DM present in the usual scotogenic scenarioMa (2006a). We provide the production cross-section for $pp\rightarrow H^{+}H^{-}$ in the following Table 3 which was generated in CalcHEP Belyaev et al. (2012) at $\sqrt{s}=14$ TeV using cteq6l1 Placakyte et al. (2011) parton distribution corresponding to a benchmark scenario with $m_{H_{1}(H_{2})}=300(350)~{}\rm{GeV}$. The respective relic density and the direct detection cross-section are also indicated in the same table.

VII Conclusion

While there is no compelling reason to believe the dark matter sector to be composed of only one type of particle, it still remains as a well motivated scenario particularly from minimality point of view. However, in certain minimal models, there can be motivations for a multi-component dark matter framework originating from inadequacies of a single component to fit for dark matter relic density as well as connections to other interesting observed phenomena, like neutrino mass. Motivated by this, we have studied a minimal extension of the Standard Model by two right handed neutrinos and two scalar doublets having non-trivial transformations under an unbroken $\mathbb{Z}_{2}\times\mathbb{Z}^{\prime}_{2}$ symmetry. While light neutrino masses arise at one loop level in a way similar to the scotogenic models consisting of particles from both the discrete $\mathbb{Z}_{2}$ sectors (playing the role of loop mediators), this minimal choice of the additional particles result in vanishing lightest neutrino mass.

The two dark matter candidates are chosen to be the neutral components of the two scalar doublets which are thermally produced in the early universe by virtue of their electroweak gauge interactions. We particularly focus on the mass range in between $W$ boson mass and approximately 500 GeV where single component scalar doublet dark matter can not satisfy correct relic abundance. While scalar doublet relic does not depend crucially upon the Yukawa couplings to leptons, turning the Yukawa interactions on has interesting consequences for relic due to additional coannihilation channels. The inter-conversion between two dark matter candidates also play instrumental role in deciding the total relic abundance. We show that the model can satisfy correct total relic abundance criteria in the intermediate mass regime in agreement with the latest data from Planck mission. In addition to this, the model remains very much predictive at ongoing direct, indirect detection experiments as well as rare decay experiments looking for charged lepton flavour violation. Unlike typical multi-component dark matter models, our model gets more restrictive due to electroweak gauge interactions as well as the non-trivial roles both the dark matter particles play in generating light neutrino masses, opening up detection prospects through lepton portals like charged lepton flavour violation. While the origin of the discrete symmetries $\mathbb{Z}_{2}\times\mathbb{Z}^{\prime}_{2}$ remains unexplained in the current work, we leave a more detailed study looking for UV completion of such models to an upcoming work.

Acknowledgements.

DB acknowledges the support from IIT Guwahati start-up grant (reference number: xPHYSUGI-ITG01152xxDB001), Early Career Research Award from DST-SERB, Government of India (reference number: ECR/2017/001873) and Associateship Programme of IUCAA, Pune. RR thank Biswajit Karmakar, Abhijit Kumar Saha and Amit Dutta Banik for various useful discusions while carrying out this work.