Flavoured gauge extension of singlet-doublet fermionic dark matter: neutrino mass, high scale validity and collider signatures

TL;DR

This paper introduces a gauged singlet-doublet fermionic dark matter model with an extended B-3L_tau symmetry, addressing neutrino masses, collider signatures, and high-scale stability, without ad-hoc symmetries.

Contribution

It proposes a novel gauged fermionic dark matter model with anomaly cancellation and neutrino mass generation, enhancing detection prospects and high-scale validity.

Findings

Model achieves anomaly cancellation with right-handed neutrinos.

Dark matter stability is naturally ensured by gauge symmetry.

Model remains valid up to the Planck scale with stable electroweak vacuum.

Abstract

We propose an Abelian gauged version of the singlet-doublet fermionic dark matter (DM) model where DM, combination of a vector like fermion doublet and a fermion singlet, is naturally stabilised by the gauge symmetry without requiring any ad-hoc discrete symmetries. In order to have good detection prospects at collider experiments like the large hadron collider (LHC) and enlarged parameter space for low mass DM, we consider the additional gauge symmetry to be based on the quantum $B-3L_{\tau}$ where the restriction to third generation of leptons is chosen to have weaker bounds from the LHC on the corresponding gauge boson. The triangle anomalies arising in this model can be cancelled by including a right handed neutrino which takes part in generating light neutrino masses through type I seesaw mechanism. Apart from DM, collider prospects and light neutrino masses, the model also offers high scale validity giving rise to a stable electroweak vacuum and perturbative couplings all the way up to the Planck scale. We constrain our model parameters from these requirements as well as existing relevant constraints related to DM and colliders.

Tables

Table 1. Table 1: Relevant particle content of the model and their charges under S ​ M × U ​ ( 1 ) B − 3 ​ L τ 𝑆 𝑀 𝑈 subscript 1 𝐵 3 subscript 𝐿 𝜏 SM\times U(1)_{B-3L_{\tau}} .

Particles	$SU{\left(3\right)}_c$	$SU\left(2\right)$	$U{\left(1\right)}_Y$	$U{\left(1\right)}_{B-3L_{\tau }}$
$q_L^i$	3	2	1/6	1/3
$u_R^i$	3	1	2/3	1/3
$d_R^i$	3	1	-1/3	1/3
$l_{1,2L}$	1	2	-1/2	0
$l_{3L}$	1	2	-1/2	-3
$e_{1,2R}$	1	1	-1	0
$e_{3R}$	1	1	-1	-3
$H$	1	2	-1/2	0
$𝒮$	1	1	0	3
$\Phi$	1	1	0	-3/2
$N_{R_{1,2}}$	1	1	0	0
$N_{R_3}$	1	1	0	-3
$\chi$	1	1	0	3/4
${\psi }^T:({\psi }^0,{\psi }^{-})$	1	2	-1/2	3/4

Table 2. Table 2: Choices of the benchmark points used for collider analysis. Masses, mixings, relic density and direct search cross-sections for the DM candidate are tabulated. In each case corresponding mass of Z B − 3 ​ L τ subscript 𝑍 𝐵 3 subscript 𝐿 𝜏 Z_{B-3L_{\tau}} is also quoted.

Benchmark	$\tilde{v}$	$g_{B-3L_{\tau }}$	$\sin \theta$	$\Delta M$	$M_{{\psi }_1}$	$\left|Y\right|$	${10}^3\hat{S}$	$\Delta T$	${\sigma }_{DD}$	$\Omega h^2$	$m_{Z_{B-3L_{\tau }}}$
Point	($v_{𝒮}/v_{\Phi }$)			(GeV)	(GeV)				$\left(c{m}^2\right)$		(TeV)
BP1	1.75	0.20	0.28	194.9	128.5	0.301	0.43	0.005	${10}^{-45.96}$	0.122	2.57
BP2	2.33	0.16	0.11	517.4	55.4	0.311	0.06	0.004	${10}^{-46.71}$	0.119	2.73
BP3	3.50	0.11	0.46	130.1	300.5	0.305	1.30	0.004	${10}^{-45.59}$	0.121	2.79
BP4	3.63	0.10	0.14	274.0	36.5	0.218	0.03	0.003	${10}^{-46.78}$	0.120	2.61
BP5	5.49	0.07	0.42	111.4	245.5	0.248	1.23	0.003	${10}^{-45.81}$	0.119	2.68

Table 3. Table 3: Production cross-section of charged VLF pairs for the chosen BPs in table  2 in presence of Z B − 3 ​ L τ subscript 𝑍 𝐵 3 subscript 𝐿 𝜏 Z_{B-3L_{\tau}} (left column) and in presence of only SM Z 𝑍 Z (right column).

Benchmark	${\sigma }_{pp\to {\psi }^{+}{\psi }^{-}}$	${\sigma }_{pp\to {\psi }^{+}{\psi }^{-}}$
Point	($Z+Z_{B-3L_{\tau }}$)	(Only $Z$)
	(fb)	(fb)
BP1	45.27	44.72
BP2	3.82	3.71
BP3	15.78	15.69
BP4	46.76	46.62
BP5	32.12	32.02

Table 4. Table 4: Variation of final state signal cross-section with MET cut for a fixed cut on H T > 250 ​ GeV subscript 𝐻 𝑇 250 GeV H_{T}>250\leavevmode\nobreak\ \rm GeV . All simulations are done at s = 14 ​ TeV 𝑠 14 TeV \sqrt{s}=14\leavevmode\nobreak\ \rm TeV .

Benchmark	${\mathrm{E̸}}_T$		${\sigma }^{\mathrm{OSD}}$
Points	(GeV)		(fb)
BP1	$>100$		0.82
	$>200$		0.21
	$>300$		0.06
BP2	$>100$		0.53
	$>200$		0.24
	$>300$		0.08
BP3	$>100$	$H_T>$ 250 GeV	0.04
	$>200$		0.006
	$>300$		0.001
BP4	$>100$		1.45
	$>200$		0.38
	$>300$		0.12
BP5	$>100$		0.05
	$>200$		0.01
	$>300$		0.002

Table 5. Table 5: Variation of final state SM background cross-section with MET cut for a fixed cut on H T > 250 ​ GeV subscript 𝐻 𝑇 250 GeV H_{T}>250\leavevmode\nobreak\ \rm GeV . All simulations are done at s = 14 ​ TeV 𝑠 14 TeV \sqrt{s}=14\leavevmode\nobreak\ \rm TeV .

Processes	${\sigma }_{\mathrm{production}}$	${\mathrm{E̸}}_T$		${\sigma }^{\mathrm{OSD}}$
	(pb)	(GeV)		(fb)
$t\overline{t}$		$>100$		0
	814.64	$>200$		0
		$>300$		0
$W^{+}{W}^{-}$		$>100$		5.99
	99.98	$>200$		0.99
		$>300$		0
$W^{+}{W}^{-}Z$		$>100$	$H_T>$ 250 GeV	0.05
	0.15	$>200$		0.02
		$>300$		0.009
$ZZ$		$>100$
	14.01	$>200$		0
		$>300$		0

Table 6. Table 6: Invisible Higgs branching ratio and invisible Z 𝑍 Z decay width for different benchmark points tabulated in table  2 . “-” stands for cases where M ψ 1 > m Z / 2 subscript 𝑀 subscript 𝜓 1 subscript 𝑚 𝑍 2 M_{\psi_{1}}>m_{Z}/2 and/or > m h 1 / 2 absent subscript 𝑚 subscript ℎ 1 2 >m_{h_{1}}/2 .

Benchmark	$B{r}_{inv}^{higgs}$	${\Gamma }_{inv}^Z$ (MeV)
Point		(MeV)
BP1	-	-
BP2	0.019	-
BP3	-	-
BP4	0.072	0.073
BP5	-	-

Equations

$$\sum _{LH}\mathrm{Tr}\left[T^a\{T^b,T^c\}\right]-\sum _{RH}\mathrm{Tr}\left[T^a\{T^b,T^c\}\right],$$

$$\displaystyle\begin{split}{\text{gravity gauge}\implies}9(-x_{d}+2x_{q}-x_{u})+(-x_{e_{R}}+2x_{\ell}-x_{N_{R_{3}}})=\leavevmode\nobreak\ &0&\\ {\text{$U(1)_{B-3L_{\tau}}^{3}$}\implies}9\left(-x_{d}^{3}+2x_{q}^{3}-x_{u}^{3}\right)+\left(-x_{e_{R}}^{3}+2x_{\ell}^{3}-x_{N_{R_{3}}}^{3}\right)=\leavevmode\nobreak\ &0&\\ {\text{$U(1)_{B-3L_{\tau}}^{2}U(1)_{Y}\implies$}}9\left(-\frac{1}{3}(-1)x_{d}^{2}+\frac{2x_{q}^{2}}{6}-\frac{2x_{u}^{2}}{3}\right)+\left((-1)x_{\ell}^{2}-(-1)x_{e_{R}}^{2}\right)=\leavevmode\nobreak\ &0&\\ {\text{$U(1)_{Y}^{2}U(1)_{B-3L_{\tau}}\implies$}}9\left(\left(-\frac{1}{3}\right)^{2}(-x_{d})+2\left(\frac{1}{6}\right)^{2}x_{q}-\left(\frac{2}{3}\right)^{2}x_{u}\right)+\left(2\left(-\frac{1}{2}\right)^{2}x_{\ell}-(-1)^{2}x_{e_{R}}\right)&=0,\end{split}$$

$$\mathcal{L}={\mathcal{L}}_{SM}+{\mathcal{L}}_{gauge}+{\mathcal{L}}_f+{\mathcal{L}}_{scalar}+{\mathcal{L}}_{yuk}.$$

$${\mathcal{L}}_{gauge}=-\frac{1}{4}{B}_{\mu \nu }^{{}^{\prime }}{B}^{{}^{\prime }\mu \nu }+ϵB_{\mu \nu }{B}^{{}^{\prime }\mu \nu },$$

$$B_{\mu \nu }^{{}^{\prime }}={\partial }_{\mu }{\left(Z_{B-3L_{\tau }}^{{}^{\prime }}\right)}_{\nu }-{\partial }_{\nu }{\left(Z_{B-3L_{\tau }}^{{}^{\prime }}\right)}_{\mu }.$$

$$\begin{array}{rl}{\mathcal{L}}_f& =\overline{\psi }/D\psi +\overline{\chi }/D\chi +\sum _{j=1}^3\overline{N_{R_j}}/D{N}_{R_j}-M_{\psi }\overline{\psi }\psi -M_{\chi }\overline{\chi }\chi -\frac{1}{2}\sum _{i,j=1,2}{M}_{ij}\overline{{\left(N_{R_i}\right)}^c}{N}_{R_j}\\ & -\sum _{i=1,2}\frac{1}{2}{y}_{i3}\overline{{\left(N_{R_3}\right)}^c}{N}_{R_i}\mathcal{S}-y_{\chi }\overline{{\left(\chi \right)}^c}\chi \Phi +h.c.,\end{array}$$

$$D_{\mu }\equiv \left({\partial }_{\mu }-i{g}_2\frac{{\tau }^a}{2}{W}_{\mu }^a-i{g}_{1}Y{B}_{\mu }-i{g}_{B-3L_{\tau }}{Y}_{B-3L_{\tau }}{Z}_{(B-3L_{\tau })\mu }\right).$$

$${\mathcal{L}}_{scalar}={\left(D_{\mu }H\right)}^{†}\left(D^{\mu }H\right)+{\left(D_{\mu }\Phi \right)}^{†}\left(D^{\mu }\Phi \right)+{\left(D_{\mu }\mathcal{S}\right)}^{†}\left(D^{\mu }\mathcal{S}\right)-V(H,\Phi ,\mathcal{S}),$$

$$\begin{array}{rl}{D}_{\mu }H& =\left({\partial }_{\mu }-i{g}_2\frac{{\tau }^a}{2}{W}_{\mu }^a-i{g}_{1}Y{B}_{\mu }\right)H,\\ D_{\mu }\Phi & =\left({\partial }_{\mu }-i{g}_{B-3L_{\tau }}{Y}_{B-3L_{\tau }}{Z}_{(B-3L_{\tau })\mu }\right)\Phi ,\\ D_{\mu }\mathcal{S}& =\left({\partial }_{\mu }-i{g}_{B-3L_{\tau }}{Y}_{B-3L_{\tau }}{Z}_{(B-3L_{\tau })\mu }\right)\mathcal{S}.\end{array}$$

$$\begin{array}{rl}V\left(H,\Phi ,\mathcal{S}\right)& ={\mu }_H^2\left(H^{†}H\right)+{\lambda }_H{\left(H^{†}H\right)}^2+{\mu }_{\Phi }^2|\Phi {|}^2+{\lambda }_{\Phi }|\Phi {|}^4+{\mu }_{\mathcal{S}}^2|\mathcal{S}{|}^2\\ & +{\lambda }_{\mathcal{S}}|\mathcal{S}{|}^4+{\lambda }_1\left(H^{†}H\right)|\Phi {|}^2+{\lambda }_{\mathcal{S}\Phi }|\Phi {|}^2|\mathcal{S}{|}^2+{\lambda }_2\left(H^{†}H\right)|\mathcal{S}{|}^2\\ & +\mu \left(\mathcal{S}\Phi \Phi +H.c\right).\end{array}$$

$$\begin{gathered} H=\left(\begin{array}{c}{G}^{+}\\ \frac{h+v_d+i{z}_1}{\sqrt{2}}\end{array}\right) \\ ,\text{leavevmode}\Phi =\frac{1}{\sqrt{2}}\left(\varphi +v_{\Phi }+i{z}_2\right),\text{leavevmode} \\ \mathcal{S}=\frac{1}{\sqrt{2}}\left(s+v_{\mathcal{S}}+i{z}_3\right). \end{gathered}$$

