A minimal model for two-component FIMP dark matter: A basic search

Weakly interacting massive particles (WIMPs) are the most popular solution to the puzzle of dark matter (DM) [1–3]. In TeV scale (LHC scale) new physics, DM particles follow the thermal scenario in which they reach thermal and chemical equilibrium with the bath particles but lose it at the freeze-out temperature (which is around m_DM/20) and experience decoupling from the Universe plasma. WIMP candidates such as the neutralino [4] and Kaluza-Klein particle [5, 6] are found in theories such as the minimal supersymmetric Standard Model (MSSM) and universal extra dimensions (UED), respectively, and also in other extensions of the Standard Model (SM) such as singlet scalar [7–10] (or fermionic [11–13]) DM. In spite of their popularity, WIMPs have not yet been detected in direct experiments.

The other viable and well-motivated hypothesis to explain the DM problem is that there is such a feeble interaction that DM particles can never be abundant enough to thermalize. In this so-called freeze-in mechanism [14–16], feebly interacting massive particles (FIMPs) have been slowly produced in the early Universe through the collisions or decays of the bath particles. FIMP candidates are motivated in various extensions of the SM [17–19] and a well-known example which arises from neutrino physics is the sterile neutrino [20–24]. It is difficult to detect FIMP particles because of their small couplings with the SM. As for the indirect searches, depending on the type of DM candidate, i.e. scalar [25], fermion [19], etc, some experiments have parameter space where they could survive but these are very borderline. For a study of the non-thermal properties of dark matter see Ref. [26].

Although a lot of attention has been dedicated to single-particle DM models, some studies have considered DM models with the contribution of more particles in the observed DM density (multi-component DM [27–31]). The simplest and the most common case is the union of the singlet scalar (fermionic) and the singlet fermionic (scalar) models which are employed in both freeze-in [25, 32] and freeze-out [33–36] solutions (or intermediate cases [37]). Nevertheless, it remains a mystery whether DM is a single particle or multi-component.

In this paper, we analyze whether the freeze-in approach can be properly used to produce the observed DM density in our Universe. We choose a minimal two-component DM model in such a way that both of the components are FIMP particles. Following our hypothesis, we consider a singlet scalar and a Dirac fermion where an accidental symmetry guarantees their stabilities and a Higgs portal enables them to interact with the SM particles. The most striking feature of our model is its simplicity, as the two candidates of DM particles do not couple with each other and the model has separate overlaps with both the singlet scalar model [17] and the singlet fermionic model [18]. In our work, all contributing processes to the relic density are assumed and supplementary phenomenological aspects are also included.

Some promising possible signatures of FIMPs which are found to be most reliable in previous works are the γ-ray excess observed from the Galactic center (GC) [32, 37], the X-ray line at 3.55 keV [32, 37], and DM self-interaction [25, 37]. To generate the gamma ray excess, the fermionic component should have a pseudoscalar coupling to the mediator in the freeze-out regime [37]. The scalar component which does also couple to the SM Higgs directly [32], should not feature large valued couplings. An X-ray signal with E_γ = 3.55 keV from the XMM-Newton telescope and a similar signal at 3.52 keV from the Andromeda galaxy (M31) and Perseus Cluster could all be interpreted by the decay [38] or the annihilation [39] of DM. However, this requires a definite decay rate and annihilation cross section which is out of reach for our scenario. Therefore, we continue our probe relying only on the DM self-interaction. This non-gravitational interaction is a well-motivated indirect search as it solves the tensions between observations and simulations of the small-scale structure of DM.

Following the aforementioned setup, our paper is organized as follows. After introducing the construction of our model and identifying its parameter space in Section 2, we solve two independent Boltzmann equations in the following section (Section 3), in order to reach the observed relic density measured by the WMAP and Planck experiments [40]. In Section 4, we study the phenomenological implications for both direct and indirect experiments, and summarize our results in Section 5.

Beyond the SM, we employ three new fields to furnish our model: two scalars (χ and S) and one Dirac fermion (ψ), which are all assumed to be singlet under the SM gauge groups. A discrete Z₂ symmetry is applied such that it reads the SM fields and the S-scalar even, and the other two fields (χ and ψ) odd. This symmetry guarantees the stability of both odd particles in a way that we do not have any terms involving both fields ψ and χ. In this way, the decays of odd particles to one another are prevented. Therefore, we can have two DM candidates in our setup by an accidental symmetry. The framework of our model is constructed by:

$$\begin{eqnarray}&&\begin{array}{l} {\mathcal L} \supset \frac{1}{2}{({\partial }_{\mu }S)}^{2}+\frac{1}{2}{({\partial }_{\mu }\chi )}^{2}+i\bar{\psi }\rlap{/}{\partial }\psi \\ \,-{m}_{\psi }\bar{\psi }\psi -{g}_{s}S\bar{\psi }\psi -{g}_{p}S\bar{\psi }{\gamma }^{5}\psi \\ \,-V(H,S,\chi ),\end{array}\end{eqnarray} \tag{ 1 }$$

where we introduced the scalar and pseudoscalar interactions with the couplings g_s and g_p, respectively, and inserted the scalar interactions in the term V(H,S,χ) as

$$\begin{eqnarray}\begin{array}{ll}V(H,S,\chi )= & -{\mu }_{H}^{2}{H}^{\dagger }H+{\lambda }_{H}{({H}^{\dagger }H)}^{2}+{\mu }_{1}S+\frac{1}{2}{\mu }_{S}^{2}{S}^{2}\\ & +\frac{1}{3!}{\alpha }_{s}{S}^{3}+\frac{1}{4!}{\lambda }_{S}{S}^{4}+\frac{1}{2}{m}_{0\chi }^{2}{\chi }^{2}+\frac{1}{4!}{\lambda }_{\chi }{\chi }^{4}\\ & +{\lambda }_{1}S{H}^{\dagger }H+{\lambda }_{2}{S}^{2}{H}^{\dagger }H+{\lambda }_{\chi H}{\chi }^{2}{H}^{\dagger }H\\ & +{\lambda }_{3}S{\chi }^{2}+{\lambda }_{4}{S}^{2}{\chi }^{2}.\end{array}\end{eqnarray} \tag{ 2 }$$

After spontaneous symmetry breaking, the SU(2) Higgs doublet is parametrized as

$$\begin{eqnarray}&&H=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ {v}_{H}+h\end{array}\right),\end{eqnarray} \tag{ 3 }$$

where v_H = 246 GeV, but for the mediator we assume that it does not acquire a vacuum expectation value, i.e. < S > = 0, which minimalizes our model too. Now, due to the interaction terms in Eq. (2), h and S mix with each other and form a mass matrix with the following eigenstates

$$\begin{eqnarray}\begin{array}{ll}{h}_{1} & =S\sin \theta +h\cos \theta,\\ {h}_{2} & =S\cos \theta -h\sin \theta,\end{array}\end{eqnarray} \tag{ 4 }$$

and the eigenvalues as

$$\begin{eqnarray}&&{m}_{{h}_{1},{h}_{2}}^{2}=\frac{{m}_{h}^{2}+{m}_{S}^{2}}{2}\pm \frac{{m}_{h}^{2}-{m}_{S}^{2}}{2}\sqrt{1+{y}^{2}},\quad {\rm{with}}\quad y=\frac{2{m}_{h,S}^{2}}{{m}_{h}^{2}-{m}_{S}^{2}},\end{eqnarray} \tag{ 5 }$$

