Non-Abelian Vector Boson as FIMP Dark Matter

TL;DR

This paper demonstrates a freeze-in mechanism for non-abelian vector boson dark matter within an extended Standard Model, exploring its production, stability, and observational constraints, and highlighting its unique collider signatures.

Contribution

It introduces a novel freeze-in dark matter scenario with a non-abelian vector boson in an extended SM, including scalar and neutrino mass generation sectors.

Findings

Scalar decay channels significantly contribute to DM production.

The model's parameter space aligns with CMB, BBN, and AMS-02 constraints.

Distinct collider signatures differentiate freeze-in from freeze-out DM scenarios.

Abstract

In this analysis we demonstrate the freeze-in realization of a non-abelian vector boson dark matter (DM). We choose to elaborate an existing $SU(2)_N$ extension ($N$ stands for neutral) of the Standard Model (SM) with an additional $U(1)=S^{'}$ global symmetry, which stabilizes the vector boson ($X,\bar{X}$) as DM through unbroken $S=T_{3N}+S^{'}$ and as lightest odd $S$ particle. The analysis reveals that the contribution to the freeze-in production of DM from the decay of a heavier scalar bidoublet $\zeta_1^{0,\pm} \to \zeta_2^{0,\pm}X$ is important even after the freeze-out of $\zeta_1^{0,\pm}$ in equilibrium with thermal bath. Moreover, the neutral component of $SU(2)_N$ scalar triplet ($\Delta$), responsible for neutrino mass generation in this framework, turns out to serve as additional DMs in the model and offers a multipartite freeze-in DM set up to explore. The allowed parameter space is obtained after estimating constraints from CMB, BBN and AMS-02 bound. This exercise nicely complements the freeze-out realization of ($X,\bar{X}$) as weakly interacting massive particle (WIMP) and distinguishes it through stable charge track signature at collider compared to leptonic signal excess as found in WIMP scenario.

Tables

Table 1. Table 1 : Relevant particle content of the model and their charges under S ​ U ​ ( 3 ) c ⊗ S ​ U ​ ( 2 ) L ⊗ U ​ ( 1 ) Y ⊗ S ​ U ​ ( 2 ) N ⊗ S ′ tensor-product tensor-product tensor-product tensor-product 𝑆 𝑈 subscript 3 𝑐 𝑆 𝑈 subscript 2 𝐿 𝑈 subscript 1 𝑌 𝑆 𝑈 subscript 2 𝑁 superscript 𝑆 ′ SU(3)_{c}\otimes SU(2)_{L}\otimes U(1)_{Y}\otimes SU(2)_{N}\otimes S^{{}^{\prime}} . S ​ U ​ ( 2 ) L 𝑆 𝑈 subscript 2 𝐿 SU(2)_{L} doublets are indicated by vertical parenthesis with T 3 ​ L = ± 1 / 2 subscript 𝑇 3 𝐿 plus-or-minus 1 2 T_{3L}=\pm 1/2 for up and down components respectively. S ​ U ​ ( 2 ) N 𝑆 𝑈 subscript 2 𝑁 SU(2)_{N} doublet is depicted by entries in horizontal parenthesis with T 3 ​ N = ± 1 / 2 subscript 𝑇 3 𝑁 plus-or-minus 1 2 T_{3N}=\pm 1/2 to left and right fields respectively.

Particles	$SU{\left(3\right)}_c$	$SU{\left(2\right)}_L$	$U{\left(1\right)}_Y$	$SU{\left(2\right)}_N$	$S^{{}^{\prime }}$
$X_{1,2,3}$	1	1	0	3	0
$n={(n_1,n_2)}_{L,R}$	1	1	0	2	1/2
$\chi \equiv ({\chi }_1,{\chi }_2)$	1	1	0	2	1/2
	1	2	-1/2	2	-1/2
	1	1	0	3	-1
$\Phi \equiv \left(\begin{array}{c}{\varphi }^{+}\\ {\varphi }^0\end{array}\right)$	1	2	1/2	1	0

Table 2. Table 2 : S 𝑆 S charge assignment for the new particles added in the model as introduced in Tab.  1 .

Particles	$S^{{}^{\prime }}$	$S=S^{{}^{\prime }}+T_{3N}$
$X\left(\overline{X}\right)$	0	+1(-1)
$X_3$		0
$n_{1L,R}$	1/2	+1
$n_{2L,R}$		0
${\chi }_1$	1/2	+1
${\chi }_2$		0
${\zeta }_1$	-1/2	-1
${\zeta }_2$		0
${\Delta }_1$		-2
${\Delta }_2$	-1	-1
${\Delta }_3$		0

Table 3. Table 3 : Spectrum of physical particles appearing at different scales assuming g N ∼ 10 − 13 similar-to subscript 𝑔 𝑁 superscript 10 13 g_{N}\sim 10^{-13} , compatible with freeze-in requirement of X ​ ( X ¯ ) 𝑋 ¯ 𝑋 X(\bar{X}) . In the dark matter analysis we will further assume f 1 ≈ f 2 subscript 𝑓 1 subscript 𝑓 2 f_{1}\approx f_{2} , such that both ζ 1 − , 0 superscript subscript 𝜁 1 0 \zeta_{1}^{-,0} are brought down to ∼ similar-to \sim TeV scale.

Scale	GeV	GeV, $f_1=f_2$
$𝒪\left({10}^{13}\mathrm{GeV}\right)$	${\Delta }_{1,2,3},{\zeta }_1^{-},{\xi }_1^0,{\chi }_2^R$	${\Delta }_{1,2,3},{\chi }_2^R$
$𝒪\left(1\mathrm{TeV}\right)$	$H^{+},{\eta }^0,{\xi }_2^0,n_1,n_2$	$H^{+},{\eta }^0,{\xi }_1^0,{\xi }_2^0,{\zeta }_1^{-},n_1,n_2$
$𝒪\left(100\mathrm{GeV}\right)$	$W^{\pm },Z,H,t$
$𝒪\left(1\mathrm{GeV}\right)$	$X,Z^{\prime },b,c,\tau$
	$e,\mu ,{\nu }_e,{\nu }_{\mu },s,u,d$

Table 4. Table 4 : Different production scenarios for Δ Δ \Delta ’s, and their dependence on f 7 , 8 subscript 𝑓 7 8 f_{7,8} couplings. Throughout this study the freeze-in of X 𝑋 X is assumed.

Degenerate $\Delta$ ($f_7=0$) Non-degenerate $\Delta$ ($f_7\ne 0$) Production mechanism of $\Delta$ Freeze-Out Freeze-In Freeze-Out Freeze-In Scalar quartic coupling $f_8$ $𝒪\left(1\right)$ $𝒪\left({10}^{-12}\right)$ $𝒪\left(1\right)$ $𝒪\left({10}^{-12}\right)$

Table 5. Table 5 : Decay branching ratio of ζ 1 ± → ζ 1 0 , ℓ ± → superscript subscript 𝜁 1 plus-or-minus superscript subscript 𝜁 1 0 superscript ℓ plus-or-minus \zeta_{1}^{\pm}\to\zeta_{1}^{0},\ell^{\pm} via W 𝑊 W and via n 1 ​ R subscript 𝑛 1 𝑅 n_{1R} . Here we have chosen m ζ 1 = 1.2 ​ TeV subscript 𝑚 subscript 𝜁 1 1.2 TeV m_{\zeta_{1}}=1.2~{}\rm TeV and δ ​ m = 100 ​ MeV 𝛿 𝑚 100 MeV \delta m=100~{}\rm MeV with RHNs of mass M = 10 5 ​ GeV 𝑀 superscript 10 5 GeV M=10^{5}~{}\rm GeV .

$\delta m$ (GeV)	$Br\left({\zeta }_1^{\pm }\to {\zeta }_1^0,{\ell }^{\pm },{\nu }_L\right)$ (via $W$)	$Br\left({\zeta }_1^{\pm }\to {\zeta }_1^0,{\ell }^{\pm },{\nu }_L\right)$ (via $n_{1R}$)
0.1	0.99	0.0007
0.5	0.99	0.0007

Table 6. Table 6 : Decay lengths for the three-body decay of ζ 1 ± superscript subscript 𝜁 1 plus-or-minus \zeta_{1}^{\pm} assuming maximum the radiative mass splitting δ ​ m 𝛿 𝑚 \delta m of 500 MeV for m ζ 1 ± = 1.2 ​ TeV subscript 𝑚 superscript subscript 𝜁 1 plus-or-minus 1.2 TeV m_{\zeta_{1}^{\pm}}=1.2~{}\rm TeV , RHN mass M = 10 5 ​ GeV 𝑀 superscript 10 5 GeV M=10^{5}~{}\rm GeV and f ζ ∼ 𝒪 ​ ( 1 ) similar-to subscript 𝑓 𝜁 𝒪 1 f_{\zeta}\sim\mathcal{O}(1) .

$\delta m$ (GeV)	$\Gamma$ (GeV)	$L=c\tau$ (m)
0.1	$3.12\times {10}^{-18}$	63.05
0.3	$7.58\times {10}^{-16}$	0.26
0.5	$9.76\times {10}^{-15}$	0.20

Equations

$$\begin{array}{r}{m}_W^2=\frac{1}{2}{g}_2^2\left(v_1^2+v_2^2\right),m_X^2=\frac{1}{2}{g}_N^2\left(u_2^2+v_2^2+2u_3^2\right),m_{Z^{\prime }}^2\simeq \frac{1}{2}{g}_N^2\left(u_2^2+v_2^2+4u_3^2\right),\end{array}$$

$$m_{Z,Z^{\prime }}^2=\frac{1}{2}\left(\begin{array}{cc}\left(g_1^2+g_2^2\right)\left(v_1^2+v_2^2\right)& -g_N\sqrt{g_1^2+g_2^2}{v}_2^2\\ -g_N\sqrt{g_1^2+g_2^2}{v}_2^2& g_N^2\left(u_2^2+v_2^2+4u_3^2\right)\end{array}\right).$$