$$\left(\begin{array}{c}{h}_1\\ h_2\\ h_3\end{array}\right)=U\left({\theta }_{12},{\theta }_{13},{\theta }_{23}\right)\left(\begin{array}{c}h\\ \varphi \\ s\end{array}\right),$$

$$\begin{array}{rl}{\mu }_H^2=\frac{1}{2}\left(-2{\lambda }_H{v}_d^2-{\lambda }_2{v}_{\mathcal{S}}^2-{\lambda }_1{v}_{\Phi }^2\right),& \\ {\mu }_{\Phi }^2=\frac{1}{2}\left(-{\lambda }_1{v}_d^2-{\lambda }_{\Phi \mathcal{S}}{v}_{\mathcal{S}}^2-3\lambda v_{\mathcal{S}}{v}_{\Phi }-2\lambda \varphi v_{\Phi }^2\right),& \\ {\mu }_{\mathcal{S}}^2=\frac{-{\lambda }_2{v}_d^2{v}_{\mathcal{S}}-2{\lambda }_{\mathcal{S}}{v}_{\mathcal{S}}^3-{\lambda }_{\Phi \mathcal{S}}{v}_{\mathcal{S}}{v}_{\Phi }^2-\lambda v_{\Phi }^3}{2v_{\mathcal{S}}}.& \end{array}$$

$$\displaystyle\begin{split}h&=h_{1}c_{12}c_{13}+h_{2}s_{12}c_{13}+h_{3}s_{13}&\\ \phi&=c_{12}(c_{23}(h_{2}-h_{1}t_{12})+s_{23}(h_{3}\frac{c_{13}}{c_{12}}-s_{13}(h_{1}+h_{2}t_{12})))\\ s&=c_{12}c_{23}(h_{1}t_{12}t_{23}-s_{13}(h_{1}+h_{2}t_{12})-h_{2}t_{23}+h_{3}c_{12}c_{13}+h_{3}s_{12}t_{12}c_{13}).\end{split},$$

$$\{m_{h_{2,3}},m_A,{\theta }_{23}\}.$$

$$\begin{array}{rl}& 2v_d^2{\lambda }_H=m_{h_1}^2{c}_{12}^2{c}_{13}^2+m_{h_2}^2{s}_{12}^2{c}_{13}^2+m_{h_3}^2{s}_{13}^2,\\ & 2v_{\varphi }^2{\lambda }_{\Phi }=m_{h_1}^2(c_{12}{s}_{13}{s}_{23}+s_{12}{c}_{23}{)}^2+m_{h_2}^2(c_{12}{c}_{23}-s_{12}{s}_{13}{s}_{23}{)}^2+m_{h_3}^2{c}_{13}^2{s}_{23}^2,\\ & 2v_S^2{\lambda }_S=m_{h_1}^2{s}_{12}^2{s}_{23}^2+c_{12}^2\left(m_{h_1}^2{s}_{13}^2{c}_{23}^2+m_{h_2}^2{s}_{23}^2\right)-s_{12}{c}_{12}{s}_{13}(1-2s_{23}^2)\left(m_{h_1}^2-m_{h_2}^2\right)\\ & \text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}+m_{h_2}^2{s}_{12}^2{s}_{13}^2{c}_{23}^2+m_{h_3}^2{c}_{13}^2{c}_{23}^2-m_A^2,\\ & v_d{v}_S{\lambda }_2=c_{13}\left(s_{12}{c}_{12}{s}_{23}\left(m_{h_1}^2-m_{h_2}^2\right)+s_{13}{c}_{23}\left(-m_{h_1}^2{c}_{12}^2-m_{h_2}^2{s}_{12}^2+m_{h_3}^2\right)\right),\\ & v_S{v}_{\varphi }^2{\lambda }_{\varphi S}=v_{\varphi }{c}_{12}^2{s}_{23}{c}_{23}\left(m_{h_1}^2{s}_{13}^2-m_{h_2}^2\right)+v_{\varphi }{s}_{12}{c}_{12}{s}_{13}(2c_{23}^2-1)\left(m_{h_1}^2-m_{h_2}^2\right)\\ & \text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}-m_{h_1}^2{v}_{\varphi }{s}_{12}^2{s}_{23}{c}_{23}+m_{h_2}^2{v}_{\varphi }{s}_{12}^2{s}_{13}^2{s}_{23}{c}_{23}+m_{h_3}^2{v}_{\varphi }{c}_{13}^2{s}_{23}{c}_{23}+2m_A^2{v}_S,\\ & \mu =-\frac{\sqrt{2}{m}_A^2{v}_S}{v_{\varphi }^2}.\end{array}$$

$$m_{Z_{B-3L_{\tau }}}^2=\frac{9}{4}{v}_{\Phi }^2{g}_{B-3L_{\tau }}^2\left(1+4{\tilde{v}}^2\right),$$

$$\begin{array}{rl}-{\mathcal{L}}_{yuk}& =\sum _{i=1}^2\sum _{\alpha =e,\mu }{y}_{\alpha i}\text{leavevmode}\overline{L_{\alpha }}\tilde{H}{N}_{R_i}+y_{{\tau }_3}\overline{L_3}\tilde{H}{N}_{R_3}+Y\overline{\psi }\tilde{H}\chi +h.c.,\end{array}$$

$${\mathcal{M}}_{VLF}=\left(\begin{array}{cc}{M}_{\chi }& m_D\\ m_D& M_{\psi }\end{array}\right),$$

$$\left(\begin{array}{c}{\psi }_1\\ {\psi }_2\end{array}\right)=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{c}\chi \\ {\psi }^0\end{array}\right),$$

$$\tan \left(2\theta \right)=\frac{2m_D}{M_{\chi }-M_{\psi }}.$$

$$M_{{\psi }_1}\simeq M_{\chi }+m_D\sin 2\theta$$

$$M_{{\psi }_2}\simeq M_{\psi }-m_D\sin 2\theta .$$

$$M_{{\psi }^{\pm }}=M_{\psi }=M_{{\psi }_1}{\sin }^2\theta +M_{{\psi }_2}{\cos }^2\theta .$$

$$Y=-\frac{(M_{{\psi }_2}-M_{{\psi }_1})\sin 2\theta }{\sqrt{2}{v}_d}=-\frac{\Delta M\sin 2\theta }{\sqrt{2}{v}_d}.$$

$${\mathcal{M}}_{p-Dirac}=\left(\begin{array}{cc}{m}_{\chi }& M_{\chi }\\ M_{\chi }& m_{\chi }\end{array}\right),$$

$$\begin{array}{rl}{\mathcal{L}}_{p-Dirac}=& \overline{\left(\chi \text{leavevmode}\text{leavevmode}{\chi }^c\right)}\left(\begin{array}{cc}{m}_{\chi }& M_{\chi }\\ M_{\chi }& m_{\chi }\end{array}\right)\left(\begin{array}{c}{\chi }^c\\ \chi \end{array}\right)\\ & \equiv \overline{\left({\chi }^a\text{leavevmode}\text{leavevmode}{\chi }^b\right)}\left(\begin{array}{cc}{M}_{\chi }-m_{\chi }& 0\\ 0& M_{\chi }+m_{\chi }\end{array}\right)\left(\begin{array}{c}{\chi }^a\\ {\chi }^b\end{array}\right)\equiv \overline{\chi }{M}_{\chi }\chi ,\end{array}$$

$$\{g_{B-3L_{\tau }},\text{leavevmode}\tilde{v},\text{leavevmode}{v}_{\Phi },\text{leavevmode}{\theta }_{23},\text{leavevmode}\sin \theta ,\text{leavevmode}{M}_{{\psi }_1},\text{leavevmode}\Delta M\},$$

$$\begin{array}{rl}{V}^{(4)}\left(H,\mathcal{S},\Phi \right)& ={\lambda }_H|H{|}^4+{\lambda }_{\Phi }|\Phi {|}^4+{\lambda }_{\mathcal{S}}|\mathcal{S}{|}^4+{\lambda }_1|H{|}^2|\Phi {|}^2+{\lambda }_2|H{|}^2|\mathcal{S}{|}^2+{\lambda }_{\mathcal{S}\Phi }|\Phi {|}^2|\mathcal{S}{|}^2.\end{array}$$

$${\lambda }_H\ge 0,\text{leavevmode}\text{leavevmode}\text{leavevmode}{\lambda }_{\Phi }\ge 0,\text{leavevmode}\text{leavevmode}\text{leavevmode}{\lambda }_{\mathcal{S}}\ge 0$$

$${\lambda }_1+2\sqrt{{\lambda }_H{\lambda }_{\Phi }}\text{leavevmode}\ge 0,\text{leavevmode}{\lambda }_2+2\sqrt{{\lambda }_H{\lambda }_{\mathcal{S}}}\text{leavevmode}\ge 0,\text{leavevmode}{\lambda }_{\mathcal{S}\Phi }+2\sqrt{{\lambda }_{\Phi }{\lambda }_{\mathcal{S}}}\text{leavevmode}\ge 0,$$

Full text

I Introduction

The fact that our present universe has a large proportion of its energy density contained in a mysterious, non-luminous and non-baryonic form of matter, popularly known as dark matter (DM), has been very well established by now. Amongst notable evidences suggesting this, the galaxy cluster observations made by Fritz Zwicky  Zwicky:1933gu in 1933, observations of galaxy rotation curves in 1970’s by Rubin and collaborators Rubin:1970zza , the observation of the bullet cluster by Chandra observatory Clowe:2006eq and the measurements of cosmic microwave background (CMB) by several cosmology experiments, the latest of which is the Planck experiment Aghanim:2018eyx . In terms of density parameter $\Omega_{\rm DM}$ and $h=\text{Hubble Parameter}/(100\;\text{km}\leavevmode\nobreak\ \text{s}^{-1}\text{Mpc}^{-1})$, the present DM abundance is conventionally reported as Aghanim:2018eyx : $\Omega_{\text{DM}}h^{2}=0.120\pm 0.001$ ghosh1992 at 68% CL. This corresponds to around $26\%$ of the present universe’s energy density, filled up by DM. While all such evidences are purely based on gravitational interactions of DM, there exists motivations to expect that DM could have some other forms of weak interactions as well. Interestingly, if DM interactions with the standard model (SM) particles are similar to those of electroweak interactions, and particle DM’s mass also remain around the electroweak scale, such a DM can be thermally produced in the early universe, followed by its freeze-out, leaving a thermal relic very close to the observed DM abundance. This remarkable coincidence is often referred to as the weakly interacting massive particle (WIMP) miracle Kolb:1990vq .

Similarly, the origin of light neutrino masses and mixing have also been a mystery in the last few decades. While experimental evidence confirms the neutrino mass squared differences to be several order of magnitudes below the electroweak scale, their large mixing, in sharp contrast with the well known quark sector, leads to another puzzle  Tanabashi:2018oca . Due to the absence of the right handed neutrino, the SM can not accommodate light neutrino masses due to absence of neutrino-Higgs coupling at the renormalisable level. One can however, introduce non-renormalisable Weinberg operator Weinberg:1979sa $(LLHH)/\Lambda,L\equiv$ lepton doublet, $\Lambda\equiv$ unknown cut-off scale, at dimension five level, to account for tiny neutrino masses. Dynamical realisation of this operator leads to beyond standard model (BSM) scenarios Minkowski:1977sc ; GellMann:1980vs ; Schechter:1980gr where introduction of heavy singlet neutrinos take part in generating light neutrino masses through type I seesaw mechanism.

