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

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.


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 [1] in 1933, observations of galaxy rotation curves in 1970's by Rubin and collaborators [2], the observation of the bullet cluster by Chandra observatory [3] and the measurements of cosmic microwave background (CMB) by several cosmology experiments, the latest of which is the Planck experiment [4]. In terms of density parameter Ω DM and h = Hubble Parameter/(100 km s −1 Mpc −1 ), the present DM abundance is conventionally reported as [4]: Ω DM h 2 = 0.120 ± 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 [5].
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 [6]. 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 [7] (LLHH)/Λ, L ≡ lepton doublet, Λ ≡ 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 [8][9][10] 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 vanishing B − 3L τ 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 [54]. 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 [55][56][57][58][59][60][61][62]. 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 2  has non-zero charge under U (1) B−3Lτ . Apart from cancelling the triangle anomalies arising due to the U (1) B−3Lτ symmetry (shown below) these right handed neutrinos help us to generate light 2 With VLF alone, it is possible to make the electroweak vacuum stable, if they have coloured charges, as shown in [63].
neutrino masses via type I seesaw mechanism as we shall discuss in detail. We also have one VLF doublet ψ T : ψ 0 ψ − and one VLF singlet χ. 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 The condition for anomaly cancellation comes from the triangle diagrams with gauge gauge bosons at the vertices: 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: where x i stands for the U (1) B−3Lτ charges for all the particles appearing in table I. The other are trivially cancelled. The U (1) B−3Lτ charges of the relevant fermions are assigned in table I 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τ 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 [30][31][32][33] 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 [65][66][67][68][69]. 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: The gauge part of the Lagrangian is written as: with The first term in Eq. (4) is the kinetic term for the new U (1) B−3Lτ gauge boson and the second term arises due to kinetic mixing between U (1) Y and U (1) B−3Lτ gauge bosons. For simplicity we choose = 0 at the scale of B − 3L τ symmetry breaking as in [70]. The Lagrangian for the new fermion sector (excluding their couplings to the SM Higgs) reads: where under SM ⊗ U (1) B−3Lτ , the covariant derivative is defined as: 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τ corresponds to B − 3L τ 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 τ charge of the scalar singlet S (+3) as N 3 must have B − 3L τ 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: where With one SM Higgs doublet and two non-standard singlets the renormalisable scalar potential takes the form: 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 The gauged U (1) B−3Lτ 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: where U (θ 12 , θ 13 , θ 23 ) is usual unitary matrix involving the mixing angles θ 12 , θ 13 , θ 23 and complex phase δ = 0. We assume h 1 to be SM-like Higgs with a mass of 125 GeV. The minimisation conditions are given by: The weak eigenstates (h, φ, s) in terms of the physical eigenstates (h 1 , h 2 , h 3 ) and the mixing angles are given by: , where c ij = cos θ ij , s ij = sin θ ij and t ij = tan θ ij with {i, j} = 1, 2, 3 and i = j. It is easy to understand from here, for θ 12 , θ 13 = 0, θ 23 = 0 we revive the SM Higgs purely from h 1 . Therefore, θ 12 and θ 13 are constrained from Higgs data, while θ 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 τ , 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 GeV.
Also, the mixing angles θ 12 and θ 13 are constrained from Higgs data which we shall elaborate in section III. Essentially the scalar sector has the following free parameters: All the quartic couplings appearing in the scalar potential can be expressed in terms of the physical masses and mixings as follows: After SSB we are also left with a new massive charge neutral gauge bosons corresponding to broken U (1) B−3Lτ . The mass of the new gauge boson is given by: 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τ 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.