$$\begin{array}{rl}V& ={\mu }_{\zeta }^{2}Tr({\zeta }^{†}\zeta )+{\mu }_{\Phi }^2{\Phi }^{†}\Phi +{\mu }_{\chi }^2{\chi }^{†}\chi +{\mu }_{\Delta }^{2}Tr({\Delta }^{†}\Delta )+({\mu }_1{\tilde{\Phi }}^{†}\zeta \chi +{\mu }_2{\tilde{\chi }}^{†}\Delta \chi +H.c.)\\ & +\frac{1}{2}{\lambda }_1[Tr({\zeta }^{†}\zeta ){]}^2+\frac{1}{2}{\lambda }_2({\Phi }^{†}\Phi {)}^2+\frac{1}{2}{\lambda }_{3}Tr({\zeta }^{†}\zeta {\zeta }^{†}\zeta )+\frac{1}{2}{\lambda }_4({\chi }^{†}\chi {)}^2+\frac{1}{2}{\lambda }_5[Tr({\Delta }^{†}\Delta ){]}^2\\ & +\frac{1}{4}{\lambda }_{6}Tr({\Delta }^{†}\Delta -\Delta {\Delta }^{†}{)}^2+f_1{\chi }^{†}{\tilde{\zeta }}^{†}\tilde{\zeta }\chi +f_2{\chi }^{†}{\zeta }^{†}\zeta \chi +f_3{\Phi }^{†}\zeta {\zeta }^{†}\Phi +f_4{\Phi }^{†}\tilde{\zeta }{\tilde{\zeta }}^{†}\Phi \\ & +f_5({\Phi }^{†}\Phi )({\chi }^{†}\chi )+f_6({\chi }^{†}\chi )Tr({\Delta }^{†}\Delta )+f_7{\chi }^{†}(\Delta {\Delta }^{†}-{\Delta }^{†}\Delta )\chi +f_8({\Phi }^{†}\Phi )Tr({\Delta }^{†}\Delta )\\ & +f_{9}Tr({\zeta }^{†}\zeta )Tr({\Delta }^{†}\Delta )+f_{10}Tr[\zeta ({\Delta }^{†}\Delta -\Delta {\Delta }^{†}){\zeta }^{†}],\end{array}$$

$${\tilde{\Phi }}^{†}=({\varphi }^0,-{\varphi }^{+}),{\tilde{\chi }}^{†}=({\chi }_2,-{\chi }_1),\tilde{\zeta }=\left(\begin{array}{cc}{\zeta }_2^{+}& -{\zeta }_1^{+}\\ -{\bar{\zeta }}_2^0& {\bar{\zeta }}_1^0\end{array}\right).$$

$$u_2^2$$

$$v_1^2$$

$$v_2$$

$$u_3$$

$$m_{{\Delta }_1}^2$$

$$m_{{\Delta }_2}^2$$

$$m_{{\Delta }_3}^2$$

$$h$$

$${\xi }_2^0$$

$${\xi }_1^0$$

$${\eta }^0$$

$$m_h^2$$

$$m_{{\xi }_2^0}^2$$

$$m_{{\eta }^0}^2$$

$$m_{{\xi }_1^0}^2$$

$$m_{{\chi }_2^R}^2$$

$$H^{+}=\frac{1}{\sqrt{v_1^2+v_2^2}}\left(v_2{\varphi }^{+}+v_1{\zeta }_2^{+}\right).$$

$$m_{{\zeta }_1^{-}}^2$$

$$m_{H^{+}}^2$$

$$\cos (\beta -\alpha )\approx (\sqrt{2}-1)\frac{v_2}{v_1}+\mathcal{O}(\frac{{\lambda }_2{v}_2^2}{{\mu }_1{u}_2})\le 1.7\times 10^{-3}$$

$$f_{\zeta }$$

$$f_{\Delta }$$

$$f_{\zeta }{v}_2{\overline{\nu }}_L{n}_{2R}-f_{\Delta }^L{u}_3{n}_{2L}{n}_{2L}-f_{\Delta }^R{u}_3{n}_{2R}{n}_{2R}+\text{h.c.}$$

$$M_{\nu }=\left(\begin{array}{ccc}0& m_D& 0\\ m_D& m_2^{\prime }& M\\ 0& M& m_2\end{array}\right),$$

$$m_{\nu }\simeq \frac{m_D^2{m}_2}{M^2}=f_{\zeta }^2{f}_{\Delta }{\left(\frac{v_2}{M}\right)}^2{u}_3.$$

$$\begin{array}{rl}{\dot{n}}_{{\Delta }_i}+3H{n}_{{\Delta }_i}=& \int d{\Pi }_{\Phi }d{\Pi }_{{\Phi }^{\ast }}d{\Pi }_{{\Delta }_i}d{\Pi }_{{\Delta }_i^{\ast }}(2\pi {)}^4{\delta }^4(p_{\Phi }+p_{{\Phi }^{\ast }}-p_{{\Delta }_i}-p_{{\Delta }_i^{\ast }})\\ \times & [|\mathcal{M}{|}_{\Phi {\Phi }^{\ast }\to {\Delta }_i{\Delta }_i^{\ast }}^2{f}_{\Phi }{f}_{{\Phi }^{\ast }}(1+f_{{\Delta }_i})(1+f_{{\Delta }_i^{\ast }})\\ & -|\mathcal{M}{|}_{{\Delta }_i{\Delta }_i^{\ast }\to \Phi {\Phi }^{\ast }}^2{f}_{{\Delta }_i}{f}_{{\Delta }_i^{\ast }}(1+f_{\Phi })(1+f_{{\Phi }^{\ast }})],\end{array}$$

$$n_i=\frac{g_i}{(2\pi {)}^3}\int d^{3}p{f}_i(p),$$

$$\begin{array}{rl}{\dot{n}}_{{\Delta }_i}+3H{n}_{{\Delta }_i}=& \frac{T}{512{\pi }^6}{\int }_{4m_{\Delta }^2}^{\infty }dsd\Omega |\mathcal{M}{|}_{\Phi {\Phi }^{\ast }\to {\Delta }_i{\Delta }_i^{\ast }}^2{P}_{\Phi {\Phi }^{\ast }}{P}_{{\Delta }_i{\Delta }_i^{\ast }}{K}_1(\sqrt{s}/T)/\sqrt{s},\end{array}$$

$$P_{ij}\equiv \frac{{\left[s-(m_i+m_j{)}^2\right]}^{1/2}{\left[s-(m_i-m_j{)}^2\right]}^{1/2}}{2\sqrt{s}}=\sqrt{\frac{s}{4}-m_i^2}.$$

$$\begin{array}{rl}\mathcal{S}\cdot {\dot{Y}}_{\Delta }=& \frac{T\cdot f_8^2}{512{\pi }^5}{\int }_{4m_{\Delta }^2}^{\infty }ds\sqrt{s-4m_{\Delta }^2}\sqrt{s-4m_{\Phi }^2}{K}_1(\sqrt{s}/T)/\sqrt{s},\end{array}$$

$$\begin{array}{rl}\frac{d{Y}_{\Delta }}{dx}=& \frac{45M_{\text{pl}}{f}_8^2{x}^3}{1024g_{\star }^S\sqrt{g_{\star }^{\rho }}1.66{\pi }^7{m}_{\Delta }^4}{\int }_{4m_{\Delta }^2}^{\infty }ds\sqrt{s-4m_{\Delta }^2}{K}_1(x\cdot \sqrt{s}/m_{\Delta }),\end{array}$$

$$\begin{array}{rl}{Y}_{{\Delta }_i}=& \frac{0.10588M_{\text{pl}}{f}_8^2}{{\pi }^7{m}_{{\Delta }_i}}{\int }_0^{\infty }dx\frac{x^2}{g_{\star }^S\sqrt{g_{\star }^{\rho }}}\cdot K_1(x{)}^2\approx 3.244\times 10^{-7}\left(\frac{M_{\text{pl}}}{m_{{\Delta }_i}}\right)f_8^2,\end{array}$$

$${\Omega }_{\Delta }\cdot h^2=3\times \frac{2970m_{\Delta }{Y}_{{\Delta }_i}^{\infty }{\text{cm}}^{-3}}{1.88\times 10^{-29}\text{g}\cdot {\text{cm}}^{-3}}\approx 10^{21}{f}_8^2$$

Full text

1 Introduction

A particle dark matter (DM) is highly motivated from astrophysical observations. However, laboratory experiments like direct search and collider searches have not detected any signals yet. It is therefore an important exercise to look for possibilities where DM interaction with the Standard Model (SM) is suppressed, yet produces correct relic density as has been dictated by anisotropies in Cosmic Microwave Background Radiation (CMBR). Among several possibilities [1, 2], freeze-in mechanism serves as an interesting alternative [3, 4, 5, 6, 7, 8]. In such a case, the DM is assumed not to be in equilibrium with thermal bath in the early universe owing to its tiny coupling to the visible sector. It is then produced then non-thermally via decay or annihilation of particles in thermal bath and freezes in once the temperature drops below DM mass to yield correct relic density ($\Omega h^{2}\sim 0.1198$ [9]). Several studies have been done in this direction to show that freeze-in can give rise to DM mass $\sim$ TeV scale, but the coupling with visible sector requires to be extremely tiny \mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10^{-10}. Therefore, such DM models are classified as feebly interacting massive particle (FIMP) (for a review see [10]), as opposed to the thermal freeze-out of weakly interacting massive particle (WIMP).

Our aim here is to demonstrate the freeze-in of a non-abelian vector boson DM [11] (see [12] for an abelian example). We choose a well motivated $SU(2)_{N}$ extension ($N$ stands for neutral111This means that the $SU(2)_{N}$ vector bosons are electromagnetic charge neutral.) of the Standard Model (SM) with an additional $U(1)=S^{{}^{\prime}}$ global symmetry, which stabilizes the lightest of the vector boson ($X,\bar{X}$) as DM through unbroken $S=T_{3N}+S^{{}^{\prime}}$ [13]. Spontaneous symmetry breaking of $SU(2)_{N}$ results in massive gauge bosons. None of the SM particles are charged under $SU(2)_{N}$ and therefore $X,\bar{X}$ do not have a direct coupling to the visible sector besides Higgs portal which is required to be tiny to avoid conflict with Higgs data. This naturally leads to the possibility of FIMP nature of $X,\bar{X}$ as DM. The model, in addition, possesses several features, for example, addressing neutrino mass generation via inverse seesaw of type-III irrespective of whether the DM is undergoing freeze-in or freeze-out and a possible unification to $SU(7)$.

