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 the 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 an enlarged parameter space for the DM, accsessible at collider experiments like the Large Hadron Collider (LHC), we consider the additional gauge symmetry to be based on the quantum B − 3Lτ. The restriction to third generation of leptons is chosen in order to have weaker bounds from the LHC on the corresponding gauge boson. The triangle anomalies arising in this model can be cancelled by the inclusion of a right handed neutrino which also takes part in generating light neutrino masses through type I seesaw mechanism. The model thus offers a potential thermal DM candidate, interesting collider signatures and correct neutrino mass along with 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.

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 singlet fermion and neutral component of a vector like SU(2) L doublet fermion, popularly known as singlet-doublet fermion DM [11][12][13][14][15][16][17][18][19][20][21][22][23]. While such singletdoublet 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 [24], 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 JHEP10(2019)275 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 [25], PandaX-II [26,27] and XENON1T [28,29] 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τ gauge charge. While gauged B −L symmetric extension of the SM [30][31][32][33] 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 [34][35][36][37][38][39][40][41][42][43][44][45]. 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 [46][47][48][49][50][51][52]). 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. 1 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 τ 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 In particular, the constraints from the requirement of vacuum stability restricts the gauge coupling of TeV scale B − 3L τ gauge symmetry to g B−3Lτ 0.25 and SM Higgs coupling with VLF to Y 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 [64].