Motivated by these two problems in the SM, we consider a popular DM scenario based on vector-like fermions along with an extended gauge symmetry. This extra gauge symmetry, on top of stabilizing the DM, also plays a role in generating light neutrino masses and provides additional incentive of enhanced detection aspects. The DM is an admixture of a vector like singlet fermion and neutral component of a vector like $SU(2)_{L}$ doublet fermion, popularly known as singlet-doublet fermion DM Mahbubani:2005pt ; DEramo:2007anh ; Enberg:2007rp ; Cohen:2011ec ; Cheung:2013dua ; Restrepo:2015ura ; Calibbi:2015nha ; Cynolter:2015sua ; Bhattacharya:2015qpa ; Bhattacharya:2017sml ; Bhattacharya:2018fus ; Bhattacharya:2018cgx ; DuttaBanik:2018emv . While such singlet-doublet fermion extension of the SM can also have some other motivations like, for example, electroweak baryogenesis, leading to the observed baryon asymmetry of the universe Egana-Ugrinovic:2017jib , we confine ourselves to the discussion of DM and its relevant phenomenology. Typically, vector like fermion (VLF) DM in such scenarios are stabilised by an additional $Z_{2}$ symmetry under which the DM is odd while all the SM particles are even. A purely singlet VLF DM does not have any renormalizable portal interaction with the SM to generate correct thermal relic abundance. The purely doublet VLF, on the other hand, by virtue of its electroweak gauge interactions, annihilates a lot to SM particles, thus reducing its thermal relic abundance unless its mass is beyond a TeV. A purely doublet VLF also faces stringent constraints from DM direct detection experiments like LUX Akerib:2016vxi , PandaX-II  Tan:2016zwf ; Cui:2017nnn and XENON1T  Aprile:2017iyp ; Aprile:2018dbl because of large DM-nucleon scattering mediated by electroweak gauge bosons. Even an admixture of a vector-like singlet and a doublet fermion faces tight constraints from direct detection experiments. In order to overcome that, a small Majorana mass term is introduced to split the vector-like mass eigenstates into two pseudo-Dirac ones. This splitting results in inelastic $Z$ boson coupling to DM which can be prevented kinematically. The admixture of a singlet and a doublet fermion, therefore, remains as an interesting scenario as it can circumvent both of these problems (namely, under/over-abundant relic and too large direct detection cross section of pure singlet/doublet DM), thus allowing the possibility of sub-TeV DM. We embed the singlet-doublet fermion DM within a gauge symmetry based on the $U(1)_{B-3L_{\tau}}$ gauge charge. While gauged $B-L$ symmetric extension of the SM  Mohapatra:1980qe ; Marshak:1979fm ; Masiero:1982fi ; Mohapatra:1982xz has been one of the most widely studied and well motivated BSM frameworks, we restrict it to the third lepton family only in order to evade strong bounds on the $B-L$ gauge boson from the Large Hadron Collider (LHC) as well as flavour physics  Khachatryan:2016qkc ; Aaboud:2017sjh ; Aaboud:2017buh ; Sirunyan:2018exx ; Chun:2018ibr ; Ma:1997nq ; Ma:1998dp ; Okada:2012sp ; Lee:2010hf ; Pal:2003ip ; Chang:2000xy ; Ma:1998dr . In addition, such family non-universal neutral gauge boson is also motivating from flavour anomalies point of view. Several proposals have appeared in the literature in order to accommodate the reported anomalies in $B$-meson systems by incorporating additional flavoured gauge bosons (see, for example Chang:2009tx ; Crivellin:2015lwa ; Altmannshofer:2015mqa ; Allanach:2015gkd ; Ko:2017lzd ; Ko:2017yrd ; Allanach:2018lvl ). The possibility of a TeV scale neutral gauge boson enhances the production cross section of different components of the vector-like fermions in comparison to the usual singlet-doublet fermion DM models. Since such neutral gauge bosons also mediate DM annihilations, henc they also affect DM parameter space compared to the scenarios with universal $B-L$ gauge bosons 111A recent summary of such $Z^{\prime}$ mediated DM scenarios can be found in Blanco:2019hah .. Also, the requirement of anomaly cancellation in this model introduces a heavy singlet right handed neutrino (RHN) which, along with one or two more right handed neutrinos having vanishing $B-3L_{\tau}$ charges, can take part in generating light neutrino masses through usual type I seesaw mechanism. Such right handed neutrino, charged under the additional gauge symmetry, can be produced resonantly in colliders and can lead to exotic signatures like displaced vertex, if sufficiently long lived Das:2019fee . Apart from offering a potential DM candidate, generating light neutrino mass, and providing tantalizing collider signatures, the model also attempts to give a solution to the metastable nature of electroweak vacuum Isidori:2001bm ; Ellis:2009tp ; EliasMiro:2011aa ; Alekhin:2012py ; Buttazzo:2013uya ; Anchordoqui:2012fq ; Tang:2013bz ; Bhattacharya:2019fgs . The negative fermionic contribution (primarily due to top quark and VLFs) to the renormalisation group (RG) running of the Higgs quartic coupling is compensated by respective contributions from additional scalars in the model 222With VLF alone, it is possible to make the electroweak vacuum stable, if they have coloured charges, as shown in Gopalakrishna:2018uxn .. In particular, the constraints from the requirement of vacuum stability restricts the gauge coupling of TeV scale $B-3L_{\tau}$ gauge symmetry to g_{B-3L_{\tau}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.25 and SM Higgs coupling with VLF to Y\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.3. For such couplings the model also remains perturbative all the way upto the Planck scale. A recent attempt to constrain a similar Abelian gauge model with inverse seesaw origin of light neutrino mass from DM and vacuum stability criteria has appeared in Das:2019pua .

This paper is organised as follows: In section II we have introduced the new particles in our model and their interaction Lagrangian. In section III we have discussed the constraints on the model parameters arising from stability, unitarity and perturbativity of the scalar potential, electroweak precision observables, LHC searches and generation of light neutrino mass requirements. The details of the parameter space scan for the DM phenomenology is elaborated in section IV where in subsection IV.1 we have illustrated the relic density allowed parameter space and in subsection IV.2 direct search is discussed. The high scale validity and perturbativity of the model is elaborated in section V, where we have also chosen some of the benchmark points satisfying all relevant constraints in order to perform the collider analysis. The collider signatures of the model, along with discovery potential in the LHC is thoroughly explained in section VI. Finally in section VII we have concluded and summarised our findings.

II The Model: fields and interactions

We first tabulate all the particles of the model in table 1 with their corresponding gauge charges. On top of the familiar SM particles we have a few additional particles. In the fermion sector, there are three right handed neutrinos (RHN): $N_{R}$ of which $N_{R_{1,2}}$ have zero $B-3L_{\tau}$ charges and $N_{R_{3}}$ has non-zero charge under $U(1)_{B-3L_{\tau}}$. Apart from cancelling the triangle anomalies arising due to the $U(1)_{B-3L_{\tau}}$ symmetry (shown below) these right handed neutrinos help us to generate light neutrino masses via type I seesaw mechanism as we shall discuss in detail. We also have one VLF doublet $\psi^{T}:\left(\psi^{0}\leavevmode\nobreak\ \psi^{-}\right)$ and one VLF singlet $\chi$. The lightest physical state that arises from the mixing of these two will serve as the DM, while the charged components can be produced at the collider to give interesting signatures. In the scalar sector, apart from the SM Higgs doublet $H$ we have two more singlet scalars: $\mathcal{S}$ and $\Phi$ both charged under $U(1)_{B-3L_{\tau}}$. The additional SM singlet scalar fields take part in spontaneous symmetry breaking (SSB) of $U(1)_{B-3L_{\tau}}$, thus giving mass to the new heavy charge neutral gauge boson. While one such scalar is sufficient for spontaneous gauge symmetry breaking, we need both of them in the set up for phenomenological reasons. One of the scalars helps in providing the required texture of RH neutrino mass matrix through its non zero vev after spontaneous breaking of $U(1)_{B-3L_{\tau}}$. The other scalar splits the Dirac VLFs into pseudo-Dirac states, required to avoid stringent direct detection bounds, as we shall explicitly show while discussing the DM phenomenology.

The condition for anomaly cancellation comes from the triangle diagrams with gauge gauge bosons at the vertices:

[TABLE]

with “Tr” standing for trace and $T^{a,b,c}$ denoting the corresponding generators. This vanishes exactly as the VLFs contribute identically to the left and right-handed representation. However, the contributions from other fermions do not vanish in general due to their chiral nature. Demanding cancellation of gauge and gravitational anomalies we end up with the following non-trivial relations:

[TABLE]

where $x_{i}$ stands for the $U(1)_{B-3L_{\tau}}$ charges for all the particles appearing in table 1. The other anomalies, namely $SU(2)^{2}U(1)_{B-3L_{\tau}},SU(3)_{c}^{2}U(1)_{B-3L_{\tau}}$ are trivially cancelled. The $U(1)_{B-3L_{\tau}}$ charges of the relevant fermions are assigned in table 1 from which it is easy to check that all the anomalies are cancelled by taking all three fermion generations into account. This is a typical feature of non-universal gauge symmetry where anomalies are not cancelled generation wise, but they vanish only for all three generations combined. The $U(1)_{B-3L_{\tau}}$ charges of the scalar fields and vector fermions will be dictated through the interaction terms in the Lagrangian as we explain in the following paragraphs. It is worth mentioning that the minimal $B-L$ gauge symmetric model Mohapatra:1980qe ; Marshak:1979fm ; Masiero:1982fi ; Mohapatra:1982xz is anomaly free if three right handed neutrinos having $B-L$ charge $-1$ each are taken into account. However, this is not the only solution to anomaly cancellation conditions. There exist exotic charges of additional chiral fermions that can give rise to vanishing triangle anomalies Montero:2007cd ; Wang:2015saa ; Patra:2016ofq ; Nanda:2017bmi ; Bernal:2018aon . Similarly, the anomaly cancellation solution mentioned here for our model is not the only possible one, there exists non-minimal solutions for the same which we do not discuss in our work.

With this particle content at our disposal, now we can proceed to write the Lagrangian for this model. The Lagrangian contains four more parts on top of the SM Lagrangian:

[TABLE]

The gauge part of the Lagrangian is written as:

[TABLE]

with

[TABLE]

The first term in Eq. (4) is the kinetic term for the new $U(1)_{B-3L_{\tau}}$ gauge boson and the second term arises due to kinetic mixing between $U(1)_{Y}$ and $U(1)_{B-3L_{\tau}}$ gauge bosons. For simplicity we choose $\epsilon=0$ at the scale of $B-3L_{\tau}$ symmetry breaking as in Okada:2018ktp . The Lagrangian for the new fermion sector (excluding their couplings to the SM Higgs) reads:

[TABLE]

where under $SM\otimes U(1)_{B-3L_{\tau}}$, the covariant derivative is defined as:

[TABLE]

where the second and third terms on the right hand side will be present only for the VLF doublet while for VLF singlet and RHNs, only the first and the last terms are present. In the last term on right hand side of Eq. (7), $Y_{B-3L_{\tau}}$ corresponds to $B-3L_{\tau}$ charge of the respective fermion. The first two terms in Eq. (6) are the kinetic terms for the VLF doublet and VLF singlet respectively. The third term corresponds to the kinetic term of the RHNs. The mass term for the VLFs is given by the fourth and fifth term. The Majorana mass terms for the RHNs is given by the sixth and seventh term. The seventh term decides the $B-3L_{\tau}$ charge of the scalar singlet $\mathcal{S}$ (+3) as $N_{3}$ must have $B-3L_{\tau}$ charge of -3 for the sake of anomaly cancellation. The last term is the Majorana term for the VLF singlet, which is responsible for the pseudo-Dirac splitting of the VLFs. For the scalar sector the Lagrangian can be expressed as:

[TABLE]

where

[TABLE]

With one SM Higgs doublet and two non-standard singlets the renormalisable scalar potential takes the form:

[TABLE]

The last term is important as it provides mass for the pseudoscalar by explicitly breaking the global symmetry of the potential, in absence of which we would have ended up with one massless Goldstone boson. Note that, this term also determines the $B-3L_{\tau}$ charge for the second scalar singlet $\Phi$ (-3/2) as the charge for $\mathcal{S}$ is already determined from RHN Majorana mass term in Eq. 6 as discussed earlier. At the same time, this term also restricts the choice for the $B-3L_{\tau}$ charge of the VLFs (+3/4), as evident from the last term in Eq. (6). Thus, the charge assignment of the new scalars and the VLFs is completely determined by the desired interaction Lagrangian as well as mass spectrum.

Now, the SM Higgs field and the $B-3L_{\tau}$ sector scalars are expanded around their respective vacuum expectation value (VEV):

[TABLE]

The gauged $U(1)_{B-3L_{\tau}}$ is spontaneously broken as the two scalar singlets acquire non-zero VEV. Then the weak eigenstates of the scalars mix with each other. Thus, in order to obtain the physical mass eigenstates, we diagonalise the mass matrix as:

[TABLE]

where $U\left(\theta_{12},\theta_{13},\theta_{23}\right)$ is usual unitary matrix involving the mixing angles $\theta_{12},\theta_{13},\theta_{23}$ and complex phase $\delta=0$. We assume $h_{1}$ to be SM-like Higgs with a mass of 125 GeV. The minimisation conditions are given by:

[TABLE]

The weak eigenstates $(h,\leavevmode\nobreak\ \phi,\leavevmode\nobreak\ s)$ in terms of the physical eigenstates $(h_{1},\leavevmode\nobreak\ h_{2},\leavevmode\nobreak\ h_{3})$ and the mixing angles are given by:

[TABLE]

where $c_{ij}=\cos\theta_{ij}$, $s_{ij}=\sin\theta_{ij}$ and $t_{ij}=\tan\theta_{ij}$ with $\{i,j\}=1,2,3$ and $i\neq j$. It is easy to understand from here, for $\theta_{12},\theta_{13}=0,\theta_{23}\neq 0$ we revive the SM Higgs purely from $h_{1}$. Therefore, $\theta_{12}$ and $\theta_{13}$ are constrained from Higgs data, while $\theta_{23}$ is a free parameter. After diagonalising, we are left with three CP-even scalars denoted by $h_{1,2,3}$. We also have three CP-odd pseudoscalars which, after diagonalising to their physical mass eigenstates, are referred to as $A$ and $G_{1,2}$, out of which $G_{1,2}$ turn out to be the Goldstone modes of $B-3L_{\tau},Z$ gauge bosons giving $m_{G_{1}}=m_{G_{2}}=0$. Now, $h_{1}$ is the lightest CP-even Higgs that has been seen at the LHC, hence $m_{h_{1}}=125\leavevmode\nobreak\ \rm GeV$. Also, the mixing angles $\theta_{12}$ and $\theta_{13}$ are constrained from Higgs data which we shall elaborate in section III. Essentially the scalar sector has the following free parameters:

[TABLE]

All the quartic couplings appearing in the scalar potential can be expressed in terms of the physical masses and mixings as follows:

[TABLE]

After SSB we are also left with a new massive charge neutral gauge bosons corresponding to broken $U(1)_{B-3L_{\tau}}$. The mass of the new gauge boson is given by:

[TABLE]

where $\tilde{v}=\frac{v_{\mathcal{S}}}{v_{\Phi}}$. Eq. (14) is important in our analysis as depending on m_{Z_{B-3L_{\tau}}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}2.4\leavevmode\nobreak\ \rm GeV 333We choose this conservative lower bound in order to be in agreement with relevant experimental constraints on flavoured gauge bosons Khachatryan:2016qkc ; Aaboud:2017sjh ; Aaboud:2017buh ; Sirunyan:2018exx ; Chun:2018ibr . we can constrain our parameter space. Finally, the Lagrangian for the Yukawa sector involving the SM Higgs reads as:

[TABLE]

where the first two terms are the interactions of the SM leptons with the RHNs. Note that, $N_{R_{3}}$ can only have interaction with the third generation SM leptons because of the $U(1)_{B-3L_{\tau}}$ charge assignment. The last term is the mixing of the two VLFs mediated by SM Higgs. This gives rise to the two physical eigenstates for the VLFs as mentioned below.

II.1 VLF mass eigenstates

The Dirac mass matrix in the basis $\{\chi,\psi^{0}\}$, containing the singlet and doublet VLF can be written as:

[TABLE]

where m_{D}=\frac{Yv_{d}}{\sqrt{2}}<<M_{\chi}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}M_{\psi}. Thus, the physical eigenstates arise as:

[TABLE]