The DM $(X,\bar{X})$ in this model is produced from the decay of a scalar bi-doublet $\zeta_{1}^{0,\pm}$, owing to $SU(2)_{N}$ interaction222If the decay is kinematically forbidden, then annihilation of $\zeta_{1}^{0,\pm}$ to produce $(X,\bar{X})$ becomes important.. Naturally, the correct relic density of $X,\bar{X}$ via freeze-in indicates that the $SU(2)_{N}$ gauge coupling needs to be $\sim\mathcal{O}(10^{-10})$. One of the salient features of this study is to show that the contribution to DM production from the decay of $\zeta_{1}^{0,\pm}$ remains significant even after the freeze-out of $\zeta_{1}^{0,\pm}$. We point out that such a feature is inevitable whenever the decay is slow enough, although not much elaborated in FIMP literature. Moreover, our formulae for computing the yield of the DM from the decay of a decoupled species is actually model-independent, and can be applied to any scenario where the freeze-in production of the DM is taking place via the decay of a heavy particle, which shared a common temperature with the thermal bath in the early universe. The subsequent change in the allowed parameter space due to the ‘late decay’ turns out to be quite noteworthy. The $SU(2)_{N}$ scalar triplet ($\Delta$) required for neutrino mass generation also provides with additional DM components in this model. The neutral components of $\Delta$ turn out to be stable at the scale of the universe life time, thanks to the small $SU(2)_{N}$ gauge coupling in the freeze-in mechanism advocated here. Therefore the model also serves as a multipartite FIMP DM set up, although the DM components do not have sizable interaction with each other. Hence, the freeze-in of each individual component remains unaffected by the presence of others, but provides with a sizable range of allowed parameter space to span the whole under abundant region in the absence of direct search constraint for these FIMP like DMs.

However, stringent bounds on the lifetime of semi-stable charged and neutral particles ($\zeta_{1,2}^{0,\pm}$) arises from Big Bang Nucleosynthesis (BBN), which essentially rules out hadronically decaying particle with lifetime $\tau>100~{}\rm sec$ [14, 15]. CMB, on the other hand, puts a lower bound on the lifetime of DM decaying to SM particles, which can potentially alter the ionization history (and hence the power spectrum) of CMB [16]. Experiments like AMS-02 [17] also puts lower bound on decaying DM from non-observation of anti proton excess. In our case, as we shall elaborate, bounds from CMB and AMS-02 are rather lose but BBN plays a crucial role by eliminating a large portion of the parameter space allowed from relic density condition.

This same model from WIMP perspective, has already been thoroughly explored in [13, 18]. The present exercise therefore provides an opportunity to compare two different realizations (freeze out versus freeze in) of the same model. In the following, we see that this provides not only a distinction in terms of $SU(2)_{N}$ gauge coupling, but also in terms of DM mass. For example, in case of freeze-in, we are bound to stick to low DM mass: m_{X}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}50~{}\rm GeV (depending on the contribution of $X$ to total relic abundance), while the WIMP scenario is valid for DM mass even upto $\sim$ TeV. Finally, the distinction between the WIMP and FIMP realization can also arise in collider signature of the model. For WIMP case, it was shown that hadronically quiet single and two lepton channels could verify the existence of such a model framework. On the contrary, in FIMP realization, the signature can arise through stable charge tracks or displaced vertices of $\zeta_{1,2}^{\pm}$ as demonstrated here.

The paper is organised as follows: in Sec. 2 we have provided the details of the model including the symmetry breaking and spectrum of the physical particles that are important for the present phenomenology. Then in Sec. 3 we have shown how light neutrino mass can be generated via inverse seesaw mechanism in this framework. Sec. 4 contains the main DM analysis under which in subsection  4.1 and 4.1.1 we have discussed in detail the yield for $X$ and $\Delta$ by solving the Boltzmann equation (BEQ) and in subsection 4.1.3 we have elaborated the impact of BBN and CMB bounds on the parameter space of this model. In Sec. 5 we have shown the possible signatures that this model may yield at the colliders. Finally in Sec. 6 we have summarized our findings.

2 The Model

We have considered a $SU(2)_{N}$ extension of the SM ($N$ stands for neutral), where the lightest of the gauge bosons acts as a DM candidate. The particle content is chosen in such a way so that the spontaneous symmetry breaking (SSB) of $SU(2)_{N}$ to yields massive gauge bosons and at the same time it is also possible to generate correct light neutrino mass successfully as proposed in [13]. All the SM fermions are singlet under the new $SU(2)_{N}$. The stability of DM is ensured by an imposed global $U(1)$ symmetry ($S^{{}^{\prime}}$), such that $S=S^{{}^{\prime}}+T_{3N}$ remains unbroken.

The new particles introduced in the model and their charges under $SU(3)_{C}\,\otimes$ $SU(2)_{L}\,\otimes$ $U(1)_{Y}\otimes$ $SU(2)_{N}\otimes\,S^{{}^{\prime}}$ are noted in Tab. 1. In the gauge sector, there are three $SU(2)_{N}$ gauge bosons $X_{1,2,3}$, where $X(\overline{X})=\frac{X_{1}\mp iX_{2}}{\sqrt{2}}$ turns out to be mass degenerate and serves as DM component(s) of the model. In the fermion sector, three families of Dirac fermion doublets $n=\left(n_{1},n_{2}\right)_{L,R}$ (under $SU(2)_{N}$) are introduced which mediate the interactions of the dark sector (non-zero $S$ charged particles as noted below) with the SM sector. The scalar sector consists of one $SU(2)_{L}$ scalar doublet $\Phi$ (which contains the 125 GeV Higgs boson), one $SU(2)_{N}$ doublet $\chi$, one scalar bi-doublet $\zeta$ and one $SU(2)_{N}$ scalar triplet $\Delta$. The Dirac fermion doublets, together with the scalar triplet ($\Delta$) are required for generating light neutrino masses, which shall be discussed in the next section. The scalar doublet and bi-doublets participate in spontaneous symmetry breaking (SSB) to generate masses for all the particles involved in the model. The minimization condition along with the physical states that appear after SSB are elaborated below. It is important to note here that $SU(2)_{N}\otimes S^{{}^{\prime}}\to S(=S^{\prime}+T_{3N})$ occurs via the non-zero vacuum expectation value (VEV) of $SU(2)_{N}$ scalar doublet: $\langle\chi_{2}\rangle=u_{2}$. In Tab. 2 we have tabulated the $S$ charge assignments for the new particles added in the model. All the SM particles have zero $S$ charge. Therefore, particles with non-zero $S$ charge will be protected from decaying into the SM. We can assume $X(\bar{X})$ to be the lightest of the particles with non-zero $S$ charge to qualify as DM candidate(s). Furthermore, neutral components of scalar triplet $\Delta_{1,2,3}$ can be stable (if $SU(2)_{N}$ coupling is assumed to be very small for successful freeze-in of $X(\bar{X})$) and be part of a multi-component DM framework.

The scalars which acquire VEV are: $\langle\chi_{2}\rangle=u_{2}$, $\langle\zeta_{2}^{0}\rangle=v_{2}$, $\langle\Delta_{3}\rangle=u_{3}$, and $\langle\phi^{0}\rangle=v_{1}$. Note that the VEV assignment is different here from the $SU(2)_{N}$ extension considered in [19], where $\langle\Delta_{1}^{0}\rangle$ is also non-zero. One should note here that this particular assignment of VEVs for different scalar multiplets of the model is the only possibility to keep the $SM\times S$ symmetry intact, so that VEVs are assigned to only those neutral scalars which have zero $S$-charge; exception of which can lead to breaking of the $S$-charge spoiling the underlying symmetry of the Lagrangian and render the vector boson DM unstable. For example, $SU(2)_{N}$ breaking is only possible through $\langle\chi_{2}\rangle$. This, in turn, ensures that only $\zeta_{2}$ can acquire a non-zero VEV due to the term $(\mu_{1}\tilde{\Phi}^{\dagger}\zeta\chi+H.c.)$ in the scalar potential (see Eq. 2.3). Similarly, it is straightforward to show from the term $\left(\mu_{2}\tilde{\chi}^{\dagger}\Delta\chi+H.c.\right)$, that only $\Delta_{3}$ can acquire a non-zero VEV to keep $S$-symmetry unbroken. One can also check by computing the Hessian333The Hessian ($\mathcal{H}$) is a multi-dimensional matrix in the field space of the scalar potential, consisting of second-order partial derivatives with respect to the fields ($\phi_{i}$): $\mathcal{H}=\frac{\partial^{2}V}{\partial\phi_{i}\partial\phi_{j}}$. of the scalar potential that with our choice of the VEVs, the potential satisfies the necessary condition to be at its local minimum.

With this VEV choices, $X_{1,2}$ bosons have equal masses in this model, and more importantly $S=S^{\prime}+T_{3N}$ global symmetry remains unbroken unlike in [19] as mentioned in the last paragraph. Also note that $U(1)_{EM}$ remains unbroken after the electroweak symmetry breaking as no scalar with EM charge receives a non-zero VEV. This ensures that the photon remains massless. The masses of the other gauge bosons are then given by:

[TABLE]

where $Z-Z^{\prime}$ mixing matrix is given by:

[TABLE]

The choice of the VEVs and couplings required in the present scenario is mainly dictated by the DM phenomenology. Since in this work we are interested to see the freeze-in aspect of the $SU(2)_{N}$ vector bosons (as mentioned in Sec. 1), we require the gauge coupling $g_{N}$ to be extremely small (\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}\mathcal{O}(10^{-10})) in order to keep the lightest gauge boson out of equilibrium. We will carefully evaluate the correct order of the coupling and explore effects of such small coupling in the DM phenomenology. An immediate consequence of this is that g_{N}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}\mathcal{O}(10^{-10}) results in a very large $u_{2}\sim\mathcal{O}(10^{10})~{}\rm GeV$ for $m_{X}\sim\mathcal{O}(\rm TeV)$. Therefore $v_{2}\ll u_{2}$ and $g_{N}\ll g_{1,2}$, which ensures small $Z$-$Z^{{}^{\prime}}$ mixing [20]. This hides $Z^{\prime}$ of this model from being observed at the LHC, and adds to the freedom of choosing $m_{Z^{\prime}}$ as a free parameter. This should again be contrasted to the case in [21], where there is a minimum limit on $M_{X_{1,2,3}}\geqslant 1$ TeV, f or the degenerate vector boson DM case to respect the bound from $Z^{{}^{\prime}}$ search data. Furthermore, owing to $u_{3}\ll u_{2}$, it is clear that the $X$ boson masses are nearly degenerate, i.e. $m_{Z^{\prime}}(m_{X_{3}})\simeq m_{X}$.