where θ is the mixing angle between h₁ and h₂ such that

$$\begin{eqnarray}&&\tan \theta =\frac{y}{1+\sqrt{1+{y}^{2}}}.\end{eqnarray} \tag{ 6 }$$

According to the definition of the mixing angle θ, h₁ can be considered as the SM-like Higgs observed at the LHC with a mass of about 125 GeV. In Eq. (5), ${m}_{h}=\sqrt{2{\lambda }_{H}}{v}_{H}$, ${m}_{S}={({\lambda }_{2}{v}_{H}^{2}+{\mu }_{S}^{2})}^{1/2}$ and ${m}_{h,S}=\sqrt{{\lambda }_{1}{v}_{H}}$.

Concerning our parameters, vacuum stability implies that the scalar potential in Eq. (2) must be bounded from below. On the other hand, perturbativity does not allow the model parameters to be too large. Eventually, these theoretical conditions can be satisfied if one has

$$\begin{eqnarray}&&-2\pi /3\lt {\lambda }_{S},{\lambda }_{\chi }\lt 2\pi /3,\\ &&\,-4\pi \lt {\lambda }_{2},{\lambda }_{\chi H}\lt 4\pi,\\ &&-8\pi \lt {\lambda }_{4},{g}_{s},{g}_{p}\lt 8\pi,\\ &&\,{\lambda }_{2}+\sqrt{{\lambda }_{H}{\lambda }_{S}}\gt 0,\\ &&{\lambda }_{\chi H}+\sqrt{{\lambda }_{H}{\lambda }_{\chi }}\gt 0,\\ &&\,2{\lambda }_{4}+\sqrt{{\lambda }_{S}{\lambda }_{\chi }}\gt 0,\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}\begin{array}{ll} & \left(\sqrt{2({\lambda }_{2}+\sqrt{{\lambda }_{H}{\lambda }_{S}})({\lambda }_{\chi H}+\sqrt{{\lambda }_{H}{\lambda }_{\chi }})(2{\lambda }_{4}+\sqrt{{\lambda }_{S}{\lambda }_{\chi }})}\right.\\ + & \sqrt{{\lambda }_{H}{\lambda }_{S}{\lambda }_{\chi }}+{\lambda }_{2}\sqrt{{\lambda }_{\chi }}+{\lambda }_{\chi H}\sqrt{{\lambda }_{S}}+2{\lambda }_{4}\sqrt{{\lambda }_{H}})\gt 0,\end{array}\end{eqnarray} \tag{ 8 }$$

where λ_H is the quartic coupling of H. Extending the SM with the new fields ψ, χ and S embeds 19 parameters in addition to the SM ones. They are m_S, m_χ, m_ψ, g_s, g_p, μ₁, α_S, λ_S, λ_χ, λ₁, λ₂, λ_{χ H}, λ₃, λ₄, m_h1, m_h2, sinθ, v_H, and m_h. However, due to the 8 model constraints, 11 independent parameters,

$$\begin{eqnarray}&&{g}_{s},\,{g}_{p},\,{\lambda }_{3},\,{\lambda }_{4},\,{\lambda }_{\chi H},\,{m}_{\psi },\,{m}_{\chi },\,{m}_{{h}_{2}},\,\sin \theta,\end{eqnarray} \tag{ 9 }$$

remain for the relic abundance and for indirect searches, λ_H and λ_χ are required. Here, we take a moment to describe the eight constraints which appear in our work. Note that, after spontaneous symmetry breaking the scalar potential given in Eq. (2) reads as V(h,S,χ). Therefore, it can be deduced from the potential that:

$$\begin{eqnarray}&&1)\quad {\left.\frac{\partial V}{\partial h}\right|}_{h=S=\chi =0}=0\Rightarrow {\mu }_{H}^{2}={\lambda }_{H}{v}_{H}^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&2)\quad {\left.\frac{\partial V}{\partial S}\right|}_{h=S=\chi =0}=0\Rightarrow {\mu }_{1}=-\frac{{\lambda }_{1}{v}_{H}^{2}}{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&3)\quad {m}_{h}^{2}=-{\mu }_{H}^{2}+3{\lambda }_{H}{v}_{H}^{2}=2{\lambda }_{H}{v}_{H}^{2},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&4)\quad {m}_{S}^{2}={\mu }_{S}^{2}+{\lambda }_{2}{v}_{H}^{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&5)\quad {m}_{\chi }^{2}={m}_{0\chi }^{2}+{\lambda }_{\chi H}{v}_{H}^{2}.\end{eqnarray} \tag{ 14 }$$

Also, the mixing between S and h produces the scalars h₁ and h₂ so one can conclude that

$$\begin{eqnarray}&&\frac{1}{2}{m}_{h}^{2}{h}^{2}+\frac{1}{2}{m}_{S}^{2}{S}^{2}+{\lambda }_{1}{v}_{H}Sh=\frac{1}{2}{m}_{{h}_{1}}^{2}{h}_{1}^{2}+\frac{1}{2}{m}_{{h}_{2}}^{2}{h}_{2}^{2}.\end{eqnarray} \tag{ 15 }$$

Substituting h₁ = Ssinθ + hcosθ and h₂ = Scosθ − hsinθ (4) and using the constraints (3)–(5), one obtains:

$$\begin{eqnarray}&&6)\quad {\lambda }_{H}=\frac{{m}_{{h}_{1}}^{2}{\cos }^{2}\theta +{m}_{{h}_{2}}^{2}{\sin }^{2}\theta }{2{v}_{H}^{2}},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&7)\quad {\lambda }_{2}=\frac{{m}_{{h}_{1}}^{2}{\sin }^{2}\theta +{m}_{{h}_{2}}^{2}{\cos }^{2}\theta -{\mu }_{S}^{2}}{{v}_{H}^{2}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&8)\quad {\lambda }_{1}=\frac{{m}_{{h}_{1}}^{2}-{m}_{{h}_{2}}^{2}}{2{v}_{H}}\sin 2\theta .\end{eqnarray} \tag{ 18 }$$

These 8 constraints reduce the 19 free parameters in the model to the 11 independent parameters. Also, the couplings α_s and λ₂ can be taken as zero without any ambiguities. However, we consider α_s ≠ 0 for future applications. In the following, we will probe our model parameter space with experimental constraints coming from the relic density, direct and indirect detections.

3.1. Boltzmann Equation

Since our model contains two DM candidates, its relic density has contributions from both fields ψ and χ. Therefore, we have to solve two Boltzmann equations for particles which will not reach equilibrium in the freeze-in mechanism where we follow the solution in Ref. [19] (following Ref. [15]). The time evolution of number density, dn_DM/dt, for the fermionic DM is given by