where $\theta$ is the mixing angle given by:

[TABLE]

The mass for the physical eigenstates are:

[TABLE]

The lightest physical charge neutral fermion from above is a viable DM candidate in this model and we choose it to be $\psi_{1}$ $(M_{\psi_{1}}<M_{\psi_{2}})$. The DM is naturally stable due to our particular choice of the $U(1)_{B-3L_{\tau}}$ charge. In the small mixing limit the charged component of the VLF doublet $\psi^{\pm}$ acquires a mass:

[TABLE]

For small $\sin\theta$ ($\approx 0$) limit, $M_{\psi^{\pm}}\simeq M_{\psi_{2}}$. From Eq. (18), we see that the VLF Yukawa $Y$ is related to the mass difference between the two physical eigenstates, and is no more a free parameter:

[TABLE]

II.2 Pseudo-Dirac mass splitting

Due to presence of the Majorana term : $y_{\chi}\overline{\chi^{c}}\chi\Phi$, the pseudo-Dirac mass matrix for the VLF singlet $\chi$ can be expressed in the basis $(\chi^{c}\leavevmode\nobreak\ \leavevmode\nobreak\ \chi)^{T}$:

[TABLE]

where $m_{\chi}=y_{\chi}\frac{v_{\Phi}}{\sqrt{2}}\leavevmode\nobreak\ \ll M_{\chi}$. The mass matrix can be expressed in terms of physical states as:

[TABLE]

where $(\chi^{a}\leavevmode\nobreak\ \leavevmode\nobreak\ \chi^{b})^{T}$ is the physical pseudo-Dirac eigenstate. Since $m_{\chi}<<M_{\chi}$, which we can always assume the last equality in above equation by taking either $y_{\chi}$ or $v_{\Phi}$ or both to be small. In the presence of the Majorana term, the singlet $\chi$ is split into two pseudo-Dirac states, $\{\chi^{a},\leavevmode\nobreak\ \chi^{b}\}$ which are propagated into the physical states $\{\psi_{1}^{a,b},\leavevmode\nobreak\ \psi_{2}^{a,b}\}$ via VLF mixing. As these two pseudo-Dirac states ($\psi^{a,b}$) are nearly degenerate i.e, $\delta m\sim\mathcal{O}\left(100\rm keV\right)$, we consider them to be a single state $\psi=\left(\psi^{a}\leavevmode\nobreak\ \leavevmode\nobreak\ \psi^{b}\right)^{T}$ with mass $M_{\psi}$, identical to the Dirac mass of $\psi$. This will not make difference in DM relic abundance calculations, while for direct detection, such splitting will play a crucial role in preventing spin independent elastic scattering mediated by neutral gauge bosons.

III Constraints on the model parameters

The phenomenology of the model is mainly dictated by following free parameters :

[TABLE]

where $\Delta M=M_{\psi_{2}}-M_{\psi_{1}}$ is the difference between heavy and light VLF mass eigenstates. All these parameters are important for both DM as well as collider phenomenology as we shall see later. But before going into the details of parameter space scan, here we would like to explain how different parameters arising in the model are already constrained from theoretical as well as existing experimental bounds. Especially existing collider bound on the mass of the neutral gauge boson $Z_{B-3L_{\tau}}$ puts stringent constraint on the model parameters. Apart from that, there are bounds from stability of the scalar potential and perturbativity of dimensionless couplings, collider bounds on non-standard scalar masses and mixings, and bounds from light neutrino mass.

III.1 Stability, perturbativity and tree-level unitarity

III.1.1 Stability

Stability of the scalar potential is mainly dictated by the quartic terms of the scalar potential, $V(H,\Phi,\mathcal{S})$ which is defined as:

[TABLE]

In order to ensure the bounded-from-below condition in any field direction, the quartic couplings of the potential (Eq.(25)) must obey the following co-positivity conditions Kannike:2012pe ; Chakrabortty:2013mha :

[TABLE]

III.1.2 Perturbativity

To prevent perturbative breakdown of the model, all quartic, Yukawa and gauge couplings should obey the following limits at any energy scale:

[TABLE]

where $j=1,2$ and $\alpha=e,\mu$.

III.1.3 Tree-level unitarity

The quartic couplings of the scalar potential which are shown in Eq. (25) are also constrained from the following tree level perturbative unitarity conditions Horejsi:2005da ; Bhattacharyya:2015nca ; Kang:2013zba :

[TABLE]

where, $x_{1,2,3}$ are the cubic roots of the polynomial equation detailed in Appendix C.

III.2 Constraint from electroweak precision observables (EWPO)

Since our model has two BSM scalars and two vector like fermions, hence there should be corrections to the SM electroweak precision observables (EWPO) i.e, $S$,$T$,$U$ parameters Peskin:1991sw ; delAguila:2008pw ; Erler:2010sk ; Cynolter:2008ea . Here we would like to estimate the effect of the BSM particles on those parameters. We have four parameters, namely $\hat{S}$, $\hat{T}$, $W$ and $Y$ Barbieri:2004qk where $\hat{S}$ is related to the Peskin-Takeuchi parameter $S$: $\hat{S}=\frac{\alpha S}{4s_{w}^{2}}$ and $\hat{T}=\alpha T$ where $\alpha$ is the fine structure constant and $s_{w}\equiv\sin\theta_{w}$ corresponds to sine of the Weinberg angle $\theta_{w}$. The parameters, $W$ and $Y$ on the other hand, are new set of parameters. $T$ parameter is more significant for small mixing in the scalar sector and the constraint on $T$-parameter is parametrised by the following data Tanabashi:2018oca : $\Delta T=T^{\rm xSM}-T^{\rm SM}=0.07\pm 0.12$, where we consider contributions from both the VLFs and the non-standard scalars to $T^{\rm xSM}$. In this situation $\Delta T$ is given by Barger:2007im ; Ghosh:2015apa :

[TABLE]

which indicates how much the oblique parameter is shifted from the SM value. Now,

[TABLE]

and $T^{\rm SM\leavevmode\nobreak\ Higgs}(h_{1})$ can be obtained by using the decoupling limits $s_{12}\to 0$ and $s_{13}\to 0$ in Eq.(30). The contribution from VLF DM is followed as,

[TABLE]

where

[TABLE]

In Eq. (30) $m_{h_{1}}$ refers to the SM Higgs boson with mass 125 GeV, while $m_{h_{2,3}}$ are the two non-standard Higgs bosons appearing in our model. $m_{W}$ and $m_{Z}$ are the masses of SM $W$ and $Z$ bosons respectively. The expression for $T$-parameter corresponding to the contribution from the VLFs is given by Eq. (31), where $g_{2}$ is the $SU(2)_{L}$ SM gauge coupling. The $\Pi$’s appearing in the expression are given as in Eq. (32) and these correspond to the gauge boson propagator correction due to the VLFs. Here ‘div’ is the usual expression that appears in dimensional regularisation: ${\rm div}=\frac{1}{\epsilon}+\ln 4\pi-\gamma_{\epsilon}$, with $\gamma_{\epsilon}=0.577$ is the Euler-Mascheroni constant ($\epsilon=4-d,d\equiv$ spacetime dimension in dimensional regularisation). Note that the divergence appearing in the last term in Eq. (31) (due to the divergences in Eq. (32)) is cancelled by the first two terms . The physical mass eigenstates appearing in this case are $M_{\psi}$, $M_{\psi_{1}}$ and $M_{\psi_{2}}$ ($M_{\psi}\equiv M_{\psi^{\pm}}$ according to Eq. (21)). Once more we would like to remind that $M_{\psi_{2}}\approx M_{\psi^{\pm}}$ under small mixing limit as apparent from Eq. (21).

The bound on $\hat{S}$ comes from a global fit: $10^{3}\hat{S}=0.0\pm 1.3$ Barbieri:2004qk . For $S$-parameter, we consider contribution only due to the VLFs444As only $T$-paramter is important for scalar extension of the SM in small mixing case Ghosh:2015apa ; Barger:2007im . as given by Bhattacharya:2018fus ; Cynolter:2008ea :

[TABLE]

where $g_{2}$ is the $SU(2)_{L}$ gauge coupling. The expression for vacuum polarization for identical masses (at $q^{2}=0$) Cynolter:2008ea :

[TABLE]

For two different masses ($m_{i}\neq m_{j}$) the expression for vacuum polarization reads Cynolter:2008ea :

[TABLE]

Note that all the divergences appearing in Eq. (34) and (35) along with the renormalization scale $\mu_{EW}$, are cancelled on substitution in Eq. (33).

We have constrained the two most important free parameters of our model, namely DM mass $M_{\psi_{1}}$ and $\Delta M$ using these constraints. These are depicted in Fig. 1. On the left hand side (LHS) of Fig. 1 we have shown the allowed values of $M_{\psi_{1}}$ and $\Delta M$ that obey the constraint from $T$-parameter given by Eq. (29). What we see from the plot in the LHS of Fig. 1, for small VLF mixing ($\sin\theta\leq 0.5$) it is always possible to get large \Delta M\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}500\leavevmode\nobreak\ \rm GeV for DM mass upto 1 TeV within the permissible range of the $T$-parameter. For \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}0.7, however, the parameter space a is constrained as $\Delta M$ as large as 1 TeV can not be achieved for low DM mass. In the RHS of Fig. 1 we have shown the allowed range of $S$-parameter in $M_{\psi_{1}}$-$\Delta M$ plane. Here we see, for small $\sin\theta\sim 0.1$ the whole parameter space is allowed (gray region), however for larger \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}0.3 one has to stick to lower DM mass. Thus, $S$ parameter constraints the DM mass for large $\sin\theta$. But in any case it is always possible to have a large $\Delta M$ within the observed range of $S$ and $T$ parameter.

III.3 Constraint on $Z_{B-3L_{\tau}}$ mass from LHC

Experimental limits from LEP II constrains such new gauge sector by putting a lower bound on the ratio of new gauge boson mass to the new gauge coupling $M_{Z^{\prime}}/g^{\prime}\geq 7$ TeV Carena:2004xs ; Cacciapaglia:2006pk . The corresponding bounds from the LHC experiment have become stronger than this by now. As the main motivation to choose family non-universal gauge boson is to have weaker collider bounds on its mass, hence it is of utmost importance to realise what choice of the free parameters can give rise to right $Z_{B-3L_{\tau}}$ mass satisfying LHC bounds.

Search for heavy neutral Higgs and $Z_{B-3L_{\tau}}$ resonances have been performed at the LHC Aaboud:2017sjh , with the assumption that the heavy resonances decay to $\tau^{+}\tau^{-}$ final states. These searches rule out $m_{Z_{B-3L_{\tau}}}<2.42\leavevmode\nobreak\ \rm TeV$ at 95% CL for sequential SM and $m_{Z_{B-3L_{\tau}}}<2.25\leavevmode\nobreak\ \rm TeV$ at 95% CL for non-universal $G\left(221\right)$ model. We choose m_{Z_{B-3L_{\tau}}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}2.5\leavevmode\nobreak\ \rm TeV for a conservative limit. This puts a bound on three parameters in our model, namely: $\{g_{B-3L_{\tau}},\tilde{v},v_{\Phi}\}$. This is shown in Fig. 2, where each contour corresponds to $m_{Z_{B-3L_{\tau}}}=2.5\leavevmode\nobreak\ \rm TeV$ and hence the region right to each of the contours is allowed from collider constraint. We have chosen three different VEVs $v_{\Phi}:\{1.0,2.0,3.0\}\leavevmode\nobreak\ \rm TeV$ corresponding to red, green and blue contours respectively. As it is seen, larger $v_{\Phi}$ allows larger gauge boson mass, which is in accordance with Eq. (14). One should note here, Eq. (14) has a consequence. It does not allow to fix either the gauge coupling or $\tilde{v}$ for a fixed $v_{\Phi}$. As a result, for the parameter space scan we have varied both $g_{B-3L_{\tau}}$ and $\tilde{v}$ for a fixed $v_{\Phi}$ to keep the $Z_{B-3L_{\tau}}$ mass in the right ballpark (\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}2.5\leavevmode\nobreak\ \rm TeV). Here we would like to mention that a combination of $B_{s}$-$\overline{B_{s}}$ mixing from flavour physics, together with ATLAS’ $Z^{{}^{\prime}}$ search puts a bound on third family hypercharge models by requiring $m_{Z^{{}^{\prime}}}>1.9\leavevmode\nobreak\ \rm TeV$ Allanach:2019mfl . However, as we are considering even more conservative bound, our models is safe from such constraints arising from flavour physics measurements.

III.4 Bounds on singlet scalar from collider

The bounds on singlet scalars typically arise due to their mixing with the SM Higgs boson. The bound on such scalar mixing angles would come from both theoretical and experimental constraints Robens:2015gla ; Chalons:2016jeu . In case of scalar singlet extension of SM, the strongest bound on scalar-SM Higgs mixing angle ($\theta_{m}$) comes form $W$ boson mass correction Lopez-Val:2014jva at NLO for $250{\rm\leavevmode\nobreak\ GeV}\lesssim m_{h_{2}}\lesssim 850$ GeV as (0.2\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}\sin\theta_{m}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.3) where $m_{h_{2}}$ is the mass of other physical Higgs. Whereas, for $m_{h_{2}}>850$ GeV, the bounds from the requirement of perturbativity and unitarity of the theory turn dominant which gives $\sin\theta_{m}\lesssim 0.2$. For lower values i.e. $m_{h_{2}}<250$ GeV, the LHC and LEP direct search Khachatryan:2015cwa ; Strassler:2006ri and measured Higgs signal strength Strassler:2006ri restrict the mixing angle $\sin\theta_{m}$ dominantly ($\lesssim 0.25$). The bounds from the measured value of EW precision parameter are mild for $m_{h_{2}}<1$ TeV. In our analysis we have two singlet scalars which we intend to keep below TeV range. Now considering all the possible bounds, we make conservative choices of the mixing angles (with SM Higgs) as $\sin\theta_{12},\sin\theta_{13}\sim 0.1$. We also fix v_{\Phi}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}1\leavevmode\nobreak\ \rm TeV which helps in keeping the perturbativity of the theory intact. The other mixing angle $\sin\theta_{23}$ is a free parameter. We keep it below 0.2.