JHEP10(2019)275
This paper is organised as follows: in section 2 we have introduced the new particles in our model and their interaction Lagrangian. In section 3 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 4 where in subsection 4.1 we have illustrated the relic density allowed parameter space and in subsection 4.2 direct search is discussed. The high scale validity and perturbativity of the model is elaborated in section 5, 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 6. Finally in section 7 we have concluded and summarised our findings.

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 τ charges and N R 3 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 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 have two more singlet scalars: S and Φ both charged under U(1) B−3Lτ . The additional SM singlet scalar fields take part in spontaneous symmetry breaking (SSB) of U(1) B−3Lτ , 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τ . 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 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 JHEP10(2019)275 following non-trivial relations: where x i stands for the U(1) B−3Lτ charges for all the particles appearing in table 1. The other anomalies, namely SU(2 are trivially cancelled. The U(1) B−3Lτ 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τ 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 JHEP10(2019)275 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: The first term in eq. (2.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. (2.7), Y B−3Lτ corresponds to B − 3L τ charge of the respective fermion. The first two terms in eq. (2.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: With one SM Higgs doublet and two non-standard singlets the renormalisable scalar potential takes the form: (2.9) 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 τ charge for the second scalar singlet Φ (-3/2) as the charge for S is already determined from RHN Majorana mass term in eq. 2.6 as discussed earlier. At the same time, this term also restricts the choice for the B −3L τ charge of the VLFs (+3/4), as evident from the last term in eq. (2.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 τ sector scalars are expanded around their respective vacuum expectation value (VEV): The gauged U(1) B−3Lτ is spontaneously broken as the two scalar singlets acquire nonzero 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 3. 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: Eq. (2.14) is important in our analysis as depending on m Z B−3Lτ 2.4 GeV 3 we can constrain our parameter space. Finally, the Lagrangian for the Yukawa sector

JHEP10(2019)275
involving the SM Higgs reads as: (2.15) 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.

VLF mass eigenstates
The Dirac mass matrix in the basis {χ, ψ 0 }, containing the singlet and doublet VLF can be written as: M χ M ψ . 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. (2.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:

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.

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.

Stability
Stability of the scalar potential is mainly dictated by the quartic terms of the scalar potential, V (H, Φ, S) which is defined as:

JHEP10(2019)275
In order to ensure the bounded-from-below condition in any field direction, the quartic couplings of the potential (eq. (3.2)) 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. (3.2) 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.

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 [76][77][78][79]. 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:Ŝ = αS 4s 2 w andT = αT where α is the fine structure constant and s w ≡ sin θ w corresponds to sine of the Weinberg angle θ 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 [6]: ∆T = T xSM − T SM = 0.07 ± 0.12, where we consider contributions JHEP10(2019)275 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. (3.7). The contribution from VLF DM is followed as, In eq. SM W and Z bosons respectively. The expression for T -parameter corresponding to the contribution from the VLFs is given by eq. (3.8), where g 2 is the SU(2) L SM gauge coupling. The Π's appearing in the expression are given as in eq. (3.9) and these correspond to the gauge boson propagator correction due to the VLFs. Here 'div' is the usual expression that appears in dimensional regularisation: div = 1 + ln 4π − γ , with γ = 0.577 is the Euler-Mascheroni constant ( = 4 − d, d ≡ spacetime dimension in dimensional regularisation). Note that the divergence appearing in the last term in eq. (3.8) (due to the divergences in eq. (3.9)) 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. (2.21)). Once more we would like to remind that M ψ 2 ≈ M ψ ± under small mixing limit as apparent from eq. (2.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]: Note that all the divergences appearing in eq. (3.11) and (3.12) along with the renormalization scale µ EW , are cancelled on substitution in eq. (3.10).
We have constrained the two most important free parameters of our model, namely DM mass M ψ 1 and ∆M using these constraints. These are depicted in figure 1. On the left hand side (l.h.s. ) of figure 1 we have shown the allowed values of M ψ 1 and ∆M that obey the constraint from T -parameter given by eq. (3.6). What we see from the plot in the l.h.s. of figure 1, for small VLF mixing (sin θ ≤ 0.5) it is always possible to get large ∆M 500 GeV for DM mass upto 1 TeV within the permissible range of the T -parameter. For sin θ 0.7, however, the parameter space a is constrained as ∆M as large as 1 TeV can not be achieved for low DM mass. In the r.h.s. of figure 1 we have shown the allowed range of S-parameter in M ψ 1 -∆M plane. Here we see, for small sin θ ∼ 0.1 the whole parameter space is allowed (gray region), however for larger sin θ 0.3 one has to stick to lower DM mass. Thus, S parameter constraints the DM mass for large sin θ. But in any case it is always possible to have a large ∆M within the observed range of S and T parameter.

Constraint on Z B−3Lτ 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 /g ≥ 7 TeV [83,84]. 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τ mass satisfying LHC bounds.
Search for heavy neutral Higgs and Z B−3Lτ resonances have been performed at the LHC [35], with the assumption that the heavy resonances decay to τ + τ − final states. These searches rule out m Z B−3Lτ < 2.42 TeV at 95% CL for sequential SM and m Z B−3Lτ < 2.25 TeV at 95% CL for non-universal G (221) model. We choose m Z B−3Lτ 2.5 TeV for a conservative limit. This puts a bound on three parameters in our model, namely: {g B−3Lτ ,ṽ, v Φ }. This is shown in figure 2, where each contour corresponds to m Z B−3Lτ = 2.5 TeV and hence the region right to each of the contours is allowed from collider constraint. We have chosen three different VEVs v Φ : {1.0, 2.0, 3.0} TeV corresponding to red, green and blue contours respectively. As it is seen, larger v Φ allows JHEP10(2019)275 larger gauge boson mass, which is in accordance with eq. (2.14). One should note here, eq. (2.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 v 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.

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 (0.2 sin θ m 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 θ m 0.2. For lower values i.e. m h 2 < 250 GeV, the LHC and LEP direct search [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.

Neutrino mass
The light neutrino mass matrix can be generated via type I seesaw mechanism as obtained from the Yukawa interaction of SM leptons in eq. (2.15) and singlet neutral fermions in eq. (2.6). Now, the neutrino mixing angles and mass squared differences are precisely measured from neutrino oscillation experiments [6]. 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. (3.13) with the usual 3 × 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 µ 1 = y µ 1 ≡ y l , , then we can produce correct order of light neutrino mass for y ∼ 0.1 and M 1 TeV. Even if we take M 10 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 S ∼ O(TeV) such that Z B−3Lτ 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 τ 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 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 τ 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. (3.13), will give rise to a phenomenologically unacceptable light neutrino mass matrix. This once again justifies the choice of singlet scalars and their B − 3L τ charges made in our model.
Here we would also like to mention that the TeV scale RHNs have decay lifetime τ N ∼ 10 −13 sec considering SM neutrino with scalar final state: N → ν, h 1 for our choice of Yukawa couplings. This shows that the RHNs do not contribute to the DM relic abundance as τ N τ universe (∼ 10 17 sec). Also, since they decay very fast ( 1 sec) to SM final states, they do not perturb the standard Big Bang Nucleosynthesis (BBN) picture and JHEP10 (2019)275 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 [54]. Before ending this subsection, we note that, although neutrino mass and mixing do not constrain the mass of B−3L τ 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τ /g B−3Lτ > 4.8 TeV [91]. Since we are considering the stronger bounds from the LHC and LEP II in our analysis, such weaker bounds are trivially satisfied.

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 2.1). In this section we have explored in detail the parameter space appearing in eq. (3.1) allowed by observed relic abundance (Ωh 2 = 0.120 ± 0.001 [4]) and direct search limits (particularly from the XENON 1T experiment [29]).

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 eff [92,93] which is given by

JHEP10(2019)275
In above equation, g eff 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 T f , where T f is the freeze out temperature. Relic density of ψ 1 one can approximately expressed as [94]: 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 [95] for which the model files are generated from LanHEP [96].
To see the behaviour of DM (ψ 1 ) relic density, we have fixed the VEV v Φ = 3.0 TeV such that the mass of Z B−3Lτ is always above the collider bound (m Z B−3Lτ > 2.  ∆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 figure 3 we have shown two different sets of {ṽ, g B−3Lτ } : {6, 0.05}; {1, 0.3} in solid black and dashed black curves respectively. This gives two different resonances at DM mass M ψ 1 ∼ 2.7 2 TeV and M ψ 1 ∼ 3.0 2 TeV due to two different masses of Z B−3Lτ . The nature of the two curves is almost identical except for two different resonances at m Z B−3Lτ /2. 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: Throughout the scan we have ensured that m Z B−3Lτ 2.5 TeV by properly adjusting g B−3Lτ andṽ, keeping v Φ fixed at 3 TeV as mentioned earlier.
The allowed parameter space from relic density requirement set by Planck experiment is shown in figure 4 in M ψ 1 -∆M plane. In the top left corner of figure 4 we have shown this parameter space for different ranges of the VLF mixing sin θ shown in red (0.01 ≤ sin θ < 0.05), green (0.05 ≤ sin θ < 0.1), blue (0.1 ≤ sin θ < 0.3) and black (0.3 ≤ sin θ < 0.5) where 0.05 ≤ g B−3Lτ ≤ 0.3. The relic abundance criteria is satisfied by moderate to large sin θ, while small sin θ's are confined near SM Z and SM Higgs resonance and near ∆M ∼ 10 GeV for M ψ 1 100 GeV. In order to understand the pattern more clearly we have chosen a fixed sin θ = 0.2 and plotted the same parameter space in the top right corner of figure 4. For DM mass around M ψ 1 ∼ 20 GeV there are only a few annihilation channels open for the DM. Now, for small ∆M co-annihilation comes into picture, increasing the effective annihilation cross-section (in eq. (4.2)). This causes the initial under abundance for small ∆M . On further increasing ∆M , co-annihilation becomes sub-dominant. As a result the DM becomes over abundant since the effective cross-section (eq. (4.2)) diminishes. For M ψ 1 40 − 70 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 ∆M region due to enhanced coannihilation and high ∆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 figure 3 as well. As mentioned earlier, the Higgs portal Yukawa Y becomes large enough in such a case (being proportional to ∆M for a fixed sin θ) resulting a net increase in the annihilation cross-section. For DM mass ∼ 1 TeV there is a huge overabundant region in the parameter space. This is due to the 1/M 2 ψ 1 suppression in the annihilation cross-section due to heavy DM mass. Correct relic abundance is still possible to reach at a very large ∆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 θ and for sin θ O(1) we can satisfy correct abundance in this region even with moderate ∆M . As we go beyond DM mass of 1 TeV, Z B−3Lτ resonance shows up (as the minimum value of m Z B−3Lτ is 2.5 TeV). Now, since both v and g B−3Lτ are being varied, m Z B−3Lτ is not fixed (according to eq. (2.14)). As a result, the resonance region is not sharp but broad due to different m Z B−3Lτ . For all possible choices ofṽ and g B−3Lτ according to eq. (4.5), m Z B−3Lτ is being varied between ∼ 2.5 TeV to ∼ 12 TeV. Because of this, the resonance band lies between M ψ 1 { 2.5 2 − 12 2 } TeV for all possible ∆M . It is seen that regions corresponding to over-abundance, under-abundance and right relic overlap on each other in the Z B−3Lτ resonance region due to different values of m Z B−3Lτ and ∆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τ .
As  relic for same DM mass are visible for M ψ 1 100 GeV. As it is observed, for DM mass say, 100 GeV (black dashed curve) the relic abundance first rises with increase in ∆M . This happens due to the fact that the co-annihilation becomes sub-dominant due to increase in ∆M , which, in turn, reduces the effective annihilation cross-section. At some point right relic density is reached (∆M ∼ 10 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 ∆M range and then again the abundance starts going downhill as the Yukawa Y becomes large enough making the net annihilation cross-section larger.  figure 4 we have shown the relic density allowed parameter space for different choices of the new gauge coupling g B−3Lτ . As we noticed earlier in figure 3, there is no strong dependence of the relic abundance on g B−3Lτ . This is evident from this plot as different coloured points (corresponding to different g B−3Lτ ) are scattered within the allowed region JHEP10(2019)275 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.

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,97]. As mentioned in section 2, 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. 2.2 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,97]: 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 [98], 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 [97], δ ∼ 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. (4.8) 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 figure 5. As we can see, in order to forbid such inelastic scattering one can choose y χ ∼ 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 θ 0.5 which is safe even if we consider the conservative bound from figure 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, JHEP10(2019)275 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 figure 6. The spin-independent (SI) direct detection cross section per nucleon is given by [99]: 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: The effective couplings (with form factors [100]) in eq.  with present bound from direct search experiment is shown in the l.h.s. figure 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 [101] where it is extremely difficult or even impossible to distinguish DM signal from the SM neutrino background (light grey region). In the l.h.s. of figure 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 ψ 1 100 GeV larger sin θ's are also allowed as the direct search cross-section has a suppression from heavy scalars: σ SI ∼ µ 2 sin 2 θ m 4 h i . Small sin θ 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 figure 7. On the r.h.s. of figure 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 ∆M for moderate sin θ. But if ∆M becomes too large 500 GeV then one has to resort to small sin θ to tame down the Yukawa Y in order to satisfy both relic abundance and direct search limits. For larger DM mass large ∆M is still allowed near the non-standard scalar resonances ∼ 100 GeV and ∼ 150 GeV. Beyond ∼ 150 GeV large ∆M is achieved with larger sin θ, while the DM remains still allowed by direct search due to suppression from heavy scalars mentioned earlier. With DM mass 400 GeV the points move towards smaller ∆M in order to reach right relic exploiting co-annihilation as we have seen earlier in figure 3. Beyond 1 TeV DM mass, the direct detection bound becomes weak as the DM mass is large, while because of Z B−3Lτ resonance there is a huge parameter space that satisfy relic abundance. As a result, for M ψ 1 1 TeV, almost all of the parameter space is allowed from direct detection for all possible ∆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 6.

High scale stability and perturbativity
In this section we will 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 [102] to extract the β 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 β 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. (3.4) and eq. (3.5)) for various couplings can be assured by analysing their runnings using the β functions. In our analysis, some of the Yukawa like couplings (y χ , y α i , y 13 , 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 JHEP10(2019)275 and scalar mixing angles 0.2. These in turn fix the magnitude of the scalar couplings which are positive and stays below ∼ 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τ 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τ plane and see 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. Recall that the EW vacuum stability (stability of Higgs potential) is dictated by the condition λ 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 v d ) can be written as [103,104], present model as provided below [105,106].
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. (3.3) 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 Y could destabilise the EW vacuum. The initial value of λ H also gets a positive shift due to the presence of additional scalars in the set up as evident from eq. (2.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 ∼ 0.02. Note that we have also considered all the other scalar couplings positive and ∼ O(0.1) in our analysis. Hence considering small order of magnitude of Yukawa like couplings (y α i , y 13 , y 23 , y χ ), 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 λ H > 0). Now we further constrain the Y − g B−L parameter space which is allowed from perturbativity criteria (figure 8) along with correct DM related observables using vacuum stability conditions. Before that in figure 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.    figure 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 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τ 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 6, 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τ

JHEP10(2019)275
Benchmarkṽ  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−3Lτ is also quoted. (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 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 [107]. Thus all our benchmark points are safe from LEP bounds.

Collider phenomenology
The detailed study of collider signature for vector like fermions can be found in [97,108]. 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 [97]. It is to be noted that the charged component of SU(2) L doublet VLF can be produced at JHEP10(2019)275 • Single lepton, with two jets plus missing energy ( ± + jj + / 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 figure 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 onshell in contrast to models where this decay takes place off-shell via SM Z boson and photon [20,97]. 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 l.h.s. of figure 12 for two different choices of {g B−3Lτ ,ṽ}. One noteworthy feature of this plot is that the production cross-section is lower for the choice {g B−3Lτ ,ṽ} = {0.2, 1.8} (black dotted curve) than for {g B−3Lτ ,ṽ} = {0.3, 0.8} (black dashed curve). This is simply attributed to the propagator suppression due to larger mass of Z B−3Lτ in the former case (m Z B−3Lτ = 3.36 TeV) over JHEP10(2019)275 the latter (m Z B−3Lτ = 2.54 TeV). On the r.h.s. of figure 12 we have illustrated how the production cross-section changes for different choices of g B−3Lτ andṽ (and hence m Z B−3Lτ ) keeping v Φ fixed at 3 TeV when both Z and Z B−3Lτ mediations are present. g B−3Lτ andṽ are chosen in such a way that m Z B−3Lτ 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 ∆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 figure 10. Although Y can be tamed down by choosing a small sin θ as per eq. (2.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.

Object reconstruction and simulation details
As already mentioned, we implemented this model in LanHEP and the parton level events are generated in CalcHEP [109]. Those events are then passed through PYTHIA [110] for showering and hadronisation. All the SM backgrounds that can mimic our final state are generated in MADGRAPH [111] and the corresponding production cross-sections are multiplied with appropriate K-factor [111] 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) [112]. 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.
• Jets (j): all the partons within ∆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 ∆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 |η| < 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 ( / 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: 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 (figure 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.

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τ . 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 figure 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 ψ 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 ∆M the signal distributions JHEP10(2019)275 Benchmark  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, ∆M . Hence for larger ∆M the DM carries away most of the p T making the distribution much flatter, while for smaller ∆M the distribution peaks up at lower value as the produced DM particles are not boosted enough. As a consequence, in the l.h.s. of top left panel of figure 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 ∆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 figure 13 we have shown the MET distribution for BP4 (red) and BP5 (green). Here we also see the same consequence: with comparatively larger ∆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 figure 13.
From the distributions one can easily see, with a MET cut of / E T 200 GeV and a H T cut of 250 GeV one can get rid off the SM backgrounds keeping most of the signal events intact. This is also shown in   The solid red and dashed red lines correspond to 3σ and 5σ discovery limits respectively.

JHEP10(2019)275
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| m + − |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 tt background is completely killed by imposing zero jet veto. Amongst other backgrounds, 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 figure 14 by choosing the final state events with / E T > 300 GeV and H T > 250 GeV. BP2 and BP4 can reach 5σ discovery for an integrated luminosity ∼ 200 fb −1 as they have the advantage of large ∆M which helps them to distinguish from the SM background as explained earlier. Due to comparatively smaller ∆M BP1 can reach a 5σ discovery at a slightly higher luminosity ∼ 500 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 ∆M 80 GeV. In that case the decay of the charged VLFs happen via off-shell W ± and the decay width can be small enough for small sin θ giving rise to charged tracks of length ∼ O(cm). This has been explored in details in [20,97] and hence we refrain from discussing it again here.

Summary and conclusion
We have proposed a flavoured gauge extension of the singlet-doublet fermionic dark matter model by considering B − 3L τ 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 τ 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 ∼ O(TeV) scale by tuning the new Yukawa couplings (y , y τ 3 ) accordingly. The family non-universal nature of this B − 3L τ gauge symmetry also helps us to avoid strong bounds from the LHC searches. The relatively lighter Z 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 θ 0.5. A substantial portion of the parameter space however merges with the neutrino floor for smaller sin θ.
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

JHEP10(2019)275
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 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σ 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 ∆M m W .
Acknowledgments 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.

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 [113] 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 [6]. 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:  2 λ 1 + 2λ 1 Y 2 + 2λ 1 (y 2 α 1 + y 2 α 2 + y 2 τ 3 + 4λ 1 y 2 χ − 16Y 2 y 2 χ ), (B.8) 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π, where, x 1,2,3 are the cubic roots of the following polynomial equation:   Figure 17. Feynmann diagrams for charged fermionic DM, ψ ± annihilation to SM particles in final states. Here a = 1, 2, 3. Figure 18. Additional Feynmann diagrams for DM, ψ i due to presence of new gauged paricle Z B−3Lτ in the model: annihilation (i = j) and Co-annihilation (i = j). Here i, j, k = 1, 2; a = 1, 2, 3 and q stand for SM quarks. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.