$$\begin{eqnarray}&&\frac{{\rm{d}}{n}_{\psi }}{{\rm{d}}t}+3H{n}_{\psi }=\frac{T}{{\pi }^{2}}\displaystyle \sum _{i=1}^{2}{m}_{{h}_{i}}^{2}{K}_{1}\left(\frac{{m}_{{h}_{i}}}{T}\right){\varGamma }_{{h}_{i}\to \bar{\psi }\psi }+\frac{T}{32{\pi }^{4}}\displaystyle \sum _{j=f,Z,W,{h}_{1},{h}_{2}}\displaystyle {\int }_{4{m}_{j}^{2}}^{\infty }{\rm{d}}s{\sigma }_{jj\to \bar{\psi }\psi }(s)(s-4{m}_{j}^{2})\sqrt{s}{K}_{1}\left(\frac{\sqrt{s}}{T}\right),\end{eqnarray} \tag{ 19 }$$

and for the scalar DM, it reads

$$\begin{eqnarray}&&\frac{{\rm{d}}{n}_{\chi }}{{\rm{d}}t}+3H{n}_{\chi }=\frac{T}{{\pi }^{2}}\displaystyle \sum _{i=1}^{2}{m}_{{h}_{i}}^{2}{K}_{1}\left(\frac{{m}_{{h}_{i}}}{T}\right){\varGamma }_{{h}_{i}\to \chi \chi }+\frac{T}{32{\pi }^{4}}\displaystyle \sum _{j=f,Z,W,{h}_{1},{h}_{2}}\displaystyle {\int }_{4{m}_{j}^{2}}^{\infty }{\rm{d}}s{\sigma }_{jj\to \chi \chi }(s)(s-4{m}_{j}^{2})\sqrt{s}{K}_{1}\left(\frac{\sqrt{s}}{T}\right).\end{eqnarray} \tag{ 20 }$$

Here H is the Hubble constant, K₁ is the modified Bessel function of order 1 and s is the center of mass energy squared. All contributions to the DM relic density are considered in the corresponding cross sections and decay widths in the two above equations. Our analytical results for the cross sections and decay widths are presented in the Appendix. The number density of DM particles is calculated as ${n}_{i}=\frac{{g}_{i}}{{(2\pi )}^{3}}\displaystyle \int {{\rm{d}}}^{3}p{f}_{i}$ [16, 41, 42] (with i = ψ, χ), where f_i is the phase space density of particle i with the g_i-internal spin degrees of freedom. As it is well-known from the freeze-in mechanism of production, the two DM candidates in the present model have negligible initial abundance (individually), thus we may set f_i = 0. Consequently, it can be derived from Eqs. (19) and (20) that the process of DM conversion, i.e. $\chi \chi \leftrightarrow \bar{\psi }\psi$, does not contribute to the total relic abundance and is suppressed in our next calculations. On the other hand, each of the DM candidates, independent of the other, can be produced or annihilated in the Universe.

By solving Eqs. (19) and (20), one can obtain the number densities (n_ψ, n_χ) scaled to the entropy of Universe $\hat{s}$, i.e. ${Y}_{\psi }={n}_{\psi }/\hat{s}$ and ${Y}_{\chi }={n}_{\chi }/\hat{s}$, as

$$\begin{eqnarray}\begin{array}{ll}{Y}_{\psi }= & \frac{1}{4{\pi }^{4}}\frac{45{M}_{pl}}{1.66{g}_{s}^{* }(T)\sqrt{{g}_{\rho }^{* }}}\left[2\displaystyle \sum _{i=1}^{2}{\varGamma }_{{h}_{i}\to \bar{\psi }\psi }{m}_{{h}_{i}}^{2}\displaystyle {\int }_{{T}_{{\rm{Now}}}}^{\infty }{\rm{d}}T\frac{{K}_{1}\left(\frac{{m}_{{h}_{i}}}{T}\right)}{{T}^{5}}\right. \\ & \left.+\displaystyle \sum _{j=f,Z,W,{h}_{1},{h}_{2}}\frac{1}{16{\pi }^{2}}\displaystyle {\int }_{{T}_{{\rm{Now}}}}^{\infty }{\rm{d}}T\frac{1}{{T}^{5}}\displaystyle {\int }_{4{m}_{j}^{2}}^{\infty }{\rm{d}}s{\sigma }_{jj\to \bar{\psi }\psi }(s)(s-4{m}_{j}^{2})\sqrt{s}{K}_{1}\left(\frac{\sqrt{s}}{T}\right)\right],\end{array}\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}\begin{array}{ll}{Y}_{\chi }= & \frac{1}{4{\pi }^{4}}\frac{45{M}_{pl}}{1.66{g}_{s}^{* }(T)\sqrt{{g}_{\rho }^{* }}}\left[2\displaystyle \sum _{i=1}^{2}{\varGamma }_{{h}_{i}\to \chi \chi }{m}_{{h}_{i}}^{2}\displaystyle {\int }_{{T}_{{\rm{Now}}}}^{\infty }{\rm{d}}T\frac{{K}_{1}\left(\frac{{m}_{{h}_{i}}}{T}\right)}{{T}^{5}}\right. \\ & \left.+\displaystyle \sum _{j=f,Z,W,{h}_{1},{h}_{2}}\frac{1}{16{\pi }^{2}}\displaystyle {\int }_{{T}_{{\rm{Now}}}}^{\infty }{\rm{d}}T\frac{1}{{T}^{5}}\displaystyle {\int }_{4{m}_{j}^{2}}^{\infty }{\rm{d}}s{\sigma }_{jj\to \chi \chi }(s)(s-4{m}_{j}^{2})\sqrt{s}{K}_{1}\left(\frac{\sqrt{s}}{T}\right)\right],\end{array}\end{eqnarray} \tag{ 22 }$$

where M_pl is the Planck mass, and ${g}_{s}^{* }$ and ${g}_{\rho }^{* }$ are the effective numbers of degrees of freedom.

3.2. Relic Abundance

The most important constraint which should be satisfied in models describing DM is the observed relic density. As the Planck experiments have measured the current amount of DM [40], our first experimental constraint is described as

$$\begin{eqnarray}&&{\varOmega }_{{\rm{DM}}}{h}^{2}={\varOmega }_{\psi }{h}^{2}+{\varOmega }_{\chi }{h}^{2}=0.1199\pm 0.0027,\end{eqnarray} \tag{ 23 }$$

where h is the Hubble parameter scaled in units of 100 km/s.Mpc. Using the yield calculated in the previous section (Eqs. (21) and (22)), we can obtain the relic density as

$$\begin{eqnarray}&&{\varOmega }_{i}{h}^{2}=2.742\times {10}^{-8}(\frac{{m}_{i}}{{\rm{GeV}}}){Y}_{i}({T}_{0}),\quad i=\psi,\chi .\end{eqnarray} \tag{ 24 }$$

First, we start with the scalar component of the model. The dependence of DM density is evaluated over the relevant parameters. The predicted relic density of our model is best behaved at m_h2 = 100 GeV and sinθ = 0.01, with the required value of 5 × 10⁻¹⁰ GeV for mediator-scalar DM coupling λ₃. Two other couplings, λ_{χ H} and λ₄, are found to have major roles in controlling the relic density. By varying the singlet scalar DM mass (inspired by Eq. (24)), we probe our parameter space in two classes: different values of λ_{χ H} and of λ₄ (see Fig. 1).

Zoom In Zoom Out Reset image size