With this particle content at our disposal, we can write the most general scalar potential as [13]:

[TABLE]

where

[TABLE]

Since $u_{3},v_{1},v_{2}\ll u_{2}$ as argued above, the minimization conditions yield the following conditions for the VEVs to the scalar potential parameters:

[TABLE]

where from the last line it follows that unless $u_{2}\ll\mu_{\Delta}$, we should have $\mu_{2}\approx u_{3}$.

The SU$(2)_{N}$ triplet scalars $\Delta_{1,2,3}$, remain complex with masses given by444See Appendix A.1 for more details.

[TABLE]

Note that since $\mu_{2}\approx u_{3}$, $\Delta_{3}$ is at the same scale as $u_{2}$, and $f_{7}$ determines the mass differences between $\Delta_{1},~{}\Delta_{2}~{}\rm{and}~{}\Delta_{3}$ as before [18]. The rest of the neutral physical scalars of our model are:

[TABLE]

with masses given by:

[TABLE]

From the two charged scalars, $\zeta_{1}^{-}$ doesn’t mix and is physical, while $H^{+}$ is an admixture of $\phi^{+}$ and $\zeta_{2}^{+}$:

[TABLE]

The masses of the charged scalars in terms of the scalar potential parameters and VEVs are as follows:

[TABLE]

The VEV $u_{2}$ being very large, dictates the scalar masses predominantly. Note that $\chi_{2}$ and $\Delta_{1,2,3}$ are at the same scale as $u_{2}$. For other scalars, dominant contribution arises from the presence of terms proportional to $\mu_{1}u_{2}v_{1}/v_{2}$ (except the SM Higgs). If we require the masses of these new scalars to be at TeV scale, we should have $m\sim 10\sqrt{\mu_{1}u_{2}}\sim 1$ TeV, which requires |\mu_{1}|\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10^{-9} GeV. In that case, $H^{+},\eta^{0},\xi_{2}^{0}$ can all be at $\mathcal{O}(1\text{ TeV})$ scale. Without any further assumptions $\Delta_{1,2,3},~{}\xi_{1}^{0},~{}\zeta_{1}^{-},~{}\chi_{2}^{R}$ are all around the $u_{2}$ scale and very heavy: $\mathcal{O}(10^{13}\text{ GeV})$.

However, as we will demonstrate the freeze-in of $X$ is mainly dictated by $\zeta_{1}^{0,-}$ decays, and a large $m_{\zeta_{1}^{0,-}}$ mass will produce over abundance of the DM. In order to avoid the overabundance of $X$, we need m_{\zeta_{1}^{0,-}}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10\text{ TeV}; which can be achieved by setting $f_{2}=f_{1}$ and cancelling the over powering $u_{2}$ term (see Eq. 2.22). In Tab. 3 we have listed the physical particles in our model (including SM) with their corresponding mass scales that fit the freeze-in requirement for the vector boson DM. It is also important to mention here that the correct DM relic density via freeze in of $X$ allows $m_{X}$ to vary in a large range between few GeV’s to few hundred GeV’s, but the heavier DM masses are constrained by BBN data. Therefore we have listed them in the range of $\mathcal{O}(1~{}\rm GeV)$ in Tab. 3. Similarly, $n_{1,2}$ masses can also be as large as $\sim\mathcal{O}(10^{5})$ GeV depending on the Yukawa coupling as we demonstrate in the next section although they have been classified to lie in $\mathcal{O}(1~{}\rm TeV)$ in the Tab. 3.

We would like to verify the alignment limit of the chosen vevs in the model. At low energies, we can identify SU(2)L scalar doublets $(\Phi,\tilde{\zeta_{2}})$ of our model with $(\Phi_{1},\Phi_{2})$ of the the usual two Higgs doublet model (2HDM) [22, 23, 24]. In our case $\tan\beta=v_{2}/v_{1}\ll 1$, while $\tan\alpha\approx-v_{1}/v_{2}$, where $\beta$ and $\alpha$ diagonalize the mass-squared matrices of the scalars and pseudoscalars respectively. Therefore, we have an approximate alignment as in Type-I case of 2HDM i.e. $\cos(\beta-\alpha)\approx 0$, in the decoupling limit ($v_{2}\ll v_{1}\ll u_{2}$). We would like to emphasize here that this alignment results in recovering a CP-even scalar mass eigenstate with similar gauge, Yukawa and self interactions at tree level as those of the SM Higgs boson555Note that, since no scalar induces a charge breaking VEV, hence the photon remains identically massless independent of the decoupling limit.. The leading contribution to this approximation is given by:

[TABLE]

which satisfies the CMS limits on type-I 2HDM safely lies within the CMS limits on Higgs couplings [25].

3 Neutrino mass

As already elaborated in [13, 18], the generation of light neutrino mass is a novel feature of this model addressed together with DM. In principle, the light neutrino mass generation mechanism is independent of DM phenomenology. However, we will take a quick tour of the neutrino sector here and advocate a subtle phenomenological connection to the freeze-in prospect of the DM, in the light of mass scales introduced in Tab. 3.

The gauge and $S$ invariant Yukawa terms responsible for neutrino mass generation in this model are given by:

[TABLE]

The lepton number is conserved in (3.1) with $n$ carrying $L=1$, and is broken to lepton parity, i.e. $(-1)^{L}$ by the $nn$ terms in (3.2). After SSB we have the following mass terms for the neutrinos:

[TABLE]

where $f_{\zeta}$ and $f_{\Delta}$ are $3\times 3$ matrices and $f_{\Delta}$ is further classified to address left handed ($f_{\Delta}^{L}$) and right handed ($f_{\Delta}^{R}$) Yukawa couplings separately. The neutrino mass matrix in the $\left(\overline{\nu}_{L},n_{2R},\overline{n}_{2L}\right)$ basis is then given by:

[TABLE]

where each entry is a $3\times 3$ matrix with $m_{D}=f_{\zeta}\,v_{2}$, $m_{2}^{\prime}=f_{\Delta}^{R}\,u_{3}$, $m_{2}={f_{\Delta}^{L}}^{*}\,u_{3}$, and $M$ is a free Dirac mass term in $M\left(\overline{n}_{2L}n_{2R}+\overline{n}_{2R}n_{2L}\right)$. Thus, the inverse seesaw neutrino mass is given in the form (assuming $f_{\Delta}^{R}\simeq f_{\Delta}^{L}\simeq f_{\Delta}$):

[TABLE]

From here we see that if we assume $u_{3}\sim\mathcal{O}\left(100~{}\rm GeV\right)$ and $M\sim\mathcal{O}\left(1~{}\rm TeV\right)$, we can generate light neutrino mass in the correct ballpark for $f_{\Delta}\sim\mathcal{O}\left(10^{-6}\right)$ and $f_{\zeta}\sim\mathcal{O}(1)$. This is shown in the LHS of Fig. 1 where the two contours correspond to light neutrino mass $m_{\nu}\sim 0.1~{}\rm eV$ for smaller values of $f_{\Delta}:\{10^{-5},10^{-6}\}$ in blue and red respectively. One can however, choose a less fine-tuned $f_{\Delta}\sim\mathcal{O}\left(1\right)$ at the expense of making the RHNs super heavy $\sim 10^{5}$ GeV. This is depicted in the RHS of Fig. 1 where we have chosen a fixed $u_{3}=100~{}\rm GeV$ and obtained contours of correct neutrino mass ($m_{\nu}\sim 0.1~{}\rm eV$) for two different choices of $M:\{10^{4},10^{5}\}~{}\rm GeV$ (red and blue curves respectively) in the plane of $f_{\zeta}-f_{\Delta}$. We would like to mention here that the second choice of heavy $n_{1,2}$ is more desirable as it does not require the Yukawa couplings to be extremely fine-tuned, i.e. $\sim\mathcal{O}(10^{-6})$ and secondly will help us in addressing the freeze-in of $X(\bar{X})$ as the only possible decay mode of $\zeta_{1}^{\pm,0}$. More interestingly, it will also distinguish the collider signature of this model from that of the WIMP example [18] (details in Sec. 5).

4 Dark Sector

In this set-up we assume $X$ to be the lightest non-zero $S$ charge particle and hence a DM candidate as stated in Sec. 2. This choice is even more natural when $g_{N}$ is considered small as we have here for the freeze-in of $X$. Because of this, $X$ is not in thermal equilibrium in the early universe and is produced via the decay or annihilation of an odd-$S$ particle that can be in thermal bath. Therefore, in our model, $X$ can be produced via freeze-in from the decays of the scalar triplet $\Delta$ (Fig. 2) and the bi-doublet scalar components $\zeta_{1}^{0,\pm}$ (Fig. 3) when kinematically accessible. It is important to note that the decay occurs before and after the decoupling of the heavier particle from the thermal bath and we carefully illustrate how the ‘late decays’ can contribute significantly to the relic density of DM. To be compatible with the required relic abundance via freeze-in, the coupling requires to be $g_{N}\sim 10^{-15}-10^{-10}$, which in turn results in a very high mass scale for the $\Delta$’s as explained in Sec. 2. The mass hierarchy among different components of the $\Delta$ (namely, $\Delta_{1,2,3}$) is controlled by the parameter $f_{7}$. If we assume $f_{7}=0$ then all three components have the same mass. $\Delta$’s can be produced from the Higgs quartic interaction, which is a function of $f_{8}$. Therefore, freeze-in production of $\Delta$ requires $f_{8}\sim 10^{-12}$, while for freeze-out: $f_{8}\sim 1$. The Feynman graph for the production of $\Delta$ is shown in Fig. 4. For simplification, we assume $f_{9,10}\ll f_{8}$ so that the $\Delta\leftrightarrow\zeta$ mixing and conversions can be neglected. We will show that $\Delta$’s are naturally stable in the freeze-in scenario (for $g_{N}\sim 10^{-15}-10^{-10},~{}f_{8}\sim 10^{-12}$), and we have $\{\Delta_{1,2,3}\}$ as long-lived relics that contribute to the DM abundance in addition to $X(\bar{X})$. Different production possibilities of $\Delta$, depending on the choice of the model parameters and possible degeneracy is classified in Tab. 4.