A. VLF mass eigenstates
The Dirac mass matrix in the basis {χ, ψ 0 }, containing the singlet and doublet VLF can be written as: where Thus, the physical eigenstates arise as: where θ is the mixing angle given by: The mass for the physical eigenstates are: The lightest physical charge neutral fermion from above is a viable DM candidate in this model and we choose it to be ψ 1 (M ψ 1 < M ψ 2 ). The DM is naturally stable due to our particular choice of the U (1) B−3Lτ charge. In the small mixing limit the charged component of the VLF doublet ψ ± acquires a mass: For small sin θ (≈ 0) limit, M ψ ± M ψ 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: B. Pseudo-Dirac mass splitting Due to presence of the Majorana term : y χ χ c χΦ, the pseudo-Dirac mass matrix for the VLF singlet χ can be expressed in the basis (χ c χ) T : The mass matrix can be expressed in terms of physical states as: where (χ a χ b ) T is the physical pseudo-Dirac eigenstate. Since m χ << M χ , which we can always assume the last equality in above equation by taking either y χ or v Φ or both to be small. In the presence of the Majorana term, the singlet χ is split into two pseudo-Dirac states, {χ a , χ b } which are propagated into the physical states {ψ a,b 1 , ψ a,b 2 } via VLF mixing. As these two pseudo-Dirac states (ψ a,b ) are nearly degenerate i.e, δm ∼ O (100keV), we consider them to be a single state ψ = ψ a ψ b T with mass M ψ , identical to the Dirac mass of ψ. 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 : where ∆M = M ψ 2 − M ψ 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τ 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.
A. Stability, perturbativity and tree-level unitarity

Stability
Stability of the scalar potential is mainly dictated by the quartic terms of the scalar potential, V (H, Φ, S) which is defined as: 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 [71,72]:

Perturbativity
To prevent perturbative breakdown of the model, all quartic, Yukawa and gauge couplings should obey the following limits at any energy scale: where j = 1, 2 and α = e, µ.

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 [73][74][75]: where, x 1,2,3 are the cubic roots of the polynomial equation detailed in Appendix C.
Here we would like to estimate the effect of the BSM particles on those parameters. We have four parameters, namelyŜ,T , W and Y [80] whereŜ is related to the Peskin-Takeuchi parameter S: T parameter is more significant for small mixing in the scalar sector and the constraint on Tparameter is parametrised by the following data [6]: ∆T = T xSM − T SM = 0.07 ± 0.12, where we consider contributions from both the VLFs and the non-standard scalars to T xSM . In this situation ∆T is given by [81,82]: which indicates how much the oblique parameter is shifted from the SM value. Now, and T SM Higgs (h 1 ) can be obtained by using the decoupling limits s 12 → 0 and s 13 → 0 in Eq. (30). The contribution from VLF DM is followed as, where 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 Π'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: 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 ψ , M ψ 1 and M ψ 2 (M ψ ≡ M ψ ± according to Eq. (21)).
Once more we would like to remind that M ψ 2 ≈ M ψ ± under small mixing limit as apparent from Eq. (21).
The bound onŜ comes from a global fit: 10 3Ŝ = 0.0 ± 1.3 [80]. For S-parameter, we consider contribution only due to the VLFs 4 as given by [21,79]: where g 2 is the SU (2) L gauge coupling. The expression for vacuum polarization for identical masses (at q 2 = 0) [79]:Π For two different masses (m i = m j ) the expression for vacuum polarization reads [79]: corresponding to red, green and blue contours respectively. As it is seen, larger v Φ 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ṽ for a fixed v Φ . As a result, for the parameter space scan we have varied both g B−3Lτ andṽ for a fixed v Φ to keep the Z B−3Lτ mass in the right ballpark ( > ∼ 2.5 TeV). Here we would like to mention that a combination of B s -B s mixing from flavour physics, together with ATLAS' Z search puts a bound on third family hypercharge models by requiring m Z > 1.9 TeV [85]. However, as we are considering even more conservative bound, our models is safe from such constraints arising from flavour physics measurements.