The behavior of the relic density Ω_χh² regarding different Higgs-scalar DM couplings, λ_{χ H}, is depicted in the logarithmic scale in Fig. 1(a). In addition to the fixed values of relevant parameters λ₃, m_h2 and sinθ, we have adopted λ_{χ H} as 10⁻¹², 5 × 10⁻¹³ and 10⁻¹³ as we fixed λ₄ = 10⁻¹³. Through Fig. 1(a), it is obvious that the resonance occurs at m_χ ∼ m_h2/2. For masses below the resonance, the relic density of the scalar component increases linearly in the log-scale, but for larger values (m_χ > m_h2/2) it seems that the relic density is independent of the mass m_χ. The difference between these two regions is due to the process h₂ → χχ which is allowed in the region below the resonance.

The complementary analysis of the scalar component is plotted in Fig. 1(b), where we have chosen λ₄ = 10⁻¹², 5 × 10⁻¹³ and 10⁻¹³. Regarding the resonance at m_χ ∼ m_h2/2, as in Fig. 1(a), it can be seen that for the region below the resonance, the relic density grows when the scalar mass m_χ increases. This part of the graph seems to be independent of the λ₄-value and it is enhanced by the h₂ → χχ process. After a significant drop at m_χ ∼ m_h2/2, the relic density seems to be independent of DM mass for the region above the resonance. It is mainly influenced by changing the quartic coupling λ₄. Finishing our investigation of the scalar component, it should be noted that this analysis has a good overlap with a singlet scalar model [17]. We continue our investigation in parallel by turning our attention to the fermionic DM in the logarithmic scale, too. As before, we consider two classes of variations defined by the effect of scalar (parameterized by g_s) and pseudoscalar (parameterized by g_p) interactions of ψ (see Fig. 2). We first look at the coupling g_s so its best effects are formed for the values of 10⁻⁸, 10⁻⁹ and 10⁻¹⁰ (Fig. 2(a)). Similar to the scalar case, the resonance position occurs at m_χ ∼ m_h2/2, so below this value one can observe the linear behavior of relic density which arises through the process h₂ → χχ. For m_χ > m_h2/2, the relic density changes by several order of magnitudes in a small interval of mass range. In this region, Ω_χh² decreases when m_χ increases. Note that the relic density is approximately independent of DM mass for large values of g_s. Also, for small enough values of g_s, decreasing g_s does not significantly change the relic density.

Zoom In Zoom Out Reset image size

We keep on probing our model parameter space by choosing the appropriate values of g_p as 10⁻⁹, 10⁻¹⁰ and 5 × 10⁻¹¹. In Fig. 2(b), we show the behavior of the relic density of the fermionic DM in terms of its mass. Here, there is a distinct point which should be expressed. As is seen, for very small values of g_p, the relic density is approximately independent of mass for massive DM. Similar to the scalar component, we can compare the fermionic component with models describing singlet fermionic DM like in Ref. [18].

4.1. Direct searches

In this section, we search for signals inspired from XENON100 [43] and LUX [44] in spin-independent elastic scattering of DM off nuclei. Our intended process includes the fundamental interaction of DM-quark which occurs via the t-channel mediated by scalars h₁ and h₂. Taking into account the contribution of each DM component and using fractions ${\xi }_{\psi }=\frac{{\varOmega }_{\psi }}{{\varOmega }_{{\rm{DM}}}}$ and ${\xi }_{\chi }=\frac{{\varOmega }_{\chi }}{{\varOmega }_{{\rm{DM}}}}$, we investigate whether the model parameter space could be affected by the experimental results in this way. To this end, we calculate the following cross sections,

$$\begin{eqnarray}&&{\sigma }_{SI}^{\psi }={\xi }_{\psi }\frac{{g}_{s}^{2}{\mu }_{m}^{2}{\sin }^{2}\theta {\cos }^{2}\theta }{\pi }{\left(\frac{1}{{m}_{{h}_{1}}^{2}}-\frac{1}{{m}_{{h}_{2}}^{2}}\right)}^{2}{\lambda }_{N}^{2},\end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray}\begin{array}{ll}{\rm{\ }}{\sigma }_{SI}^{\chi }= & {\xi }_{\chi }\frac{{\mu }_{m}^{2}}{4\pi {m}_{{h}_{1}}^{4}{m}_{{h}_{2}}^{4}{m}_{\chi }^{2}}\left[2{v}_{H}{\lambda }_{\chi H}({m}_{{h}_{1}}^{2}{\sin }^{2}\theta +{m}_{{h}_{2}}^{2}{\cos }^{2}\theta ),\right.\\ & {\left.+{\lambda }_{3}\sin 2\theta ({m}_{{h}_{1}}^{2}-{m}_{{h}_{2}}^{2})\right]}^{2}{\lambda }_{N}^{2},\end{array}\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&{\lambda }_{N}=\frac{{m}_{N}}{{v}_{H}}\left[\displaystyle \sum _{q=u,d,s}{f}_{q}+\frac{2}{9}\left(1-\displaystyle \sum _{q=u,d,s}{f}_{q}\right)\right]\approx 1.4\times {10}^{-3}.\end{eqnarray} \tag{ 27 }$$

Here, the parameter ${\mu }_{m}=\frac{{m}_{N}{m}_{{\rm{DM}}}}{{m}_{N}+{m}_{{\rm{DM}}}}$ is the reduced mass of DM-nucleon and ${m}_{\chi }={({m}_{0\chi }^{2}+{\lambda }_{\chi H}{v}_{H}^{2})}^{1/2}$ is the physical mass of the scalar DM. A cancellation effect [33] could occur when the two terms in Eq. (26) cancel each other out, giving a suppressed cross section which is not appropriate for our consideration. Generally, as was mentioned earlier, a necessary condition for our DM candidates to be nonthermal is that they have extremely small couplings (g_s, λ_{χ H}) which would yield cross sections out of the sensitivity of the aforementioned experiments by their established values of order ∼10⁻⁸−10⁻¹². Searching for other viable experiments, we consider the scattering of DM off free electrons in materials such as superconductors, semiconductors and graphene. From Ref. [45–47] it is seen that, although, these electron detectors are useful for light DM particles (${\mathscr{O}}$(MeV)), the mediator mass should also be of order ${\mathscr{O}}$(MeV), which is in conflict with the current model including m_h2 = 100 GeV. For this reason, we are not able to probe such FIMP models directly. This outcome is consistent with the lack of direct experimental signals to date.

4.2. Indirect searches

4.2.1. Invisible Higgs decays

Since the ATLAS and CMS have recorded the signature of the SM Higgs [48, 49], new searches have been prepared for DM phenomenology. This is done by considering the branching ratios of the Higgs, especially for decaying into light DM candidates,

$$\begin{eqnarray}&&{\rm{Br}}({h}_{1}\to {\rm{Invisible}})=\frac{{\varGamma }^{{\rm{inv}}}({h}_{1}\to \chi \chi )+{\varGamma }^{{\rm{inv}}}({h}_{1}\to \psi \psi )}{{\varGamma }^{{\rm{SM}}}+{\varGamma }^{{\rm{inv}}}}.\end{eqnarray} \tag{ 28 }$$

Regarding the experimental upper bound 0.23 for Br(h₁ → Invisible) [50], we see that the decays of Higgs to both DM ψ and χ are suppressed due to small couplings g_s and g_p in the former case, and small λ₄ and λ_{χ H} in the latter. However, another constraint comes from the decay of our Higgs to h₂ (if kinematically possible, i.e. m_h2 < m_h1/2) whose decay rate could be calculated as