III.5 Neutrino mass

The light neutrino mass matrix can be generated via type I seesaw mechanism

[TABLE]

where

$M_{D}=\begin{pmatrix}y_{e_{1}}v_{d}&y_{e_{2}}v_{d}&0\\ y_{\mu_{1}}v_{d}&y_{\mu_{2}}v_{d}&0\\ 0&0&y_{\tau_{3}}v_{d}\end{pmatrix}$ and $M_{R}=\begin{pmatrix}M_{11}&M_{12}&y_{13}v_{\mathcal{S}}\\ M_{12}&M_{22}&y_{23}v_{\mathcal{S}}\\ y_{13}v_{\mathcal{S}}&y_{23}v_{\mathcal{S}}&0\end{pmatrix}$

as obtained from the Yukawa interaction of SM leptons in Eq. (15) and singlet neutral fermions in Eq. (6). Now, the neutrino mixing angles and mass squared differences are precisely measured from neutrino oscillation experiments PhysRevD.98.030001 . This, in turn, puts bound on model parameters including the VEV $v_{\mathcal{S}}$ and relevent Yukawa couplings. This can be understood from the light neutrino mass matrix itself. Diagonalising the mass matrix in Eq. (36) with the usual $3\times 3$ PMNS matrix (choosing the charged lepton mass matrix diagonal) gives the light neutrino masses. We can choose: $y_{e_{1}}=y_{e_{2}}=y_{\mu_{1}}=y_{\mu_{1}}\equiv y_{l}$, $y_{13}=y_{23}\equiv y$, $M_{11}=M_{22}=M$. With this, if we assume $v_{\mathcal{S}}\sim\mathcal{O}\left(\rm TeV\right)$ and $y_{l}\simeq y_{\tau_{3}}\sim\mathcal{O}(10^{-7})$, then we can produce correct order of light neutrino mass for $y\sim 0.1$ and $M\simeq 1\leavevmode\nobreak\ \rm TeV$. Even if we take $M\simeq 10\leavevmode\nobreak\ \rm TeV$ or a different order of magnitude for $v_{S}$, correct light neutrino mass can still be obtained with Yukawa couplings of the similar order. However, in that case, the RHNs are beyond the present collider reach. We have kept $v_{\mathcal{S}}\sim\mathcal{O}(\rm TeV)$ such that $Z_{B-3L_{\tau}}$ can be produced at the coillders, which determines the collider signature for this model. And as we have shown above, this choice is not in contradiction with the neutrino mass generation. For simplicity, in our analysis we shall assume that the annihilation of the DM to RHN final state is kinematically forbidden. In that case the Yukawa couplings $y$, $y_{l}$ and $y_{\tau_{3}}$ do not play important role in the DM or collider analysis of this model, hence we can fix them to produce the neutrino mass (as well as mixing) in the right ballpark without disturbing the outcome of the DM or collider phenomenology. Note that, the requirement of generating correct neutrino mass does not put a very tight constraint on the choice of the VEV $v_{\mathcal{S}}$. Thus, in this model, the DM sector and neutrino sector are closely connected even though the bounds on dark sector from right neutrino mass requirement is not very stringent. It should be noted from the structure of $M_{R}$ that if we had considered a singlet scalar having $B-3L_{\tau}$ charge $6$ instead of $3$, $M_{R}$ will have $3-3$ element non-zero but $1-3,2-3$ elements zero. This, as can be checked by using the light neutrino mass formula in Eq. (36), will give rise to a phenomenologically unacceptable light neutrino mass matrix. This once again justifies the choice of singlet scalars and their $B-3L_{\tau}$ charges made in our model.

Here we would also like to mention that the TeV scale RHNs have decay lifetime $\tau_{N}\sim 10^{-13}\leavevmode\nobreak\ \rm sec$ considering SM neutrino with scalar final state: $N\to\nu,h_{1}$ for our choice of Yukawa couplings. This shows that the RHNs do not contribute to the DM relic abundance as $\tau_{N}<<\tau_{universe}(\sim 10^{17}\leavevmode\nobreak\ \rm sec)$. Also, since they decay very fast ($<<1\leavevmode\nobreak\ \rm sec$) to SM final states, they do not perturb the standard Big Bang Nucleosynthesis (BBN) picture and hence unconstrained from BBN data. However, for certain choices of lightest neutrino mass, RHN having such gauge interactions can be long lived enough to give interesting collider signatures like displaced vertices, as studied recently by the authors of Das:2019fee . Before ending this subsection, we note that, although neutrino mass and mixing do not constrain the mass of $B-3L_{\tau}$ gauge boson directly, constraints on neutrino non-standard interactions (NSI) can be used to set a lower bound on such gauge boson mass as $m_{Z_{B-3L_{\tau}}}/g_{B-3L_{\tau}}>4.8$ TeV Heeck:2018nzc . Since we are considering the stronger bounds from the LHC and LEP II in our analysis, such weaker bounds are trivially satisfied.

IV Dark matter phenomenology

The lightest charge neutral state, $\psi_{1}$ in VLF sector is the stable DM candidate in our model. It is naturally stable in this set-up precisely due to the $B-3L_{\tau}$ charge assignment (discussed in the subsection II.1). In this section we have explored in detail the parameter space appearing in Eq. (24) allowed by observed relic abundance ($\Omega h^{2}=0.120\pm 0.001$Aghanim:2018eyx ) and direct search limits (particularly from the XENON 1T experiment Aprile:2018dbl ).

IV.1 Relic abundance of the DM

Relic density of DM, $\psi_{1}$ is governed by SM Higgs ($h_{1}$) and heavy Higgs ($h_{2},\leavevmode\nobreak\ h_{3}$) mediated annihilation and co-annihilation types number changing processes along with SM gauge boson ($Z,\leavevmode\nobreak\ W^{\pm},\leavevmode\nobreak\ \gamma$) and additional heavy gauge boson ($Z_{B-3L_{\tau}}$) mediated annihilation and co-annihilation type processes. All the relevant Feynman graphs contributing to the DM relic abundance are listed in Appendix D. The number density of DM can be computed by solving the Boltzmann equation Kolb:1990vq of the form:

[TABLE]

with $n=n_{\psi_{1}}$ and $H$ being the Hubble parameter in radiation dominated universe. All types of DM number changing processes are taken into account inside ${\langle\sigma v\rangle}_{eff}$ Griest:1990kh ; Edsjo:1997bg which is given by

[TABLE]

In above equation, $g_{eff}$ is defined as the effective degrees of freedom, given by

[TABLE]

where $\bar{g}_{1},\leavevmode\nobreak\ \bar{g}_{2}$ and $\bar{g}_{3}$ are the internal degrees of freedom of $\psi_{1},\leavevmode\nobreak\ \psi_{2}\rm\leavevmode\nobreak\ and\leavevmode\nobreak\ \psi^{\pm}$ respectively, and $x=x_{f}=\frac{M_{\psi_{1}}}{T_{f}}$, where $T_{f}$ is the freeze out temperature.

Relic density of $\psi_{1}$ one can approximately expressed as Bhattacharya:2016ysw :

[TABLE]

where $x_{f}\approx 20$ and $g_{*}=106.7$, the degrees of freedom for all SM particles . Note here that we have not used the above approximate formula for computing DM ($\psi_{1}$) relic density. In order to calculate relic density, we have used the package MicrOmegas Belanger:2001fz for which the model files are generated from LanHEP Semenov:2014rea .

To see the behaviour of DM ($\psi_{1}$) relic density, we have fixed the VEV $v_{\Phi}=3.0\leavevmode\nobreak\ \rm TeV$ such that the mass of $Z_{B-3L_{\tau}}$ is always above the collider bound ($m_{Z_{B-3L_{\tau}}}>2.5$ TeV) for suitable choices of $\tilde{v}$ and $g_{B-3L_{\tau}}$ as shown in Fig. 2. We have also fixed the masses of all the non-standard scalars as $\{m_{h_{2}},\leavevmode\nobreak\ m_{h_{3}},\leavevmode\nobreak\ m_{A}\}=\{200,\leavevmode\nobreak\ 300,\leavevmode\nobreak\ 250\}\leavevmode\nobreak\ \rm GeV$ obeying existing collider bounds as described in subsection III.4 . All scalar mixing angles are also kept fixed: $\{\sin\theta_{12},\leavevmode\nobreak\ \sin\theta_{13},\leavevmode\nobreak\ \sin\theta_{23}\}=\{0.1,0.1,0.2\}$. We have kept fixed the above parameters throughout our analysis. For the above choice of free parameters, other dependent quartic couplings are determined by Eq. (13). With this choice of parameters, we first illustrate how the relic abundance of the DM varies with DM mass for different choices of VLF mixing $\sin\theta$, $\Delta M$, the new gauge coupling $g_{B-3L_{\tau}}$ and the ratio of the VEVs $\tilde{v}$, keeping all other parameters fixed at their values mentioned before.

In the top left panel of Fig. 3 we have shown how the relic abundance of the DM changes with its mass for different choices of the VLF mixing $\sin\theta:\{0.05,0.1,0.5\}$ in solid black, black dashed and black dot-dashed lines respectively. We have kept $\Delta M=50\leavevmode\nobreak\ \rm GeV$ fixed and chose $g_{B-3L_{\tau}}=0.2$ with $\tilde{v}=1.3$ such that $Z_{B-3L_{\tau}}$ mass satisfies the collider bound ($>2.5$ TeV). The very first feature that one should note is the presence of three major resonant drops due to SM Higgs ($M_{\psi_{1}}\sim m_{h_{1}}/2$), SM $Z\leavevmode\nobreak\ (M_{\psi_{1}}\sim m_{Z}/2)$ and new gauge boson $Z_{B-3L_{\tau}}\leavevmode\nobreak\ (M_{\psi_{1}}\sim m_{Z_{B-3L_{\tau}}}/2)$. As one can notice, with the increase in $\sin\theta$, the DM becomes more and more under-abundant as the annihilation via $Z$ and Higgs bosons $h_{1,2,3}$ become more dominant, increasing the total annihilation cross-section. $\Delta M$ in this case is large enough and we can safely ignore the effects of co-annihilation. In the top right panel of Fig. 3 we have illustrated how relic abundance varies with the DM mass for three choices of $\Delta M:\{5,100,500\}\leavevmode\nobreak\ \rm GeV$ shown in solid black, dashed black and dot-dashed black curves respectively while keeping $\sin\theta$ fixed at 0.1 along with $g_{B-3L_{\tau}}=0.2$ and $\tilde{v}=1.3$. For small $\Delta M$ the co-annihilation plays dominant role, making the DM under-abundant. On the other hand, for large $\Delta M$, co-annihilation becomes sub-dominant, and as a result the Higgs ($h_{1,2,3}$) mediated resonance peaks become more prominent. Here one can notice that $h_{1,2,3}$ mediated resonances are more visible for large value of $\Delta M$, as the corresponding annihilation processes dominate. It is interesting to note from this plot that as we increase $\Delta M$ from 100 GeV to 500 GeV, the relic abundance decreases for $M_{\psi_{1}}\gtrsim 50$ GeV. This is due to the fact that, although increased $\Delta M$ decreases the efficiency of coannihilation processes, it increases the VLF coupling with SM Higgs, thereby increasing the scalar mediated annihilation processes.

Finally, in the bottom panel of Fig. 3 we have shown two different sets of $\{\tilde{v},g_{B-3L_{\tau}}\}:\{6,\leavevmode\nobreak\ 0.05\};\{1,0.3\}$ in solid black and dashed black curves respectively. This gives two different resonances at DM mass $M_{\psi_{1}}$ $\sim\frac{2.7}{2}\leavevmode\nobreak\ \rm TeV$ and $M_{\psi_{1}}\sim\frac{3.0}{2}\leavevmode\nobreak\ \rm TeV$ due to two different masses of $Z_{B-3L_{\tau}}$. The nature of the two curves is almost identical except for two different resonances at $m_{Z_{B-3L_{\tau}}}/2$. This clearly tells the fact that the dependence of DM relic abundance on the new gauge coupling $g_{B-3L_{\tau}}$ is mild compared to the dependence on the VLF mixing $\sin\theta$ and mass difference $\Delta M$. In each of the plots the dashed red straight line corresponds to the central value of Planck limits on DM relic abundance Aghanim:2018eyx .

We now scan all the free parameters of our analysis in the following range:

[TABLE]

Here we remind the readers once again that the choices for other parameters are kept fixed in the analysis as: $\sin\theta_{12}=0.1,\leavevmode\nobreak\ \sin\theta_{13}=0.1,\leavevmode\nobreak\ \sin\theta_{23}=0.2,\leavevmode\nobreak\ m_{h_{2}}=200\;{\rm GeV},\leavevmode\nobreak\ m_{h_{2}}=300\;{\rm GeV}\leavevmode\nobreak\ ,m_{A}=250\;{\rm GeV}{\rm and}\leavevmode\nobreak\ v_{\Phi}=3.0\;{\rm TeV}$ . The gauge coupling $g_{B-3L_{\tau}}$ is varied upto $0.3$, such that the model remains valid at high scale which we shall discuss in detail later. Throughout the scan we have ensured that m_{Z_{B-3L_{\tau}}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}2.5\leavevmode\nobreak\ \rm TeV by properly adjusting $g_{B-3L_{\tau}}$ and $\tilde{v}$, keeping $v_{\Phi}$ fixed at 3 TeV as mentioned earlier.