D. 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 [86,87]. In case of scalar singlet extension of SM, the strongest bound on scalar-SM Higgs mixing angle (θ m ) comes form W boson mass correction [88] at NLO for 250 GeV m h 2 850 GeV as  [89,90] and measured Higgs signal strength [90] restrict the mixing angle sin θ m dominantly ( 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 θ 12 , sin θ 13 ∼ 0.1. We also fix v Φ > ∼ 1 TeV which helps in keeping the perturbativity of the theory intact. The other mixing angle sin θ 23 is a free parameter. We keep it below 0.2.

E. Neutrino mass
The light neutrino mass matrix can be generated via type I seesaw mechanism where 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 [91]. This, in turn, puts bound on model parameters including the VEV v 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 × 3 PMNS matrix (choosing the charged lepton mass matrix diagonal) gives the light neutrino masses. We can choose: and y l y τ 3 ∼ O(10 −7 ), then we can produce correct order of light neutrino mass for y ∼ 0.  [92]. 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, ψ 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 τ charge assignment (discussed in the subsection II A). In this section we have explored in detail the parameter space appearing in Eq. (24) allowed by observed relic abundance (Ωh 2 = 0.120 ± 0.001 [4]) and direct search limits (particularly from the XENON 1T experiment [29]).

A. Relic abundance of the DM
Relic density of DM, ψ 1 is governed by SM Higgs (h 1 ) and heavy Higgs (h 2 , h 3 ) mediated annihilation and co-annihilation types number changing processes along with SM gauge boson (Z, W ± , γ) and additional heavy gauge boson (Z B−3Lτ ) 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 [5] of the form: with n = n ψ 1 and H being the Hubble parameter in radiation dominated universe. All types of DM number changing processes are taken into account inside σv ef f [93,94] which is given by In above equation, g ef f is defined as the effective degrees of freedom, given by whereḡ 1 ,ḡ 2 andḡ 3 are the internal degrees of freedom of ψ 1 , ψ 2 and ψ ± respectively, and x = T f , where T f is the freeze out temperature. Relic density of ψ 1 one can approximately expressed as [95]: where x f ≈ 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 (ψ 1 ) relic density. In order to calculate relic density, we have used the package MicrOmegas [96] for which the model files are generated from LanHEP [97].
To interesting to note from this plot that as we increase ∆M from 100 GeV to 500 GeV, the relic abundance decreases for M ψ 1 50 GeV. This is due to the fact that, although increased ∆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  This clearly tells the fact that the dependence of DM relic abundance on the new gauge coupling g B−3Lτ is mild compared to the dependence on the VLF mixing sin θ and mass difference ∆M . In each of the plots the dashed red straight line corresponds to the central value of Planck limits on DM relic abundance [4].
We now scan all the free parameters of our analysis in the following range: points (corresponding to different g B−3Lτ ) are scattered within the allowed region of ∆M − M ψ 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.

B. Direct detection of dark matter
The presence of the Z and Z B−3Lτ mediated DM-nucleon scattering diagrams highly constrain the parameter space of singlet-doublet model by pushing sin θ to a very small value which, in turn, forces ∆M to be small [19]. This can be avoided by exploiting the pseudo-Dirac splitting of the VLFs [20,98]. 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. 24 we see that due to the presence of the Majorana term (generated by the singlet scalar Φ) and mixing between the singlet-doublet fermions, the physical mass eigenstate ψ 1 splits into two pseudo-Dirac states: {ψ i 1 , ψ j 1 }. In such a scenario, interaction of the DM with Z (Z B−3Lτ ) can be written as [20,98]: where Z ∈ {Z, Z B−3Lτ } and g z = g L sin 2 θ 2c W for SM Z and g z = 3 4 g B−3Lτ sin 2 θ for Z B−3Lτ mediation, c W ≡ cos θ W stands for cosine of the Weinberg angle. Note that, the 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 Z mediation. As pointed out in [99], such an inelastic scattering can occur only if the splitting between the two pseudo-Dirac states ψ i 1 and ψ j 1 satisfies: As computed in [98], δ ∼ 100 keV can forbid the inelastic scattering mediated by Z(Z B−3Lτ ) for a DM mass ∼ O(100 GeV) with β < ∼ 220 km/s. This is to be noted here, the splitting between the two pseudo-Dirac states is so small (∼ 100 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 ψ i 1 and ψ j 1 in terms of our model parameter is given by: From Eq. (44) one can put a bound on the VLF mixing and the Yukawa y χ for δm ∼ 100 keV such that the heavy neutral gauge boson mediated diagrams are switched off. This is depicted in Each colored region is where inelastic scattering gets disallowed.
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 θ < ∼ 0.5 which is safe even if we consider the conservative bound from Fig. 5 corresponding to v Φ = 1 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 [100]: where A is the mass number of the target nucleus, µ = M ψ 1 M N M ψ 1 +M N is the DM-nucleus reduced mass and |M| is the DM-nucleus amplitude, which reads: with where c 12 = cos θ 12 , c 13 = cos θ 12 and s 13 = sin θ 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 θ and gauge coupling g B−3Lτ . We have also shown how much of the parameter space is under the infamous neutrino floor [102] 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 θ's are allowed by direct search (0.01 < ∼ sin θ < ∼ 0.3). For ∆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 ∆M without constraining sin θ to a great extent. Larger sin θ 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 β functions. We have used PyR@Te 2.0.0 [103] to 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 v d ) can be written as [104,105], where Γ(h) = h mt γ(µ)dlnµ and γ(µ) is the anomalous dimension of the Higgs field [59]. Note that we have ignored the radiative corrections involving y χ , y α 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 λ eff H > 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 µ = m t to M P energy scales considering m t = 173.1 GeV [6], SM Higgs mass m H = 125.09 GeV [6] and strong coupling constant α s = 0.1184. We also use the initial boundary values of all the SM couplings as provided in [59]. The boundary values have been determined at µ = m t in [59] by taking various threshold corrections and mismatch between top pole mass and 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 β function of λ H includes positive contributions from the scalar couplings and negative contributions from fermionic couplings. Therefore, with y t ∼ O(1) in SM, large value of 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 λ H (λ eff H ) for two different DM relic + direct detection + perturbativity bounds satisfying points having Y ∼ 0.25 and 0.46 respectively. As it can be seen, for lower value of Y λ H (λ eff H ) remains positive throughout its running till M P energy scale thereby establishing the stability of EW vacuum. On the other hand for Y ∼ 0.46, λ H (λ eff H ) crosses zero around µ ∼ 10 15 GeV and ends with negative value at µ = 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 λ H and λ eff H are similar and they almost merge near the energy scale µ = 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 ∆M − g B−3Lτ 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τ remains more or less unaltered ( 0.25) as it does not have direct role in stability analysis.
Similar conclusion can be drawn for ∆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τ mass in correct range. These are enlisted in table II. We also include the values of relevant EW precision parameters in table II 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τ from top to bottom where BP1 has highest g B−3Lτ and BP5 has the smallest g B−3Lτ . As we shall see in section VI, the production cross-section of ψ ± will be large for small ∆M and not for large g B−3Lτ . This is due to the fact that larger g B−3Lτ results in heavier m Z B−3Lτ (for fixed VEV), which, in turn causes propagator suppression for ψ ± production (via Z B−3Lτ ) leading to decrease in cross-section. However, even if ∆M is small, but M ψ 1 is large, the is also quoted.
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 GeV at 95 % C.L.
for L ± = νW final states [108]. Thus all our benchmark points are safe from LEP bounds. The detailed study of collider signature for vector like fermions can be found in [98,109]. As we have already seen, due to the pseudo-Dirac nature of the VLFs large ∆M can be achieved satisfying both relic abundance and direct search. Such large ∆M 's are actually beneficial in order to distinguish this model at the collider from the SM background [98]. 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τ and photon mediation. The charged VLF can further decay via on-shell and/or off-shell W (depending on whether ∆M > ∼ 80 GeV or ∆M < ∼ 80 GeV) to the following final states: • Hadronically quiet opposite sign dilepton (OSD) with missing energy ( + − + / E T ).
• Four jets plus missing energy (jjjj + / E T ). 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τ can affect the coillder signatures compared to purely Z mediated scenarios studied earlier. The presence of Z B−3Lτ 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 [20,98]. Also, as there is no negative interference between the Z and Z B−3Lτ 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 σ pp→ψ + ψ − in the LHS of Fig. 12  This tells the fact that though large ∆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 ψ ± > ∼ 800 GeV is either very small (kinematically) or discarded by stability and perturbativity bound on Y as larger M ψ ± requires larger ∆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 θ 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.

A. Object reconstruction and simulation details
As already mentioned, we implemented this model in LanHEP and the parton level events are generated in CalcHEP [110]. Those events are then passed through PYTHIA [111] for showering and hadronisation. All the SM backgrounds that can mimic our final state are generated in MADGRAPH [112] and the corresponding production cross-sections are multiplied with appropriate K-factor [112] 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) [113]. 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, µ): Leptons are identified with a minimum transverse momentum p T > 20 GeV and pseudorapidity |η| < 2.5. Two leptons can be distinguished separately if their mutual distance in the η − φ plane is ∆R = (∆η) 2 + (∆φ) 2 ≥ 0.2, while the separation between a lepton and a jet needs to be ∆R ≥ 0.4.
where the sum runs over all visible objects that include the leptons, jets and the unclustered components.
• Invariant dilepton mass (m ): We can construct the invariant dilepton mass variable for two opposite sign leptons by defining: 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: 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. 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 III where we have listed the production cross-sections for our chosen BPs (table II) both in the presence and in the absence of Z B−3Lτ . As expected, in each case, the production via Z and Z B−3Lτ 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  ZZ also vanishes because of imposition of the m cut and W W is also killed by putting a hard MET cut of 300 GeV. The only remaining background is due to W W Z but that also becomes insignificant due to large MET cut.
We finally plot the signal significance for the BPs in Fig. 14  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 < ∼ 0.3 and the gauge coupling g B−3Lτ < ∼ 0.25. However, the mass difference between the heavier and lighter physical states of the VLF: M ψ 2 − M ψ 1 = ∆M can still be large enough ∼ 500 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 θ is allowed from direct search, and hence large ∆M is possible to achieve.
We finally discussed possible signatures at colliders by analysing some of the benchmark points 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 [114] with a total decay width Γ = 4.07 × 10 −3 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: Γ Z inv = 499 ± 1.5 MeV [91]. 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 ψ 1 ψ 1 pairs as other invisible decay modes are kinematically forbidden because of large ∆M . These decay widths are given by: And the singly charged two particle states for the sub-matrix M SC are as follows : |G + h , |G + z 1 , |G + φ , |G + z 2 , |G + s , |G + z 3 ; Each of distinct eigenvalues of the amplitude matrix, M will be bounded from tree level unitarity as : |λ H | ≤ 4π, |λ S | ≤ 4π, |λ 1 | ≤ 8π, |λ 2 | ≤ 8π, |λ SΦ | ≤ 8π, where, x 1,2,3 are the cubic roots of the following polynomial equation: Appendix D: Relevant Feynmann Diagrams for DM (co-)annihilation 15. Annihilation (i = j) and Co-annihilation (i = j) type number changing processes for Vector like fermionic DM in the model. Here i, j, k = 1, 2; a, b, c = 1, 2, 3 and f stands for SM fermions.
FIG. 17. Feynmann diagrams for charged fermionic DM, ψ ± annihilation to SM particles in final states.