$$\begin{eqnarray}&&\varGamma ({h}_{1}\to {h}_{2}{h}_{2})=\frac{{e}_{3}^{2}}{8\pi {m}_{{h}_{1}}}{\left(1-\frac{4{m}_{{h}_{2}}^{2}}{{m}_{{h}_{1}}^{2}}\right)}^{1/2}\Theta ({m}_{{h}_{1}}-2{m}_{{h}_{2}}),\end{eqnarray} \tag{ 29 }$$

where e₃ is the relevant vertex factor which is presented in the Appendix, see Eq. (A11). From Eq. (29), it can be seen that for sinθ = 0.01 and m_h1 = 125 GeV, the result is sensitive to the choice of α_s. Consequently, we investigate the behavior of the aforementioned decay rate regarding the mass of the Higgs h₁ and the relevant coupling α_s. The parameter space of our Higgs sector is plotted in Figs. 3(a) and 3(b) where, respectively, we have calculated Br(h₁ → h₂h₂) as a function of m_h2 and depicted the parameter space for the (α_s, m_h2)-plane which is consistent with experimental measurements.

Zoom In Zoom Out Reset image size

4.2.2. DM Self-Interactions

Of the different DM models, the collision-less cold DM (CDM) paradigm has been successful in explaining the large scale structure of the Universe. However, there are discrepancies between the CDM predictions and observations on smaller scales. The self-interacting DM (SIDM) paradigm has the potential to solve these issues (for a review of SIDM, see Ref. [51]). Although such interactions cannot be detected in experiments, we can infer bounds on σ_DM/m by evaluating the trajectory of DM in colliding galaxy clusters [52, 53]. An updated work [54] has considered a set of twelve galaxies and six clusters in order to cover different scales. Including the core sizes from dwarf to cluster (varying from 0.5 to 50 kpc), the aforementioned cross section is parametrized as

$$\begin{eqnarray}&&{\sigma }_{{\rm{DM}}}/m{\unicode{x007E}}0.1-2\,{{\rm{cm}}}^{2}{{\rm{g}}}^{-1}.\end{eqnarray} \tag{ 30 }$$

In this section, we analyze this constraint to see if it can put new limits on the parameter space of our model.

The DM self-interaction in the present model includes the processes χχ → χχ, ψψ → ψψ, χχ → ψψ, ψψ → χχ and χψ → χψ. Except for χχ → χχ and, ψψ → ψψ, the processes contain cross sections proportional to the coupling λ₄, which is very small in our work. Therefore, the specified processes do not contribute to this cosmological constraint. Concerning the processes, χχ → χχ and, ψψ → ψψ, we start first with the scalar component which has been studied in a singlet FIMP scalar model in Ref. [55]. Here, we just consider the contact interaction which is parameterized by the coupling λ_χ. Practically, we neglect the contributions from the s-channel mediated diagrams. This is due to the small couplings of the scalar DM with both the SM Higgs and the mediator and also due to the large masses which appear in the propagator. One way to vitalize the s-channel contribution might be through fine tuning by considering the scattering near resonance (similar to Ref. [55]). In this way, in the denominator of the propagator, m_h2 should be tuned such that |m_χ − m_h2| ≪ 1 GeV. Considering the values of couplings needed for the observed relic density, this scenario fails too. Therefore, following Ref. [55], we obtain the self-interaction cross section per mass m_χ, as

$$\begin{eqnarray}&&\frac{{\sigma }_{\chi }}{{m}_{\chi }}=\frac{9{\lambda }_{\chi }^{2}}{2\pi {m}_{\chi }^{3}},\end{eqnarray} \tag{ 31 }$$

where λ_χ is the quartic self-coupling of the scalar DM (see Eq. (2)). Following the theoretical constraints in Eqs. (7) and (8), we obtain an experimental upper bound of about 0.1 GeV on the mass of the scalar DM, which is depicted in Fig. 4. Going back to Figs. 1(a) and 1(b), we observe that this range of scalar mass can produce proper total relic density along with the contribution of the fermionic component.

Zoom In Zoom Out Reset image size

Another significant point, which we would like to clarify in this work, is the self-interaction of a singlet fermionic FIMP DM. In general, we have two concerns. First, it should be noted that significant self-scattering at dwarf scales requires the mediator masses to be smaller than 100 MeV [56]. In fact, following Ref. [56], the fermionic DM should satisfy the relation (m_ψ/10 GeV) (m_h2/100 MeV)² ∼ 1, with Yukawa interactions of strengths 10⁻⁵ to 1, which is in contradiction with our fermionic DM coupling and mediator mass. Second, if the mediator couples to the SM through a Higgs portal, one should make sure that the mediator decays before the start of Big Bang nucleosynthesis (BBN), so the decay products do not affect the BBN. Eventually, we require a mediator with a lifetime ∼ 1 s. One way to alleviate the second constraint in DM models with extremely weak interactions is to open a new decay channel for the mediator so it can decay faster. This is done in Refs. [57] and [58] by coupling the mediator to a light sterile neutrino (it should be noted that this new coupling does not affect the relic density). However, the first constraint (light mediator) is in conflict with our mediator of mass 100 GeV (and other usual two-component models). A promising solution seems to be to work at the resonance region to minimize this mass constraint [59]. Due to the small coupling g_p (pseudoscalar interaction type), and the fact that there is an energy (velocity) dependent correction to the width in the resonance region (as explained in Ref. [59]), we conclude that the resonance DM self-interaction scenario does not work in fermionic FIMP models.

We have constructed a minimal two-component model to analyze the implications of multi-component DM in the Universe. Using the freeze-in mechanism, we calculated the relic abundance predicted by our model and compared it to the observed relic abundance of DM. We started our investigation by proposing two DM particles: a real scalar and a Dirac fermion. Furthermore, a scalar mediator between the dark sector and the SM sector was added. The couplings for this interaction are assumed to be small as we are utilizing the freeze-in mechanism. We solved two independent Boltzmann equations in order to obtain the observed relic density with the contributions of both DM components. It should be noted here that at the time of finishing this work, a new version of micrOMEGAs [60] was presented which can compute the relic abundance of FIMP candidates. In the following, using theoretical constraints, we probed the model parameter space and compared our results with the relevant singlet models. Although it is difficult to probe FIMP particles, we looked for astrophysical probes, first considering direct detection. We considered the scattering of DM particles off nuclei and free electrons. As we explained, it is impossible to see this direct signature for our FIMP model.

In order to constrain the parameter space of our model, we also checked the limits from the invisible decay width of the Higgs. Finally, we probed the self-interaction of DM in this model. We used the bounds on non-gravitational interactions of DM in giant cluster collisions and constrained the mass of DM candidates in our model. In addition to the mentioned probes of DM, we can refer to the Big Bang nucleosynthesis and cosmic microwave background constraints regarding dark photons and dark Higgs [61, 62]. These neutral bosons mix with the SM photon (kinetically) and the SM Higgs, respectively, by accepting significant bounds on their coupling parameters.