We primarily focus on the degenerate scalar triplet scenario with $f_{7}=0$ and calculate the yield of the DM components ($\{X,\Delta_{1,2,3}\}$) via freeze-in by solving Boltzmann equations (BEQ) in subsection. 4.1.

4.1 Degenerate $\Delta$’s with $f_{8}\sim\mathcal{O}(10^{-12})$

The masses of particles in $\Delta$ triplet are degenerate if $f_{7}=0$. In this case, $\Delta_{1,2,3}$ can decay via the tree level diagram on the LHS of Fig. 2. This diagram is only possible after the EWSB when $\Delta_{3}$ mixes with the Higgs via $f_{8}\left(\Phi^{\dagger}\Phi\right)tr\left(\Delta^{\dagger}\Delta\right)$. The effective vertex for $\Delta_{3}\to t\bar{t}$ decay is $\sim f_{8}m_{t}u_{3}/m_{\Delta_{3}}^{2}$. Since $m_{\Delta_{3}}\sim\mathcal{O}(10^{13}~{}\rm GeV)$ and $f_{8}\sim\mathcal{O}(10^{-12})$ in the freeze-in scenario, the decay vertex factor is tiny. This results in a very large lifetime for $\Delta_{3}$. As an estimate, if we set $m_{\Delta_{3}}=10^{13}~{}\rm GeV$, $f_{8}=10^{-12}$ and $u_{3}=100~{}\rm GeV$, this decay width turns out to be $2\times 10^{-57}~{}\rm GeV$, which dubs into $\sim 6\times 10^{32}~{}\rm sec$ in terms of decay lifetime. This is of course much larger than lifetime of the universe, which is $\sim 10^{17}~{}\rm sec$. This makes all of the $\Delta$’s stable 666If we calculate the lifetime of $\Delta_{1}$ from the decay $\Delta_{1}\to\bar{X}\bar{X}t\bar{t}$ turns out to be $6\times 10^{34}~{}\rm sec$ for similar choices of $g_{N}=10^{-13}$, $f_{8}=10^{-12}$ and $m_{X}=5~{}\rm GeV$..

On the other hand, the loop-induced decay of $\Delta_{3}$ to SM neutrinos (RHS of Fig. 2), is possible even before the EWSB mediated by Yukawa couplings with $n_{2}$. The decay width for this process 777The amplitude for the loop diagram has been computed manually and cross-checked using Package-X [26]. is understandably small as it is proportional to the SM neutrino mass. Even if we assume the Yukawa couplings involved in this decay to be $f_{\zeta}\sim\mathcal{O}(1)$ and $f_{\Delta}\sim\mathcal{O}(10^{-6})$, the decay width turns out to be $\sim 4.75\times 10^{-65}~{}\rm GeV$ for $m_{n_{2}}\sim 1$ TeV and $m_{\zeta_{2}}\sim 200$ GeV. Therefore, with a lifetime of $10^{41}~{}\rm sec$, we conclude that $\Delta_{3}$ (and hence $\Delta_{1,2}$) is always a long-lived relic. $X$ is kinematically stable and the prime DM candidate of the model with $m_{X}<m_{\zeta_{1}},m_{\Delta_{2}}$. If we further assume $g_{N}\sim\mathcal{O}(10^{-13})$, then $\zeta_{1}^{0,-}\to\overline{X}+\zeta_{2}^{0,-}$ will produce $X$ via freeze-in mechanism (LHS of Fig. 3). Note that after spontaneous breakdown of $SU(2)_{N}$, bi-doublet $\zeta$ breaks into two doublets, i.e. $\left(\zeta_{1}^{0},\zeta_{1}^{-}\right)$ & $\left(\zeta_{2}^{0},\zeta_{2}^{-}\right)$. After the electroweak SSB, $\zeta_{2}^{0,-}$ components will decay as shown in the RHS of Fig. 3. As mentioned above, we are assuming $f_{9,10}\ll f_{8}$ for simplicity, such that $\Delta\leftrightarrow\zeta$ conversion is negligible. In the following subsection we calculate abundance for both $X$ and $\Delta$ using appropriate BEQ.

4.1.1 Computation of yield for $\Delta$ and $X$

Let us first estimate the yield of $\Delta$, which become stable due to small $f_{8}$ and contribute to DM relic density as elaborated in the last section. The rate of change of number density of $\Delta$ is governed by the following BEQ:

[TABLE]

where $d\Pi_{i}\equiv d^{3}p_{i}/(2\pi)^{3}2E_{i}$ is the phase space factor, and number density is given by

[TABLE]

where $g_{i}$ denotes the effective relativistic degrees of freedom and $H$ is the Hubble constant. Following Fig. 4, the only way of producing $\Delta$ or depleting its number density occurs through the quartic interaction with Higgs ($H$). We are focusing on the freeze-in production of $\Delta$ before EWSB and after SSB of $SU(2)_{N}$. Higgs can then produce a pair of superheavy $\Delta$ only beyond the threshold center-of-mass (C.O.M) energy $s\geq 4m_{\Delta}^{2}$.

Now, assuming negligible initial abundance for $\Delta_{i}$’s, we can set $f_{\Delta_{i}}=0$. We can also neglect the Pauli-blocking/stimulated emission effects, i.e. $f_{i}\ll 1$. Therefore, Eq. (4.1) simplifies to

[TABLE]

where, $T$ is the temperature, $s$ is the C.O.M energy of the production process $\Phi\Phi^{*}\to\Delta_{i}\Delta_{i}^{*}$ and we defined $P_{ij}$ as [27]:

[TABLE]

Since $|\mathcal{M}|^{2}_{\Phi\Phi^{*}\to\Delta_{i}\Delta_{i}^{*}}=f_{8}^{2}$, Eq. (4.3) simplifies to:

[TABLE]

where the number density is converted to comoving yield, scaled as $Y_{\Delta}=n_{\Delta}/\mathcal{S}$ and $\mathcal{S}=2\pi^{2}g_{\star}^{S}\,T^{3}/45$ is the entropy and $m_{\Phi}=0$. Using $\dot{T}\approx-T\cdot H$, and $x\equiv m_{\Delta}/T$ we rewrite Eq. (4.5) as:

[TABLE]

so the yield for $\Delta_{i}$ can be written as:

[TABLE]

where the last equality only holds if we assume a constant relativistic DOF $g_{\star}\approx 100$ during the freeze-in of $\Delta_{i}$’s. The total relic abundance of $\Delta$’s, i.e. $\Omega_{\Delta}=3\Omega_{\Delta_{i}}$ owing to its degeneracy is given by

[TABLE]

We see that for $f_{8}\in 6.19\times[10^{-13},10^{-12}]$, we get relic density in the correct ball park, when the corresponding relic density of $\Delta$ is $\Omega_{\Delta}\cdot h^{2}\approx[\%1,\%100]$ of the observed DM relic density.

We now focus on $X$. The relic abundance of $X$ is coming from the $\zeta_{1}^{0,-}$ decays (see Fig. 3). We also assume that $n_{1,2}$ is heavier than $\zeta_{1}$, so that $\zeta_{1}\to\zeta_{2}+X$ constitutes 100% decay branching fraction of $\zeta_{1}$. The contributions from the $\Delta$ decays in Fig. 2, are exceedingly slow for small $f_{8}$. This makes $\Delta$ a DM candidate (as explained before) with negligible contribution to the freeze-in production of $X$. The production of $X$ from annihilation ($\zeta_{1}\zeta_{1}\to XX$) is much more suppressed due to the presence of $g_{N}^{2}$ in the amplitude and can be neglected. Therefore, the contribution to the number density of $X$ can then be written as:

[TABLE]

where $\theta$ is the step function, $t_{D}$ is the $\zeta$ decoupling time from the SM particles, and the factor of $2$ arises from same contribution of $\zeta_{1}^{0}$ and $\zeta_{1}^{-}$ decays to $X$ production. It is important to note that when $\zeta_{1}$ is in thermal equilibrium with the hot plasma, $T=T^{\prime}$, but after $\zeta_{1}$ decouples: $T\neq T^{\prime}$. We have set $f_{\bar{X}}=0$ in accordance with freeze-in prescription and used:

[TABLE]

where $g_{\zeta_{1}}=1$. Since $\zeta$ is a cold relic then $E_{\zeta_{1}}\approx m_{\zeta_{1}}$ for $t>t_{D}$, which simplifies Eq. (4.9):

[TABLE]

where $\gamma_{\zeta_{1}}=E_{\zeta_{1}}/m_{\zeta_{1}}$. Using the comoving yield $Y_{\bar{X}}=n_{\bar{X}}/\mathcal{S}$ and $\dot{T}\approx-T\cdot H$ we have:

[TABLE]

The presence of $T^{{}^{\prime}}$ in BEQ segregates the equation to the cases (i) before decoupling and (ii) after decoupling of $\zeta_{1}$ from thermal bath. It is important to relate temperature $T^{{}^{\prime}}$ to the decoupling temperature ($T_{D}$). Throughout our calculations we use (see Appendix A.2):

[TABLE]

for cold relics. We now express $T^{\prime}$ in terms of the temperature of the hot plasma ($T$) using the conservation of entropy, i.e. $S=\mathcal{S}\cdot R^{3}(t)=g_{\star}^{S}\,T^{3}\,R^{3}(t)=\text{const.}$888$g_{\star}^{S}$ and $g_{\star}^{\rho}$ are the effective number of relativistic degrees of freedom for entropy and energy respectively.:

[TABLE]

Now, since at the decoupling $T^{\prime}_{D}=T_{D}$, we then have:

[TABLE]

Using the definitions of $\mathcal{S}$ and $H$ we get BEQ after decoupling as

[TABLE]

while for $t<t_{D}$ we have:

[TABLE]

Changing the variable from $T\to x=m_{\zeta_{1}}/T$ yields:

[TABLE]

where

[TABLE]

Using above, the BEQs before and after decoupling Eq.’s (4.16) &  (4.17) become:

[TABLE]

where we defined

[TABLE]

The total $\bar{X}$ yield from $\zeta$ decays is given by:

[TABLE]

The second integral within the parenthesis (from $x_{D}$ to $\infty$) indicates the contribution of the late decays of $\zeta_{1}$ (after its decoupling from thermal bath) to the production of DM ($X$). We would like to point out that an elaborate estimation of such contribution has not been carried out before in literature. The relic abundance of $X$ and $\bar{X}$ is then given by

[TABLE]

where the factor of two is due to the contribution from both $X$ and $\bar{X}$. Note that, Eq. 4.23 is practically independent of the model concerned (except for the couplings and masses appearing in $\Sigma_{\zeta}$) as it solely involves the number of degrees of freedom $g_{*}^{s}$ and $g_{*}^{\rho}$ and the dimensionless quantity $x_{D}$ which for all practical purposes can be taken to be $x_{D}\sim 25$. Hence this is a generic expression for computing yield in all such cases where the DM is produced via freeze-in from the decay of a heavy species before and after its decoupling from the thermal bath.

In Fig. 5 we have plotted the evolution of the relic abundances of $\Delta$ and $X$ as a function of the dimensionless parameter $x=M_{\zeta_{1}}/T$. The parameters in Fig. 5 are chosen such that $\Omega_{\Delta}+\Omega_{X}=\Omega_{\text{DM}}\cdot h^{2}=0.1198$. We have chosen four different values of $m_{\zeta}$, to generate four different relative contributions of $X$ and $\Delta$ to the total relic density as $\Omega_{X}/\Omega_{\text{DM}}=\{\%25,\%50,\%75,\%99\}$ from top-left to bottom-right respectively. The evolution of yield for $\Delta$ (shown by solid blue line in each of the plots) follows the familiar pattern of freeze-in production, which initially increases with $x$, and finally becomes constant as the temperature of the universe falls below the mass of $\Delta$. For $X$ (the solid magenta curve), note that the first point where the relic abundance function flattens is around x\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle>}}{{\sim}}~{}}}10, where the exponential suppression inside the Bessel function becomes dominant. However, the function then rises for a second time, and captures the contribution from the decay after decoupling of $\zeta$. First of all, the second rise in freeze-in curve for $X$ shows that the yield from decays after the decoupling can be as large as the yield before decoupling. This indicates the importance of the late decay to be considered for correct evaluation of freeze-in relic density. The effect of late decay can be understood in more details from Eq. 4.21. The share from the second integral ($x>x_{D}$) in Eq. 4.21 is at first suppressed by $K_{1}[x_{D}]$ and therefore it takes a long time (not until $x\sim 10^{4}$) for the yield to build up to an amount comparable to the yield before $x_{D}$. This happens because the asymptotic exponential suppression in $K_{1}(\alpha x^{2}/x_{D})$ cancels out the increasing exponential term $\exp(\alpha x^{2}/x_{D})$. This gain stops later as the exponential suppression from the decays becomes dominant around $x\sim 10\cdot m_{\zeta_{1}}/\sqrt{M_{\text{Pl}}\Gamma}\sim(5\cdot 10^{4}-10^{5})$999Assuming $m_{\zeta_{1}}\sim 300$ - $700$ GeV, $m_{\zeta_{2}}=200$ GeV, $m_{X}=5$ GeV, we have $\Gamma\sim 10^{-23}-10^{-21}$ GeV, where $g_{\star}^{\rho}\sim 10.75$. . The sole contribution of the late decays to the yield are separately plotted by dashed magenta curves in each plot in Fig. 5, to show the exact $x$ for the late rise of the yield. It is also clear that if the decay rate was much faster, the second integral would become suppressed and we wouldn’t see any effect of late production.

The behaviour of the relic abundance of $X$ (for two sample sets of values of $m_{X}$) as a function of $g_{N}$ is plotted in Fig. 6. We see that for small values of $g_{N}$, the abundance remains the same. This is because for these values of $g_{N}$ the decay rate is so small that the decays before the freeze-out of $\zeta_{1}$ can be neglected and the freeze-out yield of $\zeta_{1}$ is the same for different $g_{N}$ values. After $\zeta_{1}$ freezes out, all of them will eventually decay to $X$ and changing $g_{N}$ will only vary the time scale of these decays. However, beyond a threshold value of $g_{N}$ (characteristic to specific $m_{X}$), relic density rises sharply with $g_{N}$ as a function of $\sim g_{N}^{2}$ as the yield in is proportional to the decay width in case of freeze-in. This is in contrast to the freeze-out scenario where relic density is inversely proportional to annihilation cross-section. As we explain in the next section, smaller values of $g_{N}$ are ruled out by the constraints from the BBN (indicated by dashed curves in Fig. 6).

4.1.2 Bounds on decaying relic particles from BBN and CMB

Standard BBN may be significantly perturbed by the energy injections due to both neutral decays [14] $\zeta_{1}^{0}\to\overline{X}+\zeta_{2}^{0}$ and charged decays [15] $\zeta_{1}^{-}\to\overline{X}+\zeta_{2}^{-}$. In the case of decaying electrically charged particles with $\tau>100$ seconds, the existence of bound states between nuclei and the relics may significantly change nuclear reaction rates. Here, we analyze the decays with $0.05$ s $<\tau<100$ s using the results from Ref. [14], and rule out the ones with $\tau>100$ s for simplicity. As a result, this puts a lower bound on $g_{N}$, depending on different choices of DM masses. The hadronic branching ratio ($B_{h}$) for $\zeta_{1}^{0}$ decays is $B_{h}^{1}\approx 0.695$, and for $\zeta_{1}^{-}$ depends on parameters like $M_{n_{2}},f_{\zeta}$. We can assume for heavy $n_{2}$ and/or small $f_{\zeta}$, the weak decays $\zeta_{1}^{-}\to\bar{X}+\zeta_{2}^{-}\to\bar{X}+\zeta_{2}^{0}+W^{-}$ become dominant and so the hadronic branching ratio of the charged component is $B_{h}^{2}\approx 0.9$. Therefore, for a conservative bound we assume hadronic branching for the both cases to be $B_{h}=1$ throughout our analysis. We also assume that the mass difference $(m_{\zeta_{1}}^{-}-m_{\zeta_{1}}^{0})$ generated during the EWSB is negligible compared to their masses before EWSB. We then extrapolate the results in [14] for arbitrary $m_{\zeta_{1}}$. The total energy of the resulting quarks in these decays is taken to be the same as $m_{\zeta_{1}}$ following [14]. We use the relation in Eq. 4.25 to convert the bounds on the relic abundance ($\Omega^{\text{BBN}}_{\zeta_{1}}$) of $\zeta_{1}$ (assuming it was stable) to bounds on the DM ($X$) relic abundance (from $\zeta_{1}$), i.e. $\Omega_{\text{DM}}^{X}$. Given that $n_{\zeta_{1}}=n_{X}$, and $\Omega_{i}\propto m_{i}\cdot n_{i}$ we have:

[TABLE]

where $\Omega_{\zeta_{1}}^{\rm BBN}$ corresponds to the abundance of the decaying species ($\zeta_{1}$ in this case), which is fixed from BBN depending on its hadronic branching ratio and lifetime (or in other words $g_{N}$) [14]. The dashed portions in Fig. 6 are thus ruled out by the BBN bound since smaller values of $g_{N}$ are discarded in order to prevent prolonged decays that may disturb the standard BBN mechanism.

A decaying long-lived DM candidate can be constrained by various observations. For example, it can alter the ionization and heating history of the CMB and its power spectrum [16]. Decays to several decay modes e.g. $b\bar{b}$, $WW$, $\mu\bar{\mu}$ can also be constrained by the AMS-02 precise measurements of the antiproton/proton ($p/\bar{p}$ ) fraction [17], as no evidence of new source of antiproton has been found in these data. This in turn, results in the following bounds on the hadronic decay life time of the DM [16]:

[TABLE]

As we already have shown, the $\Delta$’s have lifetime comparable to that of the universe. Now, if we set $m_{X}=5~{}\rm GeV$, $g_{N}=10^{-12}$ and $f_{8}=10^{-11}$ then for $u_{3}=100~{}\rm GeV$ we obtain the following decay width and corresponding decay lifetime for the $\Delta$’s:

[TABLE]

All these decay rates turn out to be much longer than the bounds mentioned in Eq. 4.26, and so only the bounds from the $\zeta_{1}$ decays are relevant for final estimate of allowed parameter space. Note that, since some of the energy of $\Delta_{1}$ goes into $\bar{X}$’s, the actual bounds are smaller, and these numbers are conservative.

4.1.3 Summary of available parameter space

The relic density allowed parameter space for $X$ is plotted in $m_{\zeta_{1}}-m_{X}$ plane in Fig. 7 for different choices of $g_{N}=\{10^{-12},10^{-13},10^{-14},10^{-15}\}$ (clockwise from top left). Each plot contains several contours corresponding to a specific percentage of the total DM density coming from $X$ (numbers inside parenthesis). Therefore, each curve corresponds to a specific value of $f_{8}$ such that the rest of the DM abundance ($100\%-$ the percentage in the parenthesis) is coming from $\Delta$. This is possible because the relic density of $X$ and that of $\Delta$ are uncorrelated in the model. We see from the plots that constant $\Omega_{X}$ curves rises in $m_{\zeta_{1}}-m_{X}$ plane for larger $m_{X}$ and then takes a sharp turn to decrease the required $m_{\zeta_{1}}$ to satisfy relic density and then rises again. In order to interpret the shape of the contours in Fig. 7, let’s first take a look at the parameter space assuming the late decays are ignored. This is shown in Fig. 8 with the dashed upward moving faint lines in the background. These contours basically correspond to constant relic density ($\Omega_{X}h^{2}$), which is mainly dictated by $\Omega_{X}\propto c\approx m_{X}\cdot\Sigma_{\zeta}$, since the second integral in Eq. (4.23) is to be ignored. Therefore, it is straightforward to see that the larger is the value of $m_{X}$, the smaller $\Sigma_{\zeta}$ has to be and therefore $m_{\zeta}$ will needs to be larger as $\Sigma_{\zeta}\propto 1/m_{\zeta}^{2}$. It is also clear, for small $m_{X}$ (as we require larger $\Sigma_{\zeta}$) the second integral in Eq. (4.23) has smaller contribution due to exponential suppression from $\Sigma_{\zeta}$. This is reflected in Fig. 8 where the case with late decays and the case without them merge together at low $m_{X}$. As we trace these contours from low $m_{X}$ region (left) to high $m_{X}$ regions (right), we observe a departure for the curves where the late decays are involved. Now, with larger $m_{X}$, $\Sigma_{\zeta}$ decreases. This increases the contribution from the late decays as the second integral in Eq. (4.23) is proportional to $e^{-\Sigma_{\zeta}}$. Therefore, if we compare the case with late decay ($\Sigma_{\zeta}^{1}$) to that without late decay ($\Sigma_{\zeta}^{2}$), it is easy to see that:

[TABLE]

which is reflected in the drop of the curve beyond certain $m_{X}$.

We can estimate the value of $m_{X}$ for which this change of behaviour happens. As mentioned at the end of subsection 4.1.1, the exponential suppression in the second integral becomes notable around $x_{e}\sim m_{\zeta_{1}}/\sqrt{M_{\rm Pl}\Gamma}=1/\sqrt{\Sigma_{\zeta}}$. Since the yield from the second integral becomes relevant around $x_{obs}\sim 10^{4}$ (see Fig.7), we can use x_{e}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}x_{obs} as a criterion for this change. Let’s start from low $m_{X}$ region and trace the constant $c$ contours for which $\Omega_{X}\cdot h^{2}\approx 0.1198$. Using Eq. 4.23 and 4.24 we obtain $c\approx 10^{-7}$ for the cyan dashed curve where 100 % contribution to the relic abundance comes from $X$. Now, since $x_{e}\approx\sqrt{M_{X}/c}$, this criterion transmits to M_{X}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10 GeV. If we choose the curve for $5\%$ yield instead, we get $c\sim 10^{-8}$ and M_{X}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}1 GeV. This explains the sharp turn in all of the curves that happens around $M_{X}\approx 1-10$ GeV. Also note that as $m_{\zeta_{1}}$ is lowered and $m_{X}$ is increased, we get into a region where the decay rate is so slow that it violates the BBN bounds (dashed region). As $m_{X}$ is increased even further, all the contours converge to the kinematically forbidden boundary i.e. $M_{\zeta_{1}}=M_{X}+M_{\zeta_{2}}$. This characteristic feature is reflected in each plot in Fig. 7 for different choices of $g_{N}$ with $M_{\zeta_{2}}$ fixed at 200 GeV.

4.2 Degenerate $\Delta$’s with $2m_{X}<m_{\Delta}$ and $f_{8}\sim\mathcal{O}(1)$

In this case let us first compute the lifetime of the $\Delta$’s assuming $f_{8}\sim\mathcal{O}(1)$. We once again recall here that $m_{\Delta}\sim\mathcal{O}(10^{13})~{}\rm GeV$. Let us choose $m_{X}=5~{}\rm GeV$ and $g_{N}\sim\mathcal{O}(10^{-13})$ for which $X$ yields relic density in correct ballpark. We then find, $\Gamma_{\Delta_{3}}\to t\bar{t}=1.3\times 10^{-35}~{}\rm GeV$ which dubs into a lifetime of $\sim 5\times 10^{11}~{}\rm sec\ll\tau_{\text{universe}}$. With the same set of parameters, we find $\Gamma_{\Delta_{2}\to\overline{X}t\overline{t}}=5.5\times 10^{-36}~{}\rm GeV$, which gives rise to a lifetime of $\sim 10^{11}~{}\rm sec$. The decay width for $\Delta_{1}$ turns out ot be $\Gamma_{\Delta_{1}\to\overline{X}\overline{X}t\bar{t}}=4.5\times 10^{-37}~{}\rm GeV$, with a lifetime of $\sim 10^{12}~{}\rm sec$. Therefore, with $f_{8}\sim\mathcal{O}(1)$ we end up with a situation where all of the triplet scalar components are unstable and don’t contribute to the DM relic abundance. Therefore we end up with a single component $X$ DM set up. The initial freeze-out abundance of $\Delta$’s (before they decay to $X$) is determined by

[TABLE]

which yields a freeze-out abundance for $\Delta$ at the order of:

[TABLE]

If we use the fact that the number of $X$’s produced from the decays of $\Delta$, i.e. $n_{X}^{\rm FI}\sim n_{\Delta}^{\rm FO}$, we then see that $\Omega_{X}^{\rm FI}=(m_{X}/m_{\Delta})\Omega_{\Delta}^{\rm FO}$ is about $6$ orders of magnitude higher than the correct DM abundance, which rules out this scenario with $f_{8}\sim\mathcal{O}(1)$.

5 Collider search

Possible collider signatures of this model in the context of usual freeze-out of $X$ was elaborated in [18]. It was pointed out that, the only field connecting the dark sector (having $SU(2)_{N}$ charge) with the SM is the scalar bi-doublet. Therefore, the charged and neutral components of the scalar bi-doublet can be produced at the LHC via the diagrams shown in the top panel of Fig. 9. We should also note that the other scalars, namely the $SU(2)_{N}$ doublet ($\chi$) and triplet ($\Delta$) possess Higgs portal interactions. However, we have assumed all those portal couplings to be small to satisfy the Higgs data (see Sec. 2). In particular, the Higgs portal coupling ($f_{8}$) between the triplet ($\Delta$) and SM doublet ($\phi$) was assumed to be small to address correct freeze-in production of the triplet components, which serve as DM. Therefore, the other particles do not have any chance of being produced at the LHC. The variation of production cross-section of the charged and neutral components of the scalar bi-doublet at the LHC with respect to its mass $m_{\zeta_{1}}$ for $E_{\text{CM}}=$14 TeV is shown in the bottom panel of Fig. 9. It should be noted that charge current contribution is more pronounced than the neutral current contribution which we pointed out in our earlier analysis [18]. We have implemented the model in CalcHEP [28] and used CTEQ6l [29] as a representative parton distribution function for generating this process. The SM gauge coupling is mostly responsible for the production of the bi-doublet components at the LHC. The Higgs portal interactions are again assumed to be small to be compatible with the Higgs data.

If the additional neutrinos ($n_{1},n_{2}$) are lighter than $\zeta_{1}^{\pm,0}$ (which is quite legitimate in WIMP scenario even after addressing correct neutrino mass generation), then these bi-doublet scalars can further decay to neutrinos and SM leptons ( $\zeta_{1}^{+}\to\ell^{+}+n_{1R}$) via Yukawa interaction given in Eq. 3.2, giving rise to the following signatures in colliders [18]:

• •

Single lepton with missing energy ($1\ell^{\pm}+\not{E}_{T}$) due to charged current interaction.

• •

Opposite sign di-lepton with missing energy $(\ell^{+}\ell^{-}+\not{E}_{T})$ due to neutral current interaction.

In the WIMP-like freeze-out scenario of $X$, therefore this model may leave an imprint of leptonic signal excess at the colliders. However, in the present framework, we assume that the DM production is from the decay of $\zeta_{1}\to\zeta_{2}+X$, and the additional neutrinos ($n_{1},n_{2}$) to be heavier than $\zeta_{1}$. It is understandable that the DM analysis will be modified by the corresponding decay branching ratio if $n_{1,2}$ are lighter than $\zeta_{1}$. It is also instructive to remind that heavier $n_{1,2}$ ($\sim 10^{6}$ GeV) arises naturally when we assume the Yukawa in a legitimate ball park $f_{\Delta}\sim 0.1$ via a large $u_{3}\sim 100$ GeV to address the required neutrino mass (see Sec. 3). However, one may still assume the presence of lighter $n_{1,2}$ ($\sim$ TeV) with fine tuned $f_{\Delta}\sim 10^{-7}$ and $m_{n_{1R}}<m_{\zeta_{1}}$, yielding single lepton and di-lepton signature as before. However in that case, the DM production will be further suppressed and appropriate decay branching ratios have to be assumed.

Therefore, in the light of DM analysis performed here, the neutrinos ($n_{1,2}$) are heavier than the scalar bi-doublet $\zeta_{1,2}^{\pm},\zeta_{1,2}^{0}$. As a result, these scalars can not decay to RHN plus SM leptons. Due to radiative correction, there should be a small mass splitting between the charged and neutral components of the bi-doublet: $\delta m$. As a result, the charged scalars, once produced, can undergo the decays shown in Fig. 10 via off-shell $W$ or via off-shell $n_{1,2R}$ to $\zeta_{1}^{0}+\ell+\nu_{\ell}$.

Note that the diagram in the LHS of Fig. 10 involves the vertex $\zeta_{1}^{+}\zeta_{1}^{0}W$, which is proportional to $g_{L}(p_{1}+p_{2})_{\mu}$ where $p_{1}$ and $p_{2}$ are the momenta of the incoming and outgoing scalars, while the other vertex involves only the SM gauge coupling $g_{L}$. The diagram on the RHS of Fig. 10, on the other hand, only depends on the Yukawa coupling $f_{\zeta}$, which we assume to be $\sim\mathcal{O}(1)$. Since $n_{1R}$ is heavy, these decays are dominated by the $W$-mediated process as in the LHS of Fig. 10. In Tab. 5, we show a couple of sample points with two different mass splitting, where 99% of the branching is carried away by the $W$-mediated decay as $m_{W}\ll m_{n_{1R}}$. Also note that in the limit $\delta m\ll m_{\zeta_{1}}$, the decay width mostly depends on $\delta m$ (which controls the phase space) and not on $m_{\zeta_{1}^{\pm}}$. For example, if we fix $\delta m=0.1~{}\rm GeV$, then for $m_{\zeta_{1}^{\pm}}=1.2~{}\rm TeV$, the decay width is $3.127\times 10^{-18}$ GeV, which changes to $3.126\times 10^{-18}$ GeV for $m_{\zeta_{1}^{\pm}}=1.5~{}\rm TeV$. For $\delta m\sim\mathcal{O}(100~{}\rm MeV)$ this model can give rise to charge track and/or displaced vertex [30, 8] at the colliders that can be probed by current or future experiments [31, 32]. A heavy stable charged particle (HSCP) such as $\zeta_{1}^{+}$ in our model, will typically travel with a velocity $\beta\equiv\frac{v}{c}<1$. Hence, as it passes through the silicon detectors, it produces an ionizing track with higher ionization energy loss rate (larger $dE/dx$) compared to the SM particles [33]. Also, if the HSCP decays outside the detectors, the time of flight (TOF) measured by the muon system will be longer than that of relativistic muons. These two features can distinguish non-standard HSCPs from the SM particles.