The allowed parameter space from relic density requirement set by Planck experiment is shown in Fig. 4 in $M_{\psi_{1}}$-$\Delta M$ plane. In the top left corner of Fig. 4 we have shown this parameter space for different ranges of the VLF mixing $\sin\theta$ shown in red $(0.01\leq\sin\theta<0.05)$, green $(0.05\leq\sin\theta<0.1)$, blue $(0.1\leq\sin\theta<0.3)$ and black $(0.3\leq\sin\theta<0.5)$ where $0.05\leq g_{B-3L_{\tau}}\leq 0.3$. The relic abundance criteria is satisfied by moderate to large $\sin\theta$, while small $\sin\theta$’s are confined near SM $Z$ and SM Higgs resonance and near $\Delta M\sim 10\leavevmode\nobreak\ \rm GeV$ for M_{\psi_{1}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}100\leavevmode\nobreak\ \rm GeV. In order to understand the pattern more clearly we have chosen a fixed $\sin\theta=0.2$ and plotted the same parameter space in the top right corner of Fig. 4. For DM mass around $M_{\psi_{1}}\sim 20\leavevmode\nobreak\ \rm GeV$ there are only a few annihilation channels open for the DM. Now, for small $\Delta M$ co-annihilation comes into picture, increasing the effective annihilation cross-section (in Eq. (38)). This causes the initial under abundance for small $\Delta M$. On further increasing $\Delta M$, co-annihilation becomes sub-dominant. As a result the DM becomes over abundant since the effective cross-section (Eq. (38)) diminishes. For $M_{\psi_{1}}\simeq 40-70\leavevmode\nobreak\ \rm GeV$ there is a huge under abundant region (green points) due to Z and SM Higgs resonances. Upon further increasing DM mass, we again get under-abundant regions in low $\Delta M$ region due to enhanced coannihilation and high $\Delta M$ region due to increased scalar portal annihilations as well as the resonance of the heavy scalars ($h_{2,3}$), which we noticed while discussing the behaviour of Fig. 3 as well. As mentioned earlier, the Higgs portal Yukawa $Y$ becomes large enough in such a case (being proportional to $\Delta M$ for a fixed $\sin\theta$) resulting a net increase in the annihilation cross-section. For DM mass $\sim 1\leavevmode\nobreak\ \rm TeV$ there is a huge overabundant region in the parameter space. This is due to the $1/M_{\psi_{1}}^{2}$ suppression in the annihilation cross-section due to heavy DM mass. Correct relic abundance is still possible to reach at a very large $\Delta M$ as then the Yukawa $Y$ becomes large enough to compensate the decrease in cross-section due to mass suppression. Large $Y$ can also be achieved by increased the VLF mixing $\sin\theta$ and for \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}\mathcal{O}(1) we can satisfy correct abundance in this region even with moderate $\Delta M$. As we go beyond DM mass of 1 TeV, $Z_{B-3L_{\tau}}$ resonance shows up (as the minimum value of $m_{Z_{B-3L_{\tau}}}$ is 2.5 TeV). Now, since both $\tilde{v}$ and $g_{B-3L_{\tau}}$ are being varied, $m_{Z_{B-3L_{\tau}}}$ is not fixed (according to Eq. (14)). As a result, the resonance region is not sharp but broad due to different $m_{Z_{B-3L_{\tau}}}$. For all possible choices of $\tilde{v}$ and $g_{B-3L_{\tau}}$ according to Eq. (41), $m_{Z_{B-3L_{\tau}}}$ is being varied between $\sim 2.5\leavevmode\nobreak\ \rm TeV$ to $\sim 12\leavevmode\nobreak\ \rm TeV$. Because of this, the resonance band lies between $M_{\psi_{1}}\simeq\{\frac{2.5}{2}-\frac{12}{2}\}\leavevmode\nobreak\ \rm TeV$ for all possible $\Delta M$. It is seen that regions corresponding to over-abundance, under-abundance and right relic overlap on each other in the $Z_{B-3L_{\tau}}$ resonance region due to different values of $m_{Z_{B-3L_{\tau}}}$ and $\Delta M$. Note that if we go to even higher DM mass we will still find relic abundance allowed parameter space due to resonances from different $m_{Z_{B-3L_{\tau}}}$.

As seen from the top left plot of Fig. 4, for DM mass $M_{\psi_{1}}\gtrsim 100$ GeV (with $\sin\theta>0.1$, there exists two different $\Delta M$ for same $M_{\psi_{1}}$ which satisfies correct relic density requirement. To understand this, we plot relic density versus $\Delta M$ for different values of $M_{\psi_{1}}$ shown in the bottom left panel of Fig. 4. Here also, the two different values of $\Delta M$, giving correct relic for same DM mass are visible for $M_{\psi_{1}}\gtrsim 100$ GeV. As it is observed, for DM mass say, 100 GeV (black dashed curve) the relic abundance first rises with increase in $\Delta M$. This happens due to the fact that the co-annihilation becomes sub-dominant due to increase in $\Delta M$, which, in turn, reduces the effective annihilation cross-section. At some point right relic density is reached ($\Delta M\sim 10\leavevmode\nobreak\ \rm GeV$) as the annihilation and co-annihilation are just sufficient to produce the correct abundance. After that, the relic abundance becomes more or less constant for a small $\Delta M$ range and then again the abundance starts going downhill as the Yukawa $Y$ becomes large enough making the net annihilation cross-section larger. The cumulative result top right and bottom left panels of Fig. 4 is reflected in the top left corner of the same figure for 0.01\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}\sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.5. For completeness, in the bottom right corner of Fig. 4 we have shown the relic density allowed parameter space for different choices of the new gauge coupling $g_{B-3L_{\tau}}$. As we noticed earlier in Fig. 3, there is no strong dependence of the relic abundance on $g_{B-3L_{\tau}}$. This is evident from this plot as different coloured points (corresponding to different $g_{B-3L_{\tau}}$) are scattered within the allowed region of $\Delta M-M_{\psi_{1}}$ parameter space. In passing we would like to comment that the annihilation to RHN final states is suppressed because of heavy mass of the RHNs and also their contribution to total annihilation cross-section is negligible compared to the SM quarks.

IV.2 Direct detection of dark matter

The presence of the $Z$ and $Z_{B-3L_{\tau}}$ mediated DM-nucleon scattering diagrams highly constrain the parameter space of singlet-doublet model by pushing $\sin\theta$ to a very small value which, in turn, forces $\Delta M$ to be small Bhattacharya:2015qpa . This can be avoided by exploiting the pseudo-Dirac splitting of the VLFs Bhattacharya:2017sml ; Barman:2019tuo . As mentioned in section II, the presence of the VLF singlet Majorana term splits the Dirac states into two pseudo-Dirac states with a mass difference between the two. From Eq. II.2 we see that due to the presence of the Majorana term (generated by the singlet scalar $\Phi$) and mixing between the singlet-doublet fermions, the physical mass eigenstate $\psi_{1}$ splits into two pseudo-Dirac states: $\{\psi_{1}^{i},\psi_{1}^{j}\}$. In such a scenario, interaction of the DM with $Z$ ($Z_{B-3L_{\tau}}$) can be written as Bhattacharya:2017sml ; Barman:2019tuo :

[TABLE]

where $\mathcal{Z}\in\{Z,Z_{B-3L_{\tau}}\}$ and $g_{z}=\frac{g_{L}\sin^{2}\theta}{2c_{W}}$ for SM $Z$ and $g_{z}=\frac{3}{4}g_{B-3L_{\tau}}\sin^{2}\theta$ for $Z_{B-3L_{\tau}}$ mediation, $c_{W}\equiv\cos\theta_{W}$ stands for cosine of the Weinberg angle. Note that, the $\mathcal{Z}$-mediated interaction term is off-diagonal unlike the kinetic terms due to the pseudo-Dirac nature of the VLFs. This results in an inelastic scattering of the DM in which the DM is scattered to an excited state via $\mathcal{Z}$ mediation. As pointed out in TuckerSmith:2001hy , such an inelastic scattering can occur only if the splitting between the two pseudo-Dirac states $\psi_{1}^{i}$ and $\psi_{1}^{j}$ satisfies:

[TABLE]

As computed in Barman:2019tuo , $\delta\sim 100\leavevmode\nobreak\ \rm keV$ can forbid the inelastic scattering mediated by $Z(Z_{B-3L_{\tau}})$ for a DM mass $\sim\mathcal{O}(100\leavevmode\nobreak\ \rm GeV)$ with \beta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}220\leavevmode\nobreak\ \rm km/s. This is to be noted here, the splitting between the two pseudo-Dirac states is so small ($\sim 100\leavevmode\nobreak\ \rm keV$) that it can be ignored in determining the relic abundance of the DM, but has to be taken into account for computing the direct detection cross-section (as emphasised earlier). The mass splitting between $\psi_{1}^{i}$ and $\psi_{1}^{j}$ in terms of our model parameter is given by:

[TABLE]

From Eq. (44) one can put a bound on the VLF mixing and the Yukawa $y_{\chi}$ for $\delta m\sim 100\leavevmode\nobreak\ \rm keV$ such that the heavy neutral gauge boson mediated diagrams are switched off. This is depicted in Fig. 5. As we can see, in order to forbid such inelastic scattering one can choose $y_{\chi}\sim\mathcal{O}(10^{-8})$, then for all small mixing the inelastic scattering can be forbidden. Now, such a choice of scalar VEVs and small Yukawa is not in conflict with the neutrino mass generation. Again, in order to satisfy the direct detection bound, we need to confine \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.5 which is safe even if we consider the conservative bound from Fig. 5 corresponding to $v_{\Phi}=1\leavevmode\nobreak\ \rm TeV$, which anyway we require to keep all the masses within the reach of the ongoing collider experiment. Therefore, for all practical purposes, we can safely ignore the scattering via heavy neutral gauge bosons in order to evade the stringent direct detection exclusion limit.

The direct detection of the VLF DM in this case, therefore, takes place dominantly via the elastic scattering mediated by the scalars ($h_{1,2,3}$). This is depicted in Fig. 6. The spin-independent (SI) direct detection cross section per nucleon is given by Duerr:2015aka :

[TABLE]

where $A$ is the mass number of the target nucleus, $\mu=\frac{M_{\psi_{1}}M_{N}}{M{\psi_{1}}+M_{N}}$ is the DM-nucleus reduced mass and $\left|\mathcal{M}\right|$ is the DM-nucleus amplitude, which reads:

[TABLE]

The effective couplings (with form factors Durr:2015dna ) in Eq. (46) are:

[TABLE]

with

[TABLE]

where $c_{12}=\cos\theta_{12},c_{13}=\cos\theta_{12}$ and $s_{13}=\sin\theta_{13}$ are the scalar mixing angles, defined earlier. The parameter space satisfying right DM relic abundance in comparison to the present bound from direct search experiment is shown in the LHS Fig. 7 for different choices of VLF mixing $\sin\theta$ and gauge coupling $g_{B-3L_{\tau}}$. We have also shown how much of the parameter space is under the infamous neutrino floor Billard:2013qya where it is extremely difficult or even impossible to distinguish DM signal from the SM neutrino background (light grey region). In the LHS of Fig. 7 we see near the $Z$ and Higgs resonance, small and moderate $\sin\theta$’s are allowed by direct search (0.01\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}\sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.3). For M_{\psi_{1}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}100\leavevmode\nobreak\ \rm GeV larger $\sin\theta$’s are also allowed as the direct search cross-section has a suppression from heavy scalars: $\sigma_{SI}\sim\frac{\mu^{2}\sin^{2}\theta}{m_{h_{i}}^{4}}$. Small \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.1 are always allowed by direct search because they produce smaller scattering cross-section, but they are mostly devoured by the neutrino floor as shown by the grey region in Fig. 7. On the RHS of Fig. 7 we see the relic density allowed parameter space that also satisfies direct search bound. In the low DM mass region, specifically near $Z$ and Higgs resonances we can achieve large $\Delta M$ for moderate $\sin\theta$. But if $\Delta M$ becomes too large \mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}500\leavevmode\nobreak\ \rm GeV then one has to resort to small $\sin\theta$ to tame down the Yukawa $Y$ in order to satisfy both relic abundance and direct search limits. For larger DM mass large $\Delta M$ is still allowed near the non-standard scalar resonances $\sim 100\leavevmode\nobreak\ \rm GeV$ and $\sim 150\leavevmode\nobreak\ \rm GeV$. Beyond $\sim 150\leavevmode\nobreak\ \rm GeV$ large $\Delta M$ is achieved with larger $\sin\theta$, while the DM remains still allowed by direct search due to suppression from heavy scalars mentioned earlier. With DM mass \mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}400\leavevmode\nobreak\ \rm GeV the points move towards smaller $\Delta M$ in order to reach right relic exploiting co-annihilation as we have seen earlier in Fig. 3. Beyond 1 TeV DM mass, the direct detection bound becomes weak as the DM mass is large, while because of $Z_{B-3L_{\tau}}$ resonance there is a huge parameter space that satisfy relic abundance. As a result, for M_{\psi_{1}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}1\leavevmode\nobreak\ \rm TeV, almost all of the parameter space is allowed from direct detection for all possible $\Delta M$. We would like to remind here once more that this is the novel feature of the pseudo-Dirac states that this model offers, which helps to achieve larger $\Delta M$ without constraining $\sin\theta$ to a great extent. Larger $\sin\theta$ is required to distinguish this model at the colliders as we shall elaborate in section VI.

V High scale stability and perturbativity

In this section we will discuss discuss the high scale feature of the model. To be specific, here we will constrain the relic density, direct detection satisfied points by applying perturbativity/unitarity and vacuum stability bounds till some high energy scale. For this purpose we need to consider the RG running of associated couplings through $\beta$ functions. We have used PyR@Te 2.0.0 Lyonnet:2016xiz to extract the $\beta$ functions corresponding to the gauge couplings, relevant scalar and fermionic couplings present in the model which are listed in Appendix B. For simplicity, we show only the one loop $\beta$ functions for both SM and BSM parameters in Appendix B while in our numerical calculations, we consider the three loop beta functions for SM particles due to better precision of SM parameters.