Two-component DM is a starting point for considering multi-component configurations where DM consists of various types of fundamental particles (scalar, fermion, vector and etc). The freeze-in framework is also a well-motivated approach which may be probed more extensively by future experiments.

We are particularly grateful to Yonit Hochberg for giving us insights into the direct probes, and we would like to thank Takashi Toma, Ian Shoemaker, Madhurima Pandey and Anirban Biswas for useful discussions.

Here, we present our calculation of the fermionic DM production cross-sections which contribute to the relic density of our model:

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{\bar{f}f\to \bar{\psi }\psi }=\frac{s{\sin }^{2}2\theta }{32\pi }{N}_{c}{\left(\frac{{m}_{f}}{{v}_{H}}\right)}^{2}\left(1-\frac{4{m}_{f}^{2}}{s}\right)F,\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{ZZ\to \bar{\psi }\psi }=\frac{{\sin }^{2}2\theta }{36\pi }{\left(\frac{{m}_{Z}^{2}}{{v}_{H}}\right)}^{2}\left(2+\frac{{(s-2{m}_{Z}^{2})}^{2}}{4{m}_{Z}^{4}}\right)F,\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{{W}^{+}{W}^{-}\to \bar{\psi }\psi }=\frac{{\sin }^{2}2\theta }{36\pi }{\left(\frac{{m}_{W}^{2}}{{v}_{H}}\right)}^{2}\left(2+\frac{{(s-2{m}_{W}^{2})}^{2}}{4{m}_{W}^{4}}\right)F,\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray}&&F={\left(1-\frac{4{m}_{\psi }^{2}}{s}\right)}^{3/2}\left({g}_{\psi }^{2}+\frac{s{g}_{p}^{2}}{s-4{m}_{\psi }^{2}}\right)\left\{\frac{1}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{1}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}-\frac{2(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+2{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right\}.\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{1}{h}_{1}\to \bar{\psi }\psi }= & \frac{1}{4\pi s}{\left(1-\frac{4{m}_{\psi }^{2}}{s}\right)}^{1/2}\left\{({g}_{p}^{2}s+{g}_{\psi }^{2}(s-4{m}_{\psi }^{2}))\left[\frac{{e}_{1}^{2}{\sin }^{2}\theta }{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{2}^{2}{\cos }^{2}\theta }{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right.\right.\\ & \left.+\frac{{e}_{1}{e}_{2}\sin 2\theta ((s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}})}{({(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2})({(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2})}\right]-8{m}_{\psi }{g}_{\psi }{\sin }^{2}\theta \left[\frac{{e}_{1}\sin \theta (s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{2}\cos \theta (s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \times \left[\frac{\arctan h({x}_{1})[{g}_{p}^{2}(2{m}_{{h}_{1}}^{2}-3s)-{g}_{\psi }^{2}(2{m}_{{h}_{1}}^{2}-8{m}_{\psi }^{2}+s)]}{{x}_{1}(s-2{m}_{{h}_{1}})}+{g}_{p}^{2}-{g}_{\psi }^{2}\right]\\ & -4{\sin }^{4}\theta \left[{g}_{p}^{4}+{g}_{\psi }^{4}-2{g}_{p}^{2}{g}_{\psi }^{2}+\frac{{({m}_{{h}_{1}}^{2}({g}_{p}^{2}-{g}_{\psi }^{2})+4{g}_{\psi }^{2}{m}_{\psi }^{2})}^{2}-8{g}_{p}^{2}{g}_{\psi }^{2}{m}_{\psi }^{2}s}{{m}_{{h}_{1}}^{4}+{m}_{\psi }^{2}(s-4{m}_{{h}_{1}}^{2})}\right.\\ & -\frac{\arctan h({x}_{1})}{{x}_{1}{(s-2{m}_{{h}_{1}})}^{2}}\left({g}_{p}^{4}[6{m}_{{h}_{1}}^{4}+s(s-4{m}_{{h}_{1}}^{2})]-2{g}_{p}^{2}{g}_{\psi }^{2}[2{m}_{{h}_{1}}^{2}(3{m}_{{h}_{1}}^{2}-4{m}_{\psi }^{2})+s(s-4{m}_{{h}_{1}}^{2})]\right.\\ & \left.\left.\left.+{g}_{\psi }^{4}[2{m}_{{h}_{1}}^{2}(3{m}_{{h}_{1}}^{2}-8{m}_{\psi }^{2})+s(s-4{m}_{{h}_{1}}^{2})+16{m}_{\psi }^{2}(s-2{m}_{\psi }^{2})]\right)\right]\right\},\end{array}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{2}{h}_{2}\to \bar{\psi }\psi }= & \frac{1}{4\pi s}{\left(1-\frac{4{m}_{\psi }^{2}}{s}\right)}^{1/2}\times \left\{({g}_{p}^{2}s+{g}_{\psi }^{2}(s-4{m}_{\psi }^{2}))\left[\frac{{e}_{3}^{2}{\sin }^{2}\theta }{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}\right.\right.\\ & \left.+\frac{{e}_{4}^{2}{\cos }^{2}\theta }{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}+\frac{{e}_{3}{e}_{4}\sin 2\theta [(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}]}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right]-8{m}_{\psi }{g}_{\psi }{\cos }^{2}\theta \left[\frac{{e}_{3}\sin \theta (s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{4}\cos \theta (s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \times \left[\frac{\arctan h({x}_{2})[{g}_{p}^{2}(2{m}_{{h}_{2}}^{2}-3s)-{g}_{\psi }^{2}(2{m}_{{h}_{2}}^{2}-8{m}_{\psi }^{2}+s)]}{{x}_{2}(s-2{m}_{{h}_{2}})}+{g}_{p}^{2}-{g}_{\psi }^{2}\right]\\ & -4{\cos }^{4}\theta \left[{g}_{p}^{4}+{g}_{\psi }^{4}-2{g}_{p}^{2}{g}_{\psi }^{2}+\frac{{({m}_{{h}_{2}}^{2}({g}_{p}^{2}-{g}_{\psi }^{2})+4{g}_{\psi }^{2}{m}_{\psi }^{2})}^{2}-8{g}_{p}^{2}{g}_{\psi }^{2}{m}_{\psi }^{2}s}{2{m}_{{h}_{2}}^{4}+2{m}_{\psi }^{2}(s-4{m}_{{h}_{2}}^{2})}\right.\\ & -\frac{\arctan h({x}_{2})}{{x}_{2}{(s-2{m}_{{h}_{2}})}^{2}}\left({g}_{p}^{4}\left(6{m}_{{h}_{2}}^{4}+s(s-4{m}_{{h}_{2}}^{2})\right)-2{g}_{p}^{2}{g}_{\psi }^{2}\left(2{m}_{{h}_{2}}^{2}(3{m}_{{h}_{2}}^{2}-4{m}_{\psi }^{2})+s(s-4{m}_{{h}_{2}}^{2})\right)\right.\\ & \left.\left.\left.+{g}_{\psi }^{4}\left(2{m}_{{h}_{2}}^{2}(3{m}_{{h}_{2}}^{2}-8{m}_{\psi }^{2})+s(s-4{m}_{{h}_{2}}^{2})+16{m}_{\psi }^{2}(s-2{m}_{\psi }^{2})\right)\right)\right]\right\},\end{array}\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{1}{h}_{2}\to \bar{\psi }\psi }= & \frac{1}{4\pi s}{\left(1-\frac{4{m}_{\psi }^{2}}{s}\right)}^{1/2}\left\{[{g}_{p}^{2}s+{g}_{\psi }^{2}(s-4{m}_{\psi }^{2})]\left[\frac{{e}_{2}^{2}{\sin }^{2}\theta }{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{3}^{2}{\cos }^{2}\theta }{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right.\right.\\ & \left.+\frac{{e}_{2}{e}_{3}\sin 2\theta ((s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}})}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right]-4{m}_{\psi }{g}_{\psi }\sin 2\theta \left[\frac{{e}_{2}\sin \theta (s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{3}\cos \theta (s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \times \left[\frac{\arctan h({x}_{i})}{{x}_{i}(s-2{m}_{{h}_{i}}^{2})}[{g}_{p}^{2}({m}_{{h}_{1}}^{2}+{m}_{{h}_{2}}^{2}-3s)-{g}_{\psi }^{2}({m}_{{h}_{1}}^{2}+{m}_{{h}_{2}}^{2}-8{m}_{\psi }^{2}+s)]+{g}_{p}^{2}-{g}_{\psi }^{2}\right]\\ & -{\sin }^{2}2\theta \left[\frac{1}{2{m}_{{h}_{i}}^{4}+2{m}_{\psi }^{2}(s-4{m}_{{h}_{i}}^{2})}\left(-{g}_{p}^{4}[2{m}_{{h}_{i}}^{4}+{m}_{{h}_{i}}^{2}({m}_{{h}_{j}}^{2}-8{m}_{\psi }^{2})+2{m}_{\psi }^{2}s]\right.\right.\\ & \left.+2{g}_{p}^{2}{g}_{\psi }^{2}(2{m}_{{h}_{i}}^{4}+{m}_{{h}_{i}}^{2}({m}_{{h}_{j}}^{2}-10{m}_{\psi }^{2})-2{m}_{\psi }^{2}({m}_{{h}_{j}}^{2}-3s))\right)-{g}_{\psi }^{4}[({m}_{{h}_{i}}^{2}-4{m}_{\psi }^{2})(2{m}_{{h}_{i}}^{2}+{m}_{{h}_{j}}^{2}-4{m}_{\psi }^{2})-2{m}_{\psi }^{2}s]\\ & +\frac{\arctan h({x}_{i})}{{x}_{i}{(s-2{m}_{{h}_{i}})}^{2}}\left({g}_{p}^{4}[3{m}_{{h}_{i}}^{2}+{m}_{{h}_{i}}^{2}(4{m}_{{h}_{j}}^{2}-3s)-{m}_{{h}_{j}}^{4}-{m}_{{h}_{i}}^{2}s+{s}^{2}]\right.\\ & +2{g}_{p}^{2}{g}_{\psi }^{2}[-3{m}_{{h}_{i}}^{2}+{m}_{{h}_{i}}^{2}(3(4{m}_{\psi }^{2}+s)-4{m}_{{h}_{j}}^{4})+{m}_{{h}_{j}}^{4}+{m}_{{h}_{j}}^{2}(s-4{m}_{\psi }^{2})-{s}^{2}]\\ & \left.\left.\left.+{g}_{\psi }^{4}[3{m}_{{h}_{i}}^{2}+{m}_{{h}_{i}}^{2}(-3(8{m}_{\psi }^{2}+s)+4{m}_{{h}_{j}}^{2})-{m}_{{h}_{j}}^{4}-{m}_{{h}_{j}}^{2}(s-8{m}_{\psi }^{2})+{s}^{2}+16{m}_{\psi }^{2}(s-2{m}_{\psi }^{2})]\right)\right]\right\},\end{array}\end{eqnarray} \tag{ A7 }$$

where the auxiliary parameters and coupling constants have the following expressions

$$\begin{eqnarray}&&{x}_{i}=\frac{{(s-4{m}_{\psi }^{2})}^{\frac{1}{2}}{(s-4{m}_{{h}_{i}}^{2})}^{\frac{1}{2}}}{(s-2{m}_{{h}_{i}}^{2})},\quad \quad {\rm{with}}\quad \quad i,j=1,2\quad {\rm{and}}\quad i\ne j,\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}&&{e}_{1}=-{\sin }^{2}\theta ({\alpha }_{s}\sin \theta +6{v}_{H}{\lambda }_{2}\cos \theta )-3{\cos }^{2}\theta ({\lambda }_{1}\sin \theta +2{v}_{H}{\lambda }_{H}\cos \theta ),\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{e}_{2}=\sin 2\theta (-\frac{1}{2}{\alpha }_{s}\sin \theta +3{\lambda }_{H}{v}_{H}\cos \theta )-{\lambda }_{1}\cos \theta (1-3{\sin }^{2}\theta )+2{v}_{H}{\lambda }_{2}\sin \theta (1-3{\cos }^{2}\theta ),\end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray}&&{e}_{3}=\sin 2\theta (-\frac{1}{2}{\alpha }_{s}\cos \theta -3{\lambda }_{H}{v}_{H}\sin \theta )-{\lambda }_{1}\sin \theta (1-3{\cos }^{2}\theta )-2{v}_{H}{\lambda }_{2}\cos \theta (1-3{\sin }^{2}\theta ),\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{e}_{4}=-{\cos }^{2}\theta ({\alpha }_{s}\cos \theta -6{v}_{H}{\lambda }_{2}\sin \theta )-3{\sin }^{2}\theta ({\lambda }_{1}\cos \theta -2{v}_{H}{\lambda }_{H}\sin \theta ).\end{eqnarray} \tag{ A12 }$$

The scalar component will account for the DM phenomenology by the following annihilation cross sections:

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{\bar{f}f\to \chi \chi }=\frac{1}{16\pi }{\left(\frac{{m}_{f}}{{v}_{H}}\right)}^{2}{\left(1-\frac{4{m}_{f}^{2}}{s}\right)}^{3/2}{N}_{c}R,\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{ZZ\to \chi \chi }=\frac{1}{36\pi s}{\left(\frac{{m}_{Z}^{2}}{{v}_{H}}\right)}^{2}\left(2+\frac{{(s-2{m}_{Z}^{2})}^{2}}{4{m}_{Z}^{4}}\right){\left(1-\frac{4{m}_{Z}^{2}}{s}\right)}^{1/2}R,\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&{(\sigma {v}_{{\rm{rel}}})}_{{W}^{+}{W}^{-}\to \chi \chi }=\frac{1}{36\pi s}{\left(\frac{{m}_{W}^{2}}{{v}_{H}}\right)}^{2}\left(2+\frac{{(s-2{m}_{W}^{2})}^{2}}{4{m}_{W}^{4}}\right){\left(1-\frac{4{m}_{W}^{2}}{s}\right)}^{1/2}R\end{eqnarray} \tag{ A15 }$$

where

$$\begin{eqnarray}&&R=\frac{{e}_{5}^{2}{\cos }^{2}\theta }{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{6}^{2}{\sin }^{2}\theta }{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}-\frac{{e}_{5}{e}_{6}\sin 2\theta [(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}]}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}.\end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{1}{h}_{1}\to \chi \chi }= & \frac{1}{16\pi s}{\left(1-\frac{4{m}_{{h}_{1}}^{2}}{s}\right)}^{1/2}\left\{\left[\frac{{e}_{1}^{2}{e}_{5}^{2}}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{2}^{2}{e}_{6}^{2}}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right.\right.\\ & \left.+\frac{2{e}_{1}{e}_{2}{e}_{5}{e}_{6}[(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}]}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right]\\ & -\left(\frac{8{e}_{5}^{2}}{s-2{m}_{{h}_{1}}^{2}}F({y}_{1})+2{e}_{7}\right)\left[\frac{{e}_{1}{e}_{5}(s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{2}{e}_{6}(s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \left.+\frac{8{e}_{5}^{2}}{s-2{m}_{{h}_{1}}^{2}}F({y}_{1})\left[\frac{{e}_{5}^{2}}{s-2{m}_{{h}_{1}}^{2}}\left[\frac{1}{F({y}_{1})(1-{y}_{1}^{2})}+1\right]+{e}_{7}\right]+{e}_{7}^{2}\right\},\end{array}\end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{2}{h}_{2}\to \chi \chi }= & \frac{1}{16\pi s}{\left(1-\frac{4{m}_{{h}_{2}}^{2}}{s}\right)}^{1/2}\left\{\left[\frac{{e}_{3}^{2}{e}_{5}^{2}}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{4}^{2}{e}_{6}^{2}}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right.\right.\\ & \left.+\frac{2{e}_{3}{e}_{4}{e}_{5}{e}_{6}[(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}]}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right]\\ & -\left(\frac{8{e}_{6}^{2}}{s-2{m}_{{h}_{2}}^{2}}F({y}_{2})+2{e}_{8}\right)\left[\frac{{e}_{3}{e}_{5}(s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{4}{e}_{6}(s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \left.+\frac{8{e}_{6}^{2}}{s-2{m}_{{h}_{2}}^{2}}F({y}_{2})\left[\frac{{e}_{6}^{2}}{s-2{m}_{{h}_{2}}^{2}}(\frac{1}{F({y}_{2})(1-{y}_{2}^{2})}+1)+{e}_{8}\right]+{e}_{8}^{2}\right\},\end{array}\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}\begin{array}{ll}{(\sigma {v}_{{\rm{rel}}})}_{{h}_{1}{h}_{2}\to \chi \chi }= & \frac{1}{16\pi s}{\left(1-\frac{4{m}_{{h}_{i}}^{2}}{s}\right)}^{1/2}\left\{\left[\frac{{e}_{2}^{2}{e}_{5}^{2}}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{3}^{2}{e}_{6}^{2}}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right.\right.\\ & \left.+\frac{2{e}_{2}{e}_{3}{e}_{5}{e}_{6}[(s-{m}_{{h}_{1}}^{2})(s-{m}_{{h}_{2}}^{2})+{m}_{{h}_{1}}{m}_{{h}_{2}}{\varGamma }_{{h}_{1}}{\varGamma }_{{h}_{2}}]}{[{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}][{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}]}\right]\\ & -\left(\frac{8{e}_{5}{e}_{6}}{s-2{m}_{{h}_{i}}^{2}}F({y}_{i})+2{e}_{9}\right)\left[\frac{{e}_{2}{e}_{5}(s-{m}_{{h}_{1}}^{2})}{{(s-{m}_{{h}_{1}}^{2})}^{2}+{m}_{{h}_{1}}^{2}{\varGamma }_{{h}_{1}}^{2}}+\frac{{e}_{3}{e}_{6}(s-{m}_{{h}_{2}}^{2})}{{(s-{m}_{{h}_{2}}^{2})}^{2}+{m}_{{h}_{2}}^{2}{\varGamma }_{{h}_{2}}^{2}}\right]\\ & \left.+\frac{8{e}_{5}{e}_{6}}{s-2{m}_{{h}_{i}}^{2}}F({y}_{i})\left[\frac{{e}_{5}{e}_{6}}{s-2{m}_{{h}_{i}}^{2}}(\frac{1}{F({y}_{i})(1-{y}_{i}^{2})}+1)+{e}_{9}\right]+{e}_{9}^{2}\right\},\end{array}\end{eqnarray} \tag{ A19 }$$

where we have employed parameter y_i and function F(y_i) as

$$\begin{eqnarray}&&\,{y}_{i}=\frac{{(s-4{m}_{\chi }^{2})}^{\frac{1}{2}}{(s-4{m}_{{h}_{i}}^{2})}^{\frac{1}{2}}}{(s-2{m}_{{h}_{i}}^{2})}\\ &&F({y}_{i})=\frac{1}{{y}_{i}}\arctan h({y}_{i})\quad \quad \quad {\rm{with}}\quad i=1,2,\end{eqnarray} \tag{ A20 }$$

and also for coupling constants, we have the following parameters:

$$\begin{eqnarray}&&{e}_{5}=-[{v}_{H}{\lambda }_{\chi H}\cos \theta +{\lambda }_{3}\sin \theta ]\times 2!,\end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray}&&{e}_{6}=[{v}_{H}{\lambda }_{\chi H}\sin \theta -{\lambda }_{3}\cos \theta ]\times 2!,\end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray}&&{e}_{7}=-[\frac{1}{2}{\lambda }_{\chi H}{\cos }^{2}\theta +{\lambda }_{4}{\sin }^{2}\theta ]\times 4!,\end{eqnarray} \tag{ A23 }$$

$$\begin{eqnarray}&&{e}_{8}=-[\frac{1}{2}{\lambda }_{\chi H}{\sin }^{2}\theta +{\lambda }_{4}{\cos }^{2}\theta ]\times 4!,\end{eqnarray} \tag{ A24 }$$

$$\begin{eqnarray}&&{e}_{9}=\left[\left(\frac{1}{2}{\lambda }_{\chi H}-{\lambda }_{4}\right)\sin 2\theta \right]\times 2!.\end{eqnarray} \tag{ A25 }$$

Finally, the decay rates of scalars h_i (with i = 1,2) into fermionic and scalar DM particles are given as:

$$\begin{eqnarray}&&\varGamma ({h}_{i}\to \bar{\psi }\psi )=\frac{{S}_{i}^{2}\theta }{8\pi }{(1-4{m}_{\psi }^{2}/{m}_{{h}_{i}}^{2})}^{\frac{3}{2}}\left[{g}_{\psi }^{2}+\frac{{m}_{{h}_{i}}^{2}{g}_{p}^{2}}{{m}_{{h}_{i}}^{2}-4{m}_{\psi }^{2}}\right],\end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray}&&\varGamma ({h}_{1}\to \chi \chi )=\frac{{e}_{5}^{2}}{32\pi {m}_{{h}_{1}}}{(1-4{m}_{\chi }^{2}/{m}_{{h}_{1}}^{2})}^{\frac{1}{2}},\end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray}&&\varGamma ({h}_{2}\to \chi \chi )=\frac{{e}_{6}^{2}}{32\pi {m}_{{h}_{2}}}{(1-4{m}_{\chi }^{2}/{m}_{{h}_{2}}^{2})}^{\frac{1}{2}},\end{eqnarray} \tag{ A28 }$$

where we have defined S₁θ = sinθ and S₂θ = cosθ.