Typically, for c\tau\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}\mathcal{O}(10~{}\rm m) the searches for HSCPs are done via displaced vertex signatures, while for $c\tau>\mathcal{O}(10~{}\rm m)$ a sizable fraction will decay only after crossing the tracker and/or muon chamber [33]. We use the CMS bound on HSCP production cross-section at the LHC [34] to constrain the mass of $\zeta_{1}$ in our model. The $\zeta_{1}$ production cross-section is plotted in the LHS of Fig. 11 with respect to $\zeta_{1}$ mass. We have also shown the limit from the CMS tracker+TOF analysis [34] (for luminosity of $12.9~{}\rm fb^{-1}$) with a solid red curve. As one can see, these constraints rule out m_{\zeta_{1}^{\pm}}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}1~{}\rm TeV. We have tabulated the decay lengths for some selected values of $\delta m$ with $m_{\zeta_{1}^{\pm}}=1.2~{}\rm TeV$ in Tab. 6. The same is also shown in RHS of Fig. 11 where we have plotted the variation of the decay length with respect to mass splitting $\delta m$ for $m_{\zeta_{1}^{\pm}}=1.2~{}\rm TeV$. As mentioned, the decay length decreases as the mass splitting $\delta m$ increases. We have indicated a red dashed line, above which ($c\tau>10~{}\rm m$) the model can give rise to stable charged track and below (c\tau\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10~{}\rm m) displaced vertex signature.

6 Summary and Conclusions

In this draft, we have analysed FIMP realization of a non-abelian vector boson DM in $SU(2)_{N}$ extension of the SM. Non-abelian cases are important for several reasons, one because they require non-minimal extensions in the scalar sector for spontaneous symmetry breaking of the additional $SU(2)_{N}$ and therefore serve as an important framework to elaborate on Higgs physics in light of the present data. The model at hand also addresses neutrino mass generation and therefore neutrino mass constraint plays an important role in identifying the allowed parameters of the model together with DM constraints. For example, this exercise has led us to conclude (i) scalar triplets are super heavy (of the scale of $SU(2)_{N}$ breaking), (ii) sterile neutrinos assumed in the model can also be naturally heavier than scalar bi-doublet. Now, both of these two conclusions have immense phenomenological consequence. Therefore, the exercise performed in this analysis can serve as a benchmark to address non-thermal DM production together with neutrino mass and Higgs phenomenology.

One of the important outcomes of this analysis turns out to be out-of-equilibrium decay of a heavier particle to DM. In the present context, scalar bi-doublet $\zeta_{1}^{0,\pm}$ decays into its lighter partner $\zeta_{2}^{0,\pm}$ plus DM $X$, and this is solely responsible for non-thermal production of the DM. The production of DM occurs after $SU(2)_{N}$ symmetry breaking ($\sim 10^{12}$ GeV) and before electroweak symmetry breaking. We find that, the decay of $\zeta_{1}^{0,\pm}$ provides a significant contribution to DM relic density even after the freeze-out of $\zeta_{1}^{0,\pm}$ from thermal bath. This results in a sharp deviation of relic density contour in $M_{\zeta_{1}}-M_{X}$ allowed plane compared to the case where the late productions are neglected. The impact of this conclusion can be made in a generic and model independent way, to demand that any particle in thermal bath whose decay is slow enough ($\sim 10^{-24}$ GeV) can contribute significantly after freeze out and alter the available parameter space to a significant extent. We provide with a generic expression for the DM yield including the late decays that may serve useful in identifying such contributions for any model.

It is also important to note the connection between the dark sector and neutrino sector addressed in this model. The requirement of having a freeze-in vector boson DM makes the $SU(2)_{N}$ scalar triplet ($\Delta$), assumed for neutrino mass generation through inverse seesaw, superheavy ($\sim 10^{12}$ GeV). The decay modes of the neutral component of the triplet (to $t\bar{t}$ or to $\nu\bar{\nu}$) turns out to be extremely small, thanks to small $g_{N}~{}\rm{and}~{}f_{8}$ couplings. Therefore they are stable and serve as additional DM components in the model. It is intriguing to note that the correct relic density (or under abundance) for $\Delta$ can only be addressed if they are also produced non-thermally through Higgs quartic interaction. CMB data constrains the decay life time of long lived DM particles to hadronic final states ($\Delta\to b\bar{b}$ in our case) to be greater than $\tau(\text{DM}\to b\bar{b})\gtrapprox\,10^{24}\,s$. However, since the correct relic density is achieved for $\tau_{\Delta}\gtrapprox 10^{36}\,s$, this bound does not affect that part of the parameter space where $\Delta$ is a viable DM candidate.

There are other constraints as well. For example, bounds from AMS-02 constrains life time of hadronically decaying DM (again $\Delta\to b\bar{b}$ in our case) while BBN data constrains life time of semi-stable hadronically decaying charged and neutral particle ($\zeta_{1,2}^{0,\pm}$ in our case). The last bound crucially tames down a large allowed parameter space of our model by ruling out DM masses above $\sim 50$ GeV and g_{N}\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}10^{-14}. Constraints from CMB and AMS-02 turn out to be less sensitive due to the very long life time of scalar triplet DM.

The same model has been studied for WIMP realization as mentioned before. It is therefore important to identify the difference in their phenomenological implication. While freeze-in makes the DM insensitive to direct search, the WIMP can be detected via future direct search experiments. It is important to note that WIMP-like $X$ is allowed upto $\sim$TeV, but FIMP realization restricts it within $\sim 50$ GeV. The most crucial distinction however may arise from collider searches. While the WIMP realization could provide a signal excess in single or opposite sign di-lepton events associated with large missing energy, the FIMP case predicts stable charge track or displaced vertex signature, thanks to the production of scalar bi-doublet in the model. The decay of $\zeta_{1}^{\pm}\to\zeta_{1}^{0}$ here is restricted by the mass splitting of the order of $\delta m\sim 100$ MeV due to loop corrections. If \delta m\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}0.1 GeV, then the decay can lead to a stable charge track, while for \delta m\mbox{\raisebox{-2.58334pt}{~{}\stackrel{{\scriptstyle<}}{{\sim}}~{}}}0.5 GeV, we may see displaced vertex signature. On the contrary, in the WIMP realization, $\zeta_{1}^{\pm}$ can easily decay to $n_{1R}$ thanks to the Yukawa coupling (which is unlikely in FIMP realization due to a heavier $n_{1R}$) and serves as an interesting phenomenological consequence of the model.

Acknowledgements

M. Z. is supported by the CAS President’s International Fellowship Initiative (PIFI) grant #2019PM0110. BB and SB acknowledges DST-INSPIRE faculty grant IFA-13 PH-57. BB would also like to thank Sreemanti Chakrabarti, Rashidul Islam and Purusottam Ghosh for useful discussions. The authors would like to thank Ernest Ma for helpful discussions.

Appendix A Appendix

A.1 Scalar States

Using the minimization conditions in Eq. 2.5, we find the following massless Goldstone bosons:

[TABLE]

where the orthogonal state to A.1 and A.3 yields a heavy charged scalar field ($H^{+}$), and a physical complex scalar field ($\xi_{1}^{0}$) respectively. The orthogonal state to A.2 & A.4 yields another physical real scalar field ($\eta^{0}$).

In the Higgs sector, the $4\times 4$ mixing matrix spanning $\left(\text{Re}(\phi_{2}),\text{Re}(\zeta_{2}^{0}),\text{Re}(\chi_{2}),\text{Re}(\Delta_{3})\right)$ is given by

[TABLE]

which approximately yields a $2\times 2$ block spanning $\left(\text{Re}(\phi_{2}),\text{Re}(\zeta_{2}^{0})\right)$ given by:

[TABLE]

The resulting mass eigenstates are given in Sec. 2. We have also used the fact that the mixing term for $\left(\text{Im}(\chi_{2}),\text{Im}(\Delta_{3})\right)$ is negligible compared to $\Delta_{3}$’s mass term, and so $\Delta_{3}$ remains complex, with $m_{\Delta_{3}}^{2}\approx\,\,2u_{2}^{2}\,\mu_{2}/u_{3}$.

A.2 Evolution of chemical potential

Let $d^{3}p~{}d^{3}x$ be the phase space element at temperature $T$ and $d^{3}p^{{}^{\prime}}~{}d^{3}x^{{}^{\prime}}$ is that at temperature $T_{D}$. Since the distance scales as $R$ and the momentum (of a free particle) scales as $R^{-1}$, we can write:

[TABLE]

Since $N=n.R^{3}$ is fixed we have

[TABLE]

We consider two cases:

• (a)

Hot Relic: If particles are relativistic at the time of decoupling ($E\simeq p$), then:

[TABLE]

Hence, for hot relics: $T_{D}=\frac{R(t)}{R(t_{D})}.T$ and $T_{D}=\frac{T_{D}}{T}.\mu$.

• (b)

Cold Relic: If the particles are non-relativistic at the time of decoupling, then: $E\simeq m+\frac{p^{2}}{2m}$.

[TABLE]

Comparing $T_{D}=T(\frac{R}{R_{D}})^{2}$ with $m-\mu=(\frac{R_{D}}{R})^{2}.(m-\mu_{D})$ implies $\mu(t)=m+\left(\mu_{D}-m\right)\frac{T}{T_{D}}$.

A.3 Decay Rate for $\zeta_{1}\to\zeta_{2}+X$

Using

[TABLE]

the amplitude squared can be written as:

[TABLE]

where we used

[TABLE]

and so we have

[TABLE]

therefore the decay rate is given by

[TABLE]

A.4 Decoupling time for $\zeta_{1}\zeta_{1}^{*}\leftrightarrow$ SM

For the cold relic $\zeta_{1}$, the decoupling time $x_{D}$ can be determined using [35]:

[TABLE]

where $\Lambda$ is given by:

[TABLE]

and we assume the decoupling occurs before the EWSB so this cross section can be written as:

[TABLE]

where we use:

[TABLE]

Since $n$ is very heavy the gauge interaction dominates the cross section.

[TABLE]

[TABLE]