The non violation of perturbativity/unitarity conditions (Eq. (27) and Eq. (28)) for various couplings can be assured by analysing their runnings using the $\beta$ functions. In our analysis, some of the Yukawa like couplings ($y_{\chi},\leavevmode\nobreak\ y_{\alpha_{i}},\leavevmode\nobreak\ y_{13},\leavevmode\nobreak\ y_{23}$) are assumed to be very small. Hence they have negligible influence in the RG running of themselves and other parameters. In addition, we also fix the VEV of the singlet scalar fields within TeV range and scalar mixing angles $\lesssim 0.2$. These in turn fix the magnitude of the scalar couplings which are positive and stays below $\sim 0.1$. With these order of magnitude initial values, they are not expected to break the perturbativity conditions at high energy scale. The other two important parameters we consider are $Y$ and $g_{B-3L_{\tau}}$ which play vital role in DM phenomenology as well as in collider analysis. Largeness of these two parameters could destroy the high scale perturbativity/unitarity of the theory. Therefore we will focus on $Y-g_{B-3L_{\tau}}$ plane and see the the bounds coming from the requirement of satisfying perturbativity/unitarity criteria. While focusing on this particular plane of our interest, we also make sure that none of the other parameters violate the above mentioned criteria. First, in left panel of Fig. 8 we show the points (blue colored) in $Y-g_{B-3L_{\tau}}$ plane which satisfy relic, direct detection bounds. Then we constrain the same plane using the perturbativity criteria till $\mu=M_{P}$ (orange colored points), where $M_{P}$ is the Planck scale. It is clear that applying the perturbativity criteria significantly cuts the earlier parameter space with $Y\gtrsim 0.8$ and $g_{B-3L_{\tau}}\gtrsim 0.25$. Similar exercise has been done in $\Delta M-g_{B-3L_{\tau}}$ plane which is shown in right panel of Fig. 8.

Now the conditions of the stability or the boundedness of the scalar potential in the model till some high energy scale along various field directions are provided in Eq. (26). Recall that the EW vacuum stability (stability of Higgs potential) is dictated by the condition $\lambda_{H}>0$. However for more accurate analysis, one should consider the radiatively improved Higgs potential where the one loop correction will be provided by SM fields and other BSM fields. The radiatively corrected one loop effective Higgs potential (at high energies $h\gg v_{d}$) can be written as Casas:1994qy ; Casas:1996aq ,

[TABLE]

with $\lambda_{H}^{\rm eff}=\lambda_{H}^{\rm SM,eff}+\lambda_{H}^{(S,\Phi),\leavevmode\nobreak\ \textrm{eff}}+\lambda_{H}^{(\psi,\chi),\textrm{ eff}}$ where $\lambda_{H}^{\rm SM,eff}$ is the SM contribution to $\lambda_{H}$. The other two contributions $\lambda_{H}^{(S,\Phi),\leavevmode\nobreak\ \textrm{eff}}$ and $\lambda_{H}^{(\psi,\chi),\textrm{ eff}}$ are due to the newly added fields in the present model as provided below Khan:2012zw ; Gonderinger:2012rd .

[TABLE]

where $\Gamma(h)=\int_{m_{t}}^{h}\gamma(\mu)d\textrm{ln}\mu$ and $\gamma(\mu)$ is the anomalous dimension of the Higgs field Buttazzo:2013uya . Note that we have ignored the radiative corrections involving $y_{\chi},y_{\alpha_{i}},y_{13},y_{23}$ as they are fixed to negligibly small values in our analysis. Now with the inclusion of radiative correction to Higgs potential, the stability condition of Higgs vacuum will be modified as $\lambda_{H}^{\rm eff}>0$. The remaining co-positivity conditions in Eq. (26) will determine the boundedness of the scalar potential in different field directions.

We numerically solve the three loop RG equations for all the SM couplings and one loop RG equations for the other relevant BSM couplings in the model from $\mu=m_{t}$ to $M_{P}$ energy scales considering $m_{t}=173.1$ GeV Tanabashi:2018oca , SM Higgs mass $m_{H}=125.09$ GeV Tanabashi:2018oca and strong coupling constant $\alpha_{s}=0.1184$. We also use the initial boundary values of all the SM couplings as provided in Buttazzo:2013uya . The boundary values have been determined at $\mu=m_{t}$ in Buttazzo:2013uya by taking various threshold corrections and mismatch between top pole mass and $\overline{\rm MS}$ renormalised couplings into account. One important point is to note that during the running of couplings, we will ignore the small mass differences between the masses of heavy BSM Higgs bosons and DM particles for the sake of simplicity. The $\beta$ function of $\lambda_{H}$ includes positive contributions from the scalar couplings and negative contributions from fermionic couplings. Therefore, with $y_{t}\sim\mathcal{O}(1)$ in SM, large value of $Y$ could destabilise the EW vacuum. The initial value of $\lambda_{H}$ also gets a positive shift due to the presence of additional scalars in the set up as evident from Eq. (13). The amount of shift depends on the masses of the heavier Higgs bosons and also the corresponding mixing angles. With our choices for them as specified earlier the magnitude of the shift comes out to be $\sim 0.02$. Note that we have also considered all the other scalar couplings positive and $\sim\mathcal{O}(0.1)$ in our analysis. Hence considering small order of magnitude of Yukawa like couplings ($y_{\alpha_{i}},\leavevmode\nobreak\ y_{13},\leavevmode\nobreak\ y_{23},\leavevmode\nobreak\ y_{\chi}$), the BSM scalar couplings are expected to remain positive in their evolution, thus automatically guaranteeing the stability of the total scalar potential in the corresponding field directions (when $\lambda_{H}>0$).

Now we further constrain the $Y-g_{B-L}$ parameter space which is allowed from perturbativity criteria (Fig. 8) along with correct DM related observables using vacuum stability conditions. Before that in Fig. 9, we show running of $\lambda_{H}(\lambda_{H}^{\textrm{eff}})$ for two different DM relic + direct detection + perturbativity bounds satisfying points having $Y\sim$ 0.25 and 0.46 respectively. As it can be seen, for lower value of $Y$ $\lambda_{H}(\lambda_{H}^{\textrm{eff}})$ remains positive throughout its running till $M_{P}$ energy scale thereby establishing the stability of EW vacuum. On the other hand for $Y\sim 0.46$, $\lambda_{H}(\lambda_{H}^{\textrm{eff}})$ crosses zero around $\mu\sim 10^{15}$ GeV and ends with negative value at $\mu=M_{P}$. Hence it is clear that large values of $Y$ are disfavoured in our analysis as it could destabilise the Higgs vacuum. The plots in Fig. 10 also shows that the running of $\lambda_{H}$ and $\lambda_{H}^{\textrm{eff}}$ are similar and they almost merge near the energy scale $\mu=M_{P}$. Finally in Fig. 10, we constrain left panels of Fig. 8, using both perturbativity and the vacuum stability criteria in both $Y-g_{B-L}$ and $\Delta M-g_{B-3L_{\tau}}$ planes. Now when we compare Fig. 8 with Fig. 10, it clearly shows that the upper limit on $Y$ is significantly reduced from 0.8 to 0.3 due to the application of vacuum stability criteria till energy scale $M_{P}$. However upper limit on $g_{B-3L_{\tau}}$ remains more or less unaltered ($\lesssim 0.25$) as it does not have direct role in stability analysis. Similar conclusion can be drawn for $\Delta M$ also. As before, all the points in these plots satisfy DM related bounds.

Before we move on further, let us first choose a few benchmark points (BPs) which we shall be using for the collider study. Note that, all these BPs need to satisfy correct relic abundance, direct detection bound, vacuum stability and perturbativity constraints and on top of that should give rise to $Z_{B-3L_{\tau}}$ mass in correct range. These are enlisted in table 2. We also include the values of relevant EW precision parameters in table 2 for all the benchmark points which show they fall within correct experimental range. Another point is to note that these BPs are selected in the decreasing order to $g_{B-3L_{\tau}}$ from top to

bottom where BP1 has highest $g_{B-3L_{\tau}}$ and BP5 has the smallest $g_{B-3L_{\tau}}$. As we shall see in section VI, the production cross-section of $\psi^{\pm}$ will be large for small $\Delta M$ and not for large $g_{B-3L_{\tau}}$. This is due to the fact that larger $g_{B-3L_{\tau}}$ results in heavier $m_{Z_{B-3L_{\tau}}}$ (for fixed VEV), which, in turn causes propagator suppression for $\psi^{\pm}$ production (via $Z_{B-3L_{\tau}}$) leading to decrease in cross-section. However, even if $\Delta M$ is small, but $M_{\psi_{1}}$ is large, the production cross-section may still be small. All the BPs satisfy the invisible SM Higgs and SM $Z$ decay constraint as shown in the Appendix. A. Finally we would like to highlight that LEP has set a lower limit on pair-produced charged heavy vector-like leptons: $m_{L}>101.2\leavevmode\nobreak\ \rm GeV$ at 95 % C.L. for $L^{\pm}=\nu W$ final states Achard:2001qw . Thus all our benchmark points are safe from LEP bounds.

VI Collider phenomenology

The detailed study of collider signature for vector like fermions can be found in Bahrami:2016has ; Barman:2019tuo . As we have already seen, due to the pseudo-Dirac nature of the VLFs large $\Delta M$ can be achieved satisfying both relic abundance and direct search. Such large $\Delta M$’s are actually beneficial in order to distinguish this model at the collider from the SM background Barman:2019tuo . It is to be noted that the charged component of $SU(2)_{L}$ doublet VLF can be produced at the LHC via SM $Z$, $Z_{B-3L_{\tau}}$ and photon mediation. The charged VLF can further decay via on-shell and/or off-shell $W$ (depending on whether \Delta M\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}80\leavevmode\nobreak\ \rm GeV or \Delta M\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}80\leavevmode\nobreak\ \rm GeV) to the following final states:

• •

Hadronically quiet opposite sign dilepton (OSD) with missing energy $\left(\ell^{+}\ell^{-}+\not{E_{T}}\right)$.

• •

Single lepton, with two jets plus missing energy $\left(\ell^{\pm}+jj+\not{E_{T}}\right)$.

• •

Four jets plus missing energy $\left(jjjj+\not{E_{T}}\right)$.

We shall focus only on the leptonic final states as they are much cleaner compared to others, the Feynman diagram for which is depicted in Fig. 11. To be more specific, we shall only look into the hadronically quiet dilepton final states as we are interested to see how the presence of $Z_{B-3L_{\tau}}$ can affect the coillder signatures compared to purely $Z$ mediated scenarios studied earlier. The presence of $Z_{B-3L_{\tau}}$ significantly increases the production cross-section of the charged VLF pairs at the collider. This is due to the fact that the decay of the BSM neutral gauge boson to the charged VLFs now happens on-shell in contrast to models where this decay takes place off-shell via SM $Z$ boson and photon Bhattacharya:2017sml ; Barman:2019tuo . Also, as there is no negative interference between the $Z$ and $Z_{B-3L_{\tau}}$ mediated charged VLF production channels, hence the addition of new channel always improves the production cross-section. In order to illustrate this improvement, we first show the variation of production cross-section $\sigma_{pp\to\psi^{+}\psi^{-}}$ in the LHS of Fig. 12 for two different choices of $\{g_{B-3L_{\tau}},\tilde{v}\}$. One noteworthy feature of this plot is that the production cross-section is lower for the choice $\{g_{B-3L_{\tau}},\tilde{v}\}=\{0.2,1.8\}$ (black dotted curve) than for $\{g_{B-3L_{\tau}},\tilde{v}\}=\{0.3,0.8\}$ (black dashed curve). This is simply attributed to the propagator suppression due to larger mass of $Z_{B-3L_{\tau}}$ in the former case ($m_{Z_{B-3L_{\tau}}}=3.36\leavevmode\nobreak\ \rm TeV$) over the latter ($m_{Z_{B-3L_{\tau}}}=2.54\leavevmode\nobreak\ \rm TeV$). On the RHS of Fig. 12 we have illustrated how the production cross-section changes for different choices of $g_{B-3L_{\tau}}$ and $\tilde{v}$ (and hence $m_{Z_{B-3L_{\tau}}}$) keeping $v_{\Phi}$ fixed at 3 TeV when both $Z$ and $Z_{B-3L_{\tau}}$ mediations are present. $g_{B-3L_{\tau}}$ and $\tilde{v}$ are chosen in such a way that $m_{Z_{B-3L_{\tau}}}$ is always above the LHC lower bound. In the same plot we have also shown our chosen BPs appearing in table 2. Note that, the production cross-section for BP2 is the least, while it is highest for BP1 and BP4 (overlapped on each other). This tells the fact that though large $\Delta M$ is necessary in distinguishing the model at the collider (as we shall see) but at the same time we need to compromise with the production cross-section. Again, the production cross section for M_{\psi^{\pm}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}800\leavevmode\nobreak\ \rm GeV is either very small (kinematically) or discarded by stability and perturbativity bound on $Y$ as larger $M_{\psi^{\pm}}$ requires larger $\Delta M$, which in turn makes $Y$ large and that is constrained from Fig. 10. Although $Y$ can be tamed down by choosing a small $\sin\theta$ as per Eq. (22) but the production cross-section will still remain small. Therefore, we have overlooked all such benchmarks. In both the plots the production cross-section decrease with the increase in charged VLF mass showing the unitarity of the cross-section remains valid.

VI.1 Object reconstruction and simulation details

As already mentioned, we implemented this model in LanHEP and the parton level events are generated in CalcHEP Belyaev:2012qa . Those events are then passed through PYTHIA Sjostrand:2006za for showering and hadronisation. All the SM backgrounds that can mimic our final state are generated in MADGRAPH Alwall:2011uj and the corresponding production cross-sections are multiplied with appropriate $K$-factor Alwall:2011uj in order to match with the next to leading order (NLO) cross-sections. For all cases we have used CTEQ6l as the parton distribution function (PDF) Placakyte:2011az . Now, in order to re-create the collider environment, all the leptons, jets and unclustered objects have been reconstructed using the following set of criteria:

• •

Lepton ($l=e,\mu$): Leptons are identified with a minimum transverse momentum $p_{T}>20$ GeV and pseudorapidity $|\eta|<2.5$. Two leptons can be distinguished separately if their mutual distance in the $\eta-\phi$ plane is $\Delta R=\sqrt{\left(\Delta\eta\right)^{2}+\left(\Delta\phi\right)^{2}}\geq 0.2$, while the separation between a lepton and a jet needs to be $\Delta R\geq 0.4$.

• •

Jets ($j$): All the partons within $\Delta R=0.4$ from the jet initiator cell are included to form the jets using the cone jet algorithm PYCELL built in PYTHIA. We demand $p_{T}>20$ GeV for a clustered object to be considered as jet. Jets are isolated from unclustered objects if $\Delta R>0.4$.

• •

Unclustered Objects: All the final state objects which are neither clustered to form jets, nor identified as leptons, belong to this category. Particles with $0.5<p_{T}<20$ GeV and $|\eta|<5$ are considered as unclustered. Although unclustered objects do not interfere with our signal definition but they are important in constructing the missing energy of the events.

• •

Missing Energy ($\not{E}_{T}$): The transverse momentum of all the missing particles (those are not registered in the detector) can be estimated from the momentum imbalance in the transverse direction associated to the visible particles. Missing energy (MET) is thus defined as:

[TABLE]

where the sum runs over all visible objects that include the leptons, jets and the unclustered components.

• •

Invariant dilepton mass $\left(m_{\ell\ell}\right)$: We can construct the invariant dilepton mass variable for two opposite sign leptons by defining:

[TABLE]

Invariant mass of OSD events, if created from a single parent, peak at the parent mass, for example, $Z$ boson. As the signal events (Fig. 11) do not arise from a single parent particle, invariant mass cut plays key role in eliminating the $Z$ mediated SM background.

• •

$H_{T}$: $H_{T}$ is defined as the scalar sum of all isolated jets and lepton $p_{T}$’s:

[TABLE]

For our signal the sum only includes the two leptons that are present in the final state.

We shall use different cuts on these observables to separate the signal from the SM backgrounds and predict the significance as a function of the integrated luminosity. This is shown in the next section.

VI.2 Event rates and signal significance

Here we would first like to show how the presence of new charge neutral gauge boson mediation can affect the pair production cross-section of the charged VLFs. This is explicitly tabulated in table 3 where we have listed the production cross-sections for our chosen BPs (table 2) both in the presence and in the absence of $Z_{B-3L_{\tau}}$. As expected, in each case, the production via $Z$ and $Z_{B-3L_{\tau}}$ together is larger than that of only $Z$ mediation. The improvement, however, is not significant enough due to the reasons mentioned earlier.

In Fig. 13 we have shown the distribution of normalised number of events with respect to MET (upper panel) and $H_{T}$ (lower panel) for all the chosen BPs. In the same plot we have also shown the distribution from dominant SM backgrounds that can mimic our signal. For the SM the only source of MET are the SM neutrinos, which are almost massless with respect to centre of mass energy of the collider. As a result, the MET and $H_{T}$ distribution for SM peaks up at a lower value, while for the model MET arises from the DM $\psi_{1}$ (on top of the SM neutrinos) which is massive, and hence corresponding distribution for the signals are much flattened. Noteworthy feature here is that, for larger $\Delta M$ the signal distributions are well separated from that of the background. This is due to the fact that the peak of the MET distribution is determined by how much of $p_{T}$ is being carried away by the missing particle (i.e, the DM), which in turn depends on the mass difference of charged and neutral component of the VLF i.e, $\Delta M$. Hence for larger $\Delta M$ the DM carries away most of the $p_{T}$ making the distribution much flatter, while for smaller $\Delta M$ the distribution peaks up at lower value as the produced DM particles are not boosted enough. As a consequence, in the LHS of top left panel of Fig. 13 we see BP3 (in blue) is completely submerged in the SM background, while BP1 (in red) can still be distinguished to some extent. BP2, because of large $\Delta M$ has a rather flattened distribution (in green) and therefore can be easily distinguished from the background with judicious choice of cuts. On the top right panel of Fig. 13 we have shown the MET distribution for BP4 (red) and BP5 (green). Here we also see the same consequence: with comparatively larger $\Delta M$ BP4 can be separated from the background, while BP5 shows no excess over the SM background. This trend is similar for $H_{T}$ distribution, which we have shown in the bottom panel of Fig. 13.

From the distributions one can easily see, with a MET cut of \not{E_{T}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}200\leavevmode\nobreak\ \rm GeV and a $H_{T}$ cut of \mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle>}}{{\sim}}\leavevmode\nobreak\ }}250\leavevmode\nobreak\ \rm GeV one can get rid off the SM backgrounds keeping most of the signal events intact. This is also shown in table 4 where we have demonstrated the cut-flow i.e, how the number of events vary with MET cut, while the $H_{T}$ cut is kept fixed $H_{T}>250\leavevmode\nobreak\ \rm GeV$. On top of that we have also imposed the invariant mass cut over $Z$-window such that no events lie in the range: |m_{Z}-15|\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}m_{\ell^{+}\ell^{-}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}|m_{Z}+15| in order to reduce the SM $Z$ background as explained earlier. Corresponding cut-flow for dominant SM backgrounds are also tabulated in table 5. As one can see, the $t\bar{t}$ background is completely killed by imposing zero jet veto. Amongst other backgrounds, $ZZ$ also vanishes because of imposition of the $m_{\ell\ell}$ cut and $WW$ is also killed by putting a hard MET cut of 300 GeV. The only remaining background is due to $WWZ$ but that also becomes insignificant due to large MET cut.

We finally plot the signal significance for the BPs in Fig. 14 by choosing the final state events with $\not{E_{T}}>300\leavevmode\nobreak\ \rm GeV$ and $H_{T}>250\leavevmode\nobreak\ \rm GeV$. BP2 and BP4 can reach $5\sigma$ discovery for an integrated luminosity $\sim 200\leavevmode\nobreak\ \rm fb^{-1}$ as they have the advantage of large $\Delta M$ which helps them to distinguish from the SM background as explained earlier. Due to comparatively smaller $\Delta M$ BP1 can reach a 5$\sigma$ discovery at a slightly higher luminosity $\sim 500\leavevmode\nobreak\ \rm fb^{-1}$. BP3 and BP5 can only be probed at the very high luminosity (HL-LHC). Here we would also like to emphasize the fact that this model may also be probed at the collider via stable charged track signature for \Delta M\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}80\leavevmode\nobreak\ \rm GeV. In that case the decay of the charged VLFs happen via off-shell $W^{\pm}$ and the decay width can be small enough for small $\sin\theta$ giving rise to charged tracks of length $\sim\mathcal{O}(\rm cm)$. This has been explored in details in Bhattacharya:2017sml ; Barman:2019tuo and hence we refrain from discussing it again here.

VII Summary and conclusion

We have proposed a flavoured gauge extension of the singlet-doublet fermionic dark matter model by considering $B-3L_{\tau}$ as the additional gauge quantum number which naturally stabilises the DM without the need of additional discrete symmetries. The model also requires the existence of a singlet right handed neutrino (RHN) in order to be anomaly free. This RHN, along with another one or two singlet RHNs (having zero $B-3L_{\tau}$ charges) can take part in generating light neutrino masses via type I seesaw mechanism. The neutrino sector and the DM sector, however, are not very closely related as the bounds on the VEVs of the non-standard scalars from correct neutrino mass requirement is rather lose. This is attributed to the fact that light neutrino mass in the right ballpark can be generated keeping the scalar VEVs $\sim\mathcal{O}(\rm TeV)$ scale by tuning the new Yukawa couplings $(y_{\ell},y_{\tau_{3}})$ accordingly. The family non-universal nature of this $B-3L_{\tau}$ gauge symmetry also helps us to avoid strong bounds from the LHC searches. The relatively lighter $Z^{\prime}$ boson plays a role in generating dark matter relic abundance, leading to an enlarged parameter space satisfying DM related constraints compared to the purely singlet-doublet model. The pseudo-Dirac nature of the DM forbids the inelastic scattering via heavy neutral gauge boson allowing the DM to live over a huge parameter space satisfying both relic abundance and direct search constraints. Thus the parameter space remians valid upto DM mass of a few TeV for singlet-doublet VLF mixing of \sin\theta\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.5. A substantial portion of the parameter space however merges with the neutrino floor for smaller $\sin\theta$.

Apart from the motivations from dark matter and neutrino mass generation, the model also provides a solution to the electroweak vacuum metastability problem due to extended scalar sector. The model not only gives rise to a stable electroweak vacuum but also keeps it perturbative all the way upto the Planck scale. The requirements of vacuum stability and perturbativity however, significantly constrains the parameter space allowed purely from dark matter related constraints. As an effect of cumulative bound from relic abundance, direct search and stability and perturbativity of the scalar potential, the model substantially constraints the singlet-doublet Yukawa Y\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.3 and the gauge coupling g_{B-3L_{\tau}}\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}0.25. However, the mass difference between the heavier and lighter physical states of the VLF: $M_{\psi_{2}}-M_{\psi_{1}}=\Delta M$ can still be large enough $\sim 500\leavevmode\nobreak\ \rm GeV$ providing opportunity for the model to be probed at the LHC via hadronically quiet dilepton final states with missing energy excess. This is again attributed to the pseudo-Dirac nature of the VLFs due to which larger $\sin\theta$ is allowed from direct search, and hence large $\Delta M$ is possible to achieve.

We finally discussed possible signatures at colliders by analysing some of the benchmark points of the model which satisfy all theoretical and experimental bounds. We particularly focus on purely leptonic final states with missing energy and show that with judicious choice of cuts on different kinematical variables (eg. MET, $H_{T}$ etc) it is indeed possible to attain a $5\sigma$ discovery potential for the model at the high luminosity LHC. Apart from leptonic final states, the model may also be probed via displaced vertex signature due to the off-shell decay of the charged VLF to SM leptons and neutrino for \Delta M\mbox{\raisebox{-2.58334pt}{\leavevmode\nobreak\ \stackrel{{\scriptstyle<}}{{\sim}}\leavevmode\nobreak\ }}m_{W}.

Acknowledgements.

DB acknowledges the support from IIT Guwahati start-up grant (reference number: xPHYSUGI-ITG01152xxDB001), Early Career Research Award from DST-SERB, Government of India (reference number: ECR/2017/001873) and Associateship Programme of Inter University Centre for Astronomy and Astrophysics (IUCAA), Pune. BB and PG would like to thank Triparno Bandyopadhyay and Subhaditya Bhattacharya for useful discussions during very early stage of this work. BB would also like to thank Krishnanjan Pramanik for computational helps.

Appendix A Invisible Higgs and $Z$ decays

The combination of SM channels yields an observed (expected) upper limit on the SM Higgs branching fraction of 0.24 at 95 % CL Khachatryan:2016whc with a total decay width $\Gamma=4.07\times 10^{-3}\leavevmode\nobreak\ \rm GeV$. This gives rise to an allowed Higgs invisible decay branching fraction of 0.24(0.23). SM $Z$ boson, on the other hand, can also decay to invisible final states, and hence constrained from observation: $\Gamma_{inv}^{Z}=499\pm 1.5\leavevmode\nobreak\ \rm MeV$ PhysRevD.98.030001 . So, if $Z$ is allowed to decay invisibly, the decay width should not be more than 1.5 MeV. In our case, the decays of SM Higgs and SM $Z$ are only possible to $\psi_{1}\psi_{1}$ pairs as other invisible decay modes are kinematically forbidden because of large $\Delta M$. These decay widths are given by:

[TABLE]

and

[TABLE]

with $c_{W}=m_{W}/m_{Z}$ is the Weinberg angle where $m_{W(Z)}$ is the mass of SM $W(Z)$ boson.

Now, for our chosen BPs we would like to see whether these bounds are applicable or not. First note that, Higgs invisible decay is possible for BP2 and BP4, while $Z$ can decay invisibly to DM pairs only for BP4. Rest of the benchmarks are safe from such bounds as DM mass is much above than SM Higgs or SM $Z$ mass. In table 6 we have tabulated the invisible branching fraction ($Br_{inv}^{higgs}$) for SM Higgs (left column) and invisible $Z$ decay width $\Gamma_{inv}^{Z}$ for SM $Z$ (right column).

Appendix B RG equations at one loop

[TABLE]

Appendix C Unitarity

Quartic couplings of the scalar potential which are shown in Eq. (25) are also constrained from tree level perturbative unitarity. The unitarity bounds are related with scattering amplitude as: $|\mathcal{M}|\leq 8\pi$ . $\mathcal{M}$ be the scattering amplitude for any $2=2$ process which can be expressed in terms of partial waves as follows:

[TABLE]

where $P_{l}(\cos\theta)$ is the Legendre polynomial of order $l$, $a_{l}$ be the partial wave amplitude and $\theta$ be the scattering angle. To implement unitarity bound in our case we form an amplitude matrix $M={\mathcal{M}}_{i=j}$ where $i$ and $j$ correspond to all possible two particle state. And each eigenvalue of this amplitude matrix, $M$ should lie within $8\pi$ (i.e. $|e_{i}|\leq 8\pi$)in order to maintain unitarity. The amplitude matrix $M$ is decomposed of 22 neutral charged (NC) and 6 singly charged(SC) two particles state which is given by:

[TABLE]

The charged neutral two particles staTes which are formed the sub-matrix, $M^{NC}$ are given by:

[TABLE]

[TABLE]

And the singly charged two particle states for the sub-matrix $M^{SC}$ are as follows :

[TABLE]

[TABLE]

Each of distinct eigenvalues of the amplitude matrix, $M$ will be bounded from tree level unitarity as :

[TABLE]

where, $x_{1,2,3}$ are the cubic roots of the following polynomial equation:

[TABLE]

Appendix D Relevant Feynmann Diagrams for DM (co-)annihilation