Naturally Stable Right-Handed Neutrino Dark Matter

We point out that a class of non-supersymmetric models based on the gauge group $SU(3)_C \times SU(2)_L\times SU(2)_R\times U(1)_{Y_L}\times U(1)_{Y_R}$ possesses an automatic, exact $Z_{2 }$ symmetry under which the fermions in the $SU(2)_R\times U(1)_{Y_R}$ sector (called $R$-sector) are odd and those in the standard model sector (called $L$-sector) are even. This symmetry, which is different from the usual parity symmetry of the left-right symmetric models, persists in the lepton sector even after the gauge symmetry breaks down to $SU(3)_C \times U(1)_{\rm EM}$. This keeps the lightest right-handed neutrino naturally stable, thereby allowing it to play the role of dark matter (DM) in the Universe. There are several differences between the usual left-right models and the model presented here: (i) our model can have two versions, one which has no parity symmetry so that the couplings and masses in the $L$ and $R$ sectors are unrelated, and another which has parity symmetry so that couplings are related; (ii) the $R$-sector fermions are chosen much heavier than the $L$-sector ones in both scenarios; and finally (iii) both light and heavy neutrinos are Majorana fermions with the light neutrino masses arising from a pure type-II seesaw mechanism. We discuss the DM relic density, direct and indirect detection prospects and associated collider signatures of the model. Comparing with current collider and direct detection constraints, we find a lower bound on the DM mass of order of 1 TeV. We also point out a way to relax the DM unitarity bound in our model for much larger DM masses by an entropy dilution mechanism. An additional feature of the model is that the DM can be made very long lived, if desired, by allowing for weak breaking of the above $Z_{2}$ symmetry. Our model also predicts the existence of long-lived colored particles which could be searched for at the LHC.


Introduction
The existence of dark matter (DM) constituting about 26% of the energy budget of the Universe is by now well established from astrophysical and cosmological observations [1]. It is also well known that understanding it requires the existence of new particle(s) beyond the Standard Model (SM) with very specific properties [2]. For instance, since the DM must be either absolutely stable or very long lived (i.e. lifetime longer than the age of the Universe), it implies that there must be an exact symmetry (or a very weakly broken symmetry) in the extended model under which all the SM particles are even while the DM particle is odd. 1 This symmetry not only makes the DM stable but also allows it to annihilate only in pairs to create its observed relic density in the Universe. If evidence for the DM decaying emerges in future data, one can accomplish this using a soft breaking of make it stable; see e.g. Refs. [30][31][32][33][34][35][36][37][38][39][40] for such RHN DM models based on the SU (2) L × U (1) Y × U (1) B−L gauge group. (ii) A second possibility in seesaw models is to have one of the RHNs decouple from the seesaw formula and remain light (with keV mass) as a result of a discrete flavor symmetry, e.g. L e − L µ − L τ [41][42][43], A 4 [44,45] or Q 6 [46], and thus it becomes a warm DM [14] in the Universe. One could of course have a RHN in a type-I seesaw framework with keV mass as in the case of νMSM [10,11] and have it as a DM as noted above.
A key model building issue in all DM models is whether there is a symmetry that keeps the DM naturally stable, and if this symmetry needs to be imposed by hand (like R-parity in case of MSSM) or is automatic by the gauge structure and matter content of the model. For example, in the supersymmetric case, extension of the MSSM to include U (1) B−L gauge symmetry can make R-parity an automatic symmetry [47][48][49]. In this paper, we discuss a class of non-supersymmetric gauge models based on the gauge group SU (3) C × SU (2) L × SU (2) R × U (1) Y L × U (1) Y R where there appears an automatic Z 2 symmetry that keeps the lightest RHN stable, thereby making it a natural DM candidate. Similar models for RHN DM have been discussed before in terms of its subgroup SU (3) C × SU (2) L × SU (2) R × U (1) Y [50,51]. However in these models, the DM is unstable and its decay to SM leptons was in fact used to explain the anomalous positron excess at GeV scale [50] or the PeV neutrino excess at IceCube [51]. The Z 2 symmetry appears in these models only when some parameters allowed by the gauge symmetry are set to zero, whereas in the model presented here, the Z 2 symmetry is present for all values of parameters allowed by the extended gauge symmetry of the model. In that sense, the Z 2 symmetry that stabilizes the DM in our model is an automatic symmetry. Another difference in the model presented here from those of Refs. [50,51] is that we keep the DM mass in the multi-TeV range, which has distinct implications for cosmology, direct and indirect detection, as well as collider searches. In fact, we note that combining collider and direct detection bounds, we can put a lower bound on the DM mass of order of 1 TeV. There is also an upper limit from unitarity arguments, which ranges from multi-TeV to PeV in our model, depending on the dilution mechanism. Another important feature of this class of models is that they exclusively lead to a type-II seesaw [52][53][54][55] for neutrino masses, unlike in the minimal LR models, where both type-I and type-II seesaw mechanisms are inherently present and one needs to do some fine-tuning to switch off either of them.
This paper is organized as follows: in Section 2, we present the details of the model. In Section 3 we calculate the relic density of RHN DM. In Section 4, we discuss the direct detection constraints, and in Section 5 the indirect detection constraints. In Section 6, we point out some LHC signals of the model. In Section 7, we briefly discuss other implications of our model, summarize our results and conclude. In Appendix A, we give the analytic expressions for calculating the annihilation cross sections for the relic density.

Model
In this section, we present two versions of the model, one where parity is a good symmetry relating the Yukawa couplings of the light and heavy sectors, and another where representation field representation (1, 1, 2, 0, 1) parity symmetry is explicitly absent from the beginning. In our subsequent discussion of phenomenology, we will focus only on the model without parity.

Model without parity symmetry
The model is based on the gauge group SU In addition to the SM fermions which transform under the SU (2) L ×U (1) Y L symmetry, we have a heavy analog which transforms nontrivially under the symmetry SU (2) R × U (1) Y R . In particular, the light fermions consist of the SM doublets Q L , ψ L and the singlets u R , d R , e R , which are collected in the upper left column of Table 2.1, with the symmetry SU (2) L × U (1) Y L being clearly the standard electroweak (EW) gauge symmetry. We will call this the L (light or left-handed) sector of the model. In the R (right-handed) sector, we have the heavy analog of the SM fermions: the SU (2) R doublets Q R , Ψ R and the singlets U L , D L , E L , which are collected in the upper right column of Table 2.1. It is clear that both the SM and heavy fermionic sectors share the common SU (3) C symmetry, and the electric charge formula is given by Here we do not have parity symmetry for reasons that will be discussed later and consider it as a UV-complete model for RHN DM. The Higgs sector of the model consists of SU (2) L doublet χ L and triplet ∆ L as well as their right-handed "partners" χ R and ∆ R , which are displayed in the lower part of Table 2.1. The doublets break the SM SU (2) L and the new SU (2) R gauge symmetries, as in the conventional LR models with vector-like fermions [56][57][58][59][60][61][62][63], after they obtain the non-vanishing vacuum expectation values (VEV) with v EW = 174 GeV and v R in the TeV range or higher. The doublets also give mass to all the fermions (except the heavy and light neutrinos) through the Yukawa interactions: where χ L = iσ 2 χ * L and similarly for χ R (σ 2 being the second Pauli matrix), f L, R stand for the SM chiral fermions and F L, R the corresponding heavy partners. Exact parity symmetry would have implied that the Yukawa couplings y f = y f (and also f = f in Eq. (2.6) below). If the SU (2) R × U (1) Y R symmetry breaking is in the few TeV range, this would imply new quarks in the few GeV range and would be inconsistent with observations, given the current LHC bounds on vector-like quark masses in the 1-1.5 TeV range [64][65][66]. Since in our model there is no parity symmetry, the Yukawa couplings in the heavy R-sector are independent of those in the SM sector and by choosing the y couplings to be of order one, we can have all the heavy fermions close to or above the TeV scale if v R few TeV. It is worth noting that at this stage unlike the standard minimal LR seesaw model, the lower bound on M W R [67][68][69] from low-energy flavor changing effects such as K 0 − K 0 mixing does not apply to our model. The triplet scalars ∆ L, R are used to generate the masses of light and heavy neutrinos via type-II seesaw mechanism [52][53][54][55]. To get the VEVs for the triplets, we choose the Higgs potential of the form: In the above expression, we have omitted the quartic terms of the form (χ † χ) 2 etc. since they do not affect our results and shown only the terms relevant for heavy and light neutrino masses after spontaneous symmetry breaking: We choose the parameters of the model such that w L ∼ eV (corresponding to the parameter m L ∼ 10 −6 GeV) and w R ∼ v R ∼ 10 TeV. Given the Yukawa interactions the ∆ R term gives masses to the RH neutrinos of order of few to 10 TeV whereas w L gives masses to the left handed neutrinos via the type-II seesaw mechanism. In the LR scenario without parity at the TeV scale, the f and f couplings are independent parameters, i.e. f = f , thus the RHN sector and the DM phenomenology in this paper are completely independent of the light neutrino sector. There are enough free parameters in the fcouplings that all neutrino oscillation parameters can be fitted without modifying the DM phenomenology discussed here. The f and f couplings however could be made equal in parity symmetric models, with slight extension of the Higgs sector; see Section 2.2.
It is remarkable that at the level of dimension four interactions, i.e. the Yukawa and gauge interactions and ignoring non-perturbative effects, the model has a large global symmetry: U (1) B, L × U (1) B, R × Z 2 , L × Z 2 , R even after the symmetry breaking VEVs are turned on. Here U (1) B,L is defined as the baryon number of SM quark fields Q, u, d and U (1) B,R the baryon number of the heavy quark fields of the SU (2) R sector. In the leptonic sector, the residual symmetries are two discrete symmetries, i.e. Z 2 , L defined as (−1) n where n is the lepton number of the SM leptons and similarly Z 2 , R for the heavy leptons of the SU (2) R sector. One of the most important implications of these symmetries is that the lightest of the heavy baryons made out of QQQ, QQq, Qqq is absolutely stable and the lightest of the heavy leptons is also absolutely stable. We will see later in Section 2.5 how these heavy baryons can be depleted during the cosmological evolution by introduction of new terms in the Lagrangian. We further note that if we choose the lightest lepton of the heavy sector to be the lightest of the RHNs, it will remain absolutely stable and can play the role of cold DM of the Universe. In the rest of the paper, we study the phenomenological implications of this model for DM and collider signals.
For the phenomenological purpose of avoiding the heavy lightest baryons also becoming DM and affecting the evolution of the Universe, we add two more Higgs singlets Σ 1,2 , as shown in Table 2.1. Both these fields have non-vanishing VEVs and connect the heavy and light singlet quarks via the following terms: 4 Once the Σ fields acquire VEVs, new heavy-light quark mass mixing terms δ U = λ U Σ 1 and δ D = λ D Σ 2 appear in the Lagrangian which break the global symmetry down to U (1) B × Z 2 , L × Z 2 , R at the tree level. Note that the Z 2 symmetry responsible for the stability of RHN DM is a subgroup of this symmetry under which the L-sector leptons are even and the R-sector leptons are odd, and thus it is a symmetry of the theory that remains at the renormalizable level. 5 These singlet VEVs also break the U (1) Y L × U (1) Y R gauge symmetries down to U (1) Y below the Σ scale. This is the version of the model discussed in Ref. [50,51]. We note that Y quantum numbers are different from the SM hypercharge (which we denote as usual by Y ) and B − L of the LR models.

Parity symmetric version of the model
The minimal model discussed above cannot be made parity symmetric for RH scale in the 1-10 TeV range. The reason is that parity symmetry would imply the Yukawa couplings in the two sectors to be equal. As a result, the lightest heavy quark mass in the 1 GeV range makes the minimal parity conserving theory phenomenologically untenable. The lowest right handed scale that would make this minimal theory acceptable is v R ∼ 10 7 − 10 8 GeV. 4 These multiplets are similar to those introduced in Ref. [70]. However, unlike in Ref. [70], we do not have the leptophilic scalar ΣE, so that the Z 2 , L × Z 2 , R symmetry remains exact. 5 Note that if non-renormalizable terms such as ELeRΣ 2 1 Σ2/Λ 2 terms are included, this symmetry will be broken and will make the RHN DM unstable; however a proper choice of Λ will make the DM long lived [51]. We do not discuss this possibility here.
A simple extension of the minimal model that can make it a viable theory for O(10 TeV) RH scale, is to double the number of electroweak doublets χ L,R in both sectors. The resulting Yukawa couplings for this case can be written as Note that the Yukawa couplings in the L and R sectors are now equal due to parity symmetry, unlike in Eq. (2.3). We then arrange the soft breaking mass terms for the χ L,R so that χ L,1 0 and χ L,2 = v EW with y u,d,e;1 1 whereas y u,2 : y c,2 : y t,2 = m u : m c : m t and similarly for the down-type quark and charged lepton Yukawa couplings y f,2 . In the RH sector, if we choose the VEVs of χ R,a ∼ v R 1−10 TeV, then the R-sector quark and lepton masses come almost entirely from the y u,d,e;2 couplings, and as a result, the R-sector heavy quark and lepton masses are in the TeV range as required by current LHC limits [64][65][66]. To give masses to the light and heavy neutrinos, we introduce only one triplet ∆ L,R with Y L,R = 2 as in the parity broken model discussed in Section 2.1. Parity symmetry now makes the masses in the two sectors proportional and mixings equal. Similar models with exact parity relating the heavy and light sectors were considered in Refs. [56,71,72] for the strong CP problem.
In what follows, we only consider the minimal parity broken model of Section 2.1.

Effective theory at the
In discussions of the phenomenological consequences below, the presence of higher gauge symmetry SU (2) L × SU (2) R × U (1) Y L × U (1) Y R is not important, but rather the reduced symmetry SU (2) L × SU (2) R × U (1) Y that emerges after the Σ-like fields acquire VEVs. Though the Y numbers of SM fermions are the same as the SM hypercharges, U (1) Y is not the SM gauge group U (1) Y , with the latter a combination of the former and the U (1) subgroup of SU (2) R after symmetry breaking. At this level, the gauge interactions become essentially that of LR models with quark seesaw [56,61,62], whereas the Yukawa interactions become different due to lack of parity symmetry. In addition, there are heavy and light quark mixings induced by which arises from the Lagrangian in Eq. (2.7). The global symmetries discussed above remain valid at this level. We assume that the heavy gauge boson Z LR which emerges from the breaking of Thus in the neutral gauge boson sector we can consider only the mixing involving the SM Z boson and the conventional Z R boson in LR models (see Section 2.4 below), which is similar as in Refs. [56,61,62] except that in our model the triplet ∆ R also contributes to Z R and W R masses.

Mixing between the heavy and light sectors
In this subsection we illustrate explicitly how the RHN DM particle in the heavy sector interacts with the SM fields, which will pave the way for all the phenomenological discus-sions below on the DM annihilation in the early Universe, its direct and indirect detection, and collider searches. In absence of any Higgs and fermion mixings connecting the heavy and light sectors, the RHN DM could only interact with the SM fermions through the heavy Z R boson, which couples directly to the SM quarks and leptons through the U (1) Y interaction. Once the mixing terms are included, in the lowest order, the RHN DM could talk to the SM fields via the scalar mixing (h − ∆ 0 R ) and the neutral gauge boson mixing (Z − Z R ) at the tree level. Regarding the physical scalars, only two are directly relevant to the DM phenomenology, i.e. the SM Higgs h and the neutral CP-even ∆ 0 R from the triplet ∆ R . The SM Higgs h is assumed to be predominantly from the doublet χ L , as in the conventional LR models [56,61,62], while ∆ 0 R couples directly to the RHNs. The mixing of the two scalars could be induced from the quartic couplings of form Here for simplicity we have neglected the effects from other scalars. After spontaneous symmetry breaking, the h − ∆ 0 R mixing reads with w R the non-vanishing VEV of the RH triplet ∆ R as defined in Eq. (2.5). Note that the scalar mixing ζ S contributes to the SM Higgs mass square, at the order of which requires a larger quartic coupling than in SM and thus could possibly help to improve the stability of the EW vacuum, as compared to the SM [61].
Regarding the Z − Z R mixing, the neutral gauge boson mass matrix reads as follows, after symmetry breaking at the RH scale and the EW scale, in the basis of (W 3L , W 3R , B) with B the gauge boson for the U (1) Y symmetry: where g L , g R and g Y are the SU (2) L , SU (2) R and U (1) Y gauge couplings, respectively. The VEVs in the mass matrix (2.11) have the hierarchical structure w L v EW v R , w R , and in this limit, the matrix can be diagonalized by the unitarity rotation where θ W is the Weinberg angle and tan φ ≡ g Y /g R . Obviously A is the massless photon, and we are left with two massive states, i.e. the SM Z boson and the heavy Z R boson, with masses respectively given by 13) and the mixing angle at the leading order given by where the parameter To zeroth order in ξ Z , the Z couplings are the same as in the SM, as they should be.
There are small corrections of order ξ Z 1. It should be noted that when the gauge coupling g R approaches the theoretical lower limit g L tan θ W (which is independent of the symmetry breaking pattern [73]), the ξ Z parameter, and therefore the Z −Z R mixing, could be significantly enhanced.
In the minimal version of our model, we do not have the bi-multiplet scalars which transform non-trivially under both SU (2) L and SU (2) R , like the bi-doublet (1, 2, 2, 0) in the conventional LR models; thus the charged gauge bosons W and W R can not mix at the tree level. Only when the SM and heavy fermions talk to each other via the mixing terms in Eq. (2.9), can the W − W R mixing be generated at the 1-loop level. The largest contribution stems from the mixing of third generation quarks and their heavy partners, and the corresponding Feynman diagram is presented in Fig. 1. At the leading order the charged gauge boson mixing parameter reads with m t, b the masses of SM top and bottom quarks, and M T, B the masses of heavy top and bottom partner fermions, respectively. Since the δ U, D terms break the global symmetry U (1) B, L × U (1) B, R , they are expected to be small, and therefore, the induced W − W R mixing is also a small number: for instance, ζ W ∼ 10 −12 δ U δ D GeV −2 for TeV-scale heavy partners. The W L − W R mixing and the mixing between light and heavy quarks are crucial to deplete the heavy hadronic states (see Section 2.5 below). One point to note is that these mixings do not affect the stability of the lightest heavy lepton, as it is odd under the Z 2 symmetry, and can be the DM candidate.

Depleting the heavy hadronic states
This model has two classes of baryons and mesons: (i) light baryons (qqq) and mesons (qq) that are part of the SM and (ii) heavy baryons and mesons which arise due to the fact that QCD is shared by the heavy quarks. In this subsection, we focus on the heavy baryons and mesons, and discuss the constraints on their properties from cosmology. Let us first note that there will be three kinds of baryons (QQQ, QQq, Qqq) and two kinds of mesons (Qq,QQ) involving heavy quarks. In the absence of any mixing between the heavy and light quark sector, the lightest of all these five kinds of states will be stable and their abundance in the early Universe will be determined by their masses, as in the case of strongly interacting DM [74][75][76][77][78][79]. Rough estimates in Ref. [79] give that Ω Q /Ω B ∼ 10 M Q /m p where M Q stands for the mass of lightest baryonic or mesonic state involving the heavy quark Q and m p is the mass of proton. This means that stable dark baryons or mesons above the GeV scale already over-close the Universe. Our goal is however to have the RHN as the only DM. Furthermore, heavy colored particles with masses up to a few TeV are incompatible with bounds on anomalous nuclei [80][81][82]. Masses heavier than this also seem to be excluded by considerations of DM-cosmic ray interactions producing gamma rays [83] and DM capture and self-annihilation in Earth's core producing internal heat flow [84,85]. So we would like to provide a mechanism to deplete the heavy quark bound states. The simplest way to do that is to introduce the heavy-light quark mixings given in Eq. (2.7). In fact requiring that they are depleted by the time of QCD phase transition temperature T QCD ≈ 200 MeV imposes lower bounds on the magnitudes of δ U, D which we estimate below.
To calculate the lower limits on the δ U, D , we note that they generate heavy-light quark mixings which can be determined by analyzing the following mass matrix, in the basis of (q R , Q R ) to (q L , Q L ), Diagonalizing this mass matrix, one can find that the heavy-light quark mixing in the righthanded sector is given by β R δ Q /y v R and in the left-handed sector by β L yv L δ Q /y 2 v 2 R . Thus the RH sector mixing is expected to be much larger. This mixing will allow states with heavy quarks to decay to light states plus SM gauge bosons. Typically for the lightest meson in the heavy sector, the decay goes likeQq →qq + Z. This decay process is much like the decay of free heavy vector-like quarks. The decay width can be estimated to be where M Q is the mass of the heavy vector-like quark with M Q M Z . Equating the decay rate (2.18) to the Hubble rate 1.66 g 1/2 * T 2 /M Pl (where g * is the effective relativistic degrees of freedom at temperature T and M Pl is the Planck mass) and taking the decay temperature T ∼ 1 GeV, we get a lower limit on the heavy-light quark mixing parameter, δ Q 10 −6 GeV above which value the heavy mesons will decay before QCD phase transition and not survive as the cold DM of the Universe. For the baryons, the decay width depends on the number of heavy quarks via (δ Q /v R ) 2N Q . Thus, it follows that if we choose the value of δ Q 10 −6 GeV, all heavy hadrons will disappear from the Universe before QCD phase transition, thus preserving the success of the Big Bang Nucleosynthesis.

DM relic density
In order to determine the relic density of RHN DM, we note that in the early Universe, all the particles were in equilibrium with the light SM sector particles due to the common SU (3) C color and U (1) Y interactions. As the Universe cools, the particles of the heavy sector being heavier than the DM N , slowly annihilate away leaving the N 's in the primordial plasma. As the temperature falls below M N , the DM density goes down and freezes out for T M N /20, as in case of a generic cold DM candidate. The primary annihilation channels to the SM particles proceed via particles that connect the two sectors, i.e. the SM Z boson and the heavy Z R boson which mix at the tree level, as well as the h − ∆ 0 R Higgs portal. Both the Z and Higgs portals are suppressed by the small mixing angles ζ Z and ζ S connecting the light and heavy sector, which are respectively of order λv EW /v R and v 2 EW /v 2 R . On the other hand, though the couplings of Z R to N and the SM fermions are of order one, the Z R portal is however suppressed by the large Z R mass, except near the resonance 2M N M Z R . Since the Z R channel is very important for the RHN DM relic density calculation, we will first present the current limits on Z R mass in our model from direct collider searches.

Limits on Z R mass
A heavy Z boson could decay promptly into the SM quarks and leptons, and the √ s = 13 TeV LHC dilepton searches require that M Z 3 TeV [86] for a sequential Z model with SM-like couplings [87]. However, the exact limit in LR models depends on the specific value of g R /g L [88,89], and moreover, in our case, the corresponding limits will be slightly weaker than in conventional LR models if we have the additional decay mode Z R → N N kinematically open. In this subsection, we reinterpret the √ s = 13 TeV LHC limits on Z mass for our Z R scenario.
To do this, we use the heavy neutral boson Z SSM in sequential SM [87] as a benchmark scenario, and rescale the cross section times branching ratio to electrons and muons and the corresponding mass limit as reported in Ref. [86]. In our model, the couplings of Z R to the SM fermions, as well as to the DM, are essentially proportional to their quantum numbers Y = Y L + Y R , rescaled by the gauge mixing sin φ; cf. the rotation matrix in Eq. (2.12). For the gauge coupling g R = g L in our model, the Z R production rate is 5 times smaller than a Z SSM ; on the other hand, the branching ratio to e + e − and µ + µ − is 1/4, which is much larger than that of Z SSM , if the decay mode of Z R to DM pair is not kinematically allowed. Then in our model the cross section times branching ratio of Z R is 0.80 times that of a  Model predictions for the Z R production cross section times branching ratio to dileptons at the √ s = 13 TeV LHC for three different values of g R /g L . The solid (dashed) curves are without (with) including the Z R decay to N N for its total width. The gray curve shows the current 95% CL upper limit from ATLAS [86]. Table 2. The lower limit on M Z R and the corresponding RH VEV w R (assuming v R = w R ) from LHC dilepton constraints [86]. The values in brackets assume that the decay mode Z R → N N is also kinematically allowed and M N M Z R . Z SSM boson. Including the Z R → N N decay mode, this limit goes down slightly. When g R is different from g L , due to the gauge coupling dependence of Z R production rate (as well as the lepton branching ratio if the DM mode is open), the Z R mass limits will be weaker (stronger) if g R becomes smaller (larger) than g L . For the sake of comparison, we show in Fig. 2 our model predictions for the Z R production cross section times dilepton branching ratio as a function of the Z R mass for three representative values of g R /g L and compare it with the latest 95% CL from ATLAS [86] to derive the bounds on M Z R and w R , as listed in Table 2.

Dominant annihilation channels
For completeness, we will consider both the scalar and gauge boson mediated channels, i.e.
Note that in current model the SM Higgs h and ∆ 0 R could in principle mix with the other scalars from the singlets Σ 1, 2 , the doublet χ R and the triplet ∆ L , and these scalars also contribute to the annihilation of N . However, we can always choose the quartic couplings in the scalar potential such that the mixings of h and ∆ 0 R to these scalars are negligible. The thermally averaged DM annihilation cross section times velocity σv in various channels are collected in Appendix A, where we list explicitly the coefficients a and b in Taylor expansion σv = a + b v 2 + O(v 4 ). Combining all these channels, we find that the annihilation of DM in our model depends on the gauge coupling g R , quartic coupling λ, as well as the RH scale w R and the masses and widths of ∆ 0 R and Z R , in addition to DM mass M N . The key Yukawa coupling f is related to the DM mass via the relation M N = 2f w R , which implies that for a light DM with f = M N /2w R 1 the scalar portal is further suppressed by the small Yukawa coupling f .
Once the coefficients a and b are known for all the available channels, the relic density of the DM can be calculated using the general formula [90] where x F = M N /T F 20 (with T F being the freeze-out temperature), g * = 106.75 the relativistic degrees of freedom at T F , and a and b the annihilation coefficients summing up all the available channels. An example is given in Fig. 3, where we set explicitly the quartic coupling λ = 1 and the heavy scalar masses M ∆ 0 R = 2 TeV with width 30 GeV. Three different value of the gauge coupling g R are chosen: g R /g L = 0.6, 1 and 1.5. The corresponding Z R mass and the RH scale v R = w R are set to their current experimental constraints listed in Table 2 without the Z R → N N decay mode. The horizontal dashed line shows the observed relic density, as measured by Planck [1]. The various peaks in Fig. 3 are respectively (from left to right) due to the SM Z and Higgs bosons, ∆ 0 R and Z R . We find that a TeV scale RHN could accommodate the observed DM relic density, with the annihilation dominated by the heavy scalar and/or Z R bosons, depending largely on the quartic coupling λ, the gauge coupling g R and the R-sector VEVs v R , w R . When the RHN is light, i.e. M N O(100 GeV), it could easily overclose the Universe even at the SM h (Z) resonance, as the coupling of h (Z) to DM is heavily suppressed by the mixing angle ζ S (ζ Z ). One should note in Fig. 3 for a fixed Z R mass, the DM N can not be arbitrarily heavy, as M N and M Z R are both proportional to the RH scales and M N /M Z R ∼ f /g R . In the very high mass limit, the annihilation cross section scales as ∼ M −2 N , until it hits the unitarity bound at sufficiently high energy scale, as discussed in Section 3.3.

Going beyond the unitarity bound
For generic range of parameters of the model, as the DM mass is increased beyond O(100 TeV) or so, the well-known partial wave unitarity limit kicks in [91]. To see this explicitly in our model, we write down the thermal averaged annihilation cross section (cf. Appendix A) at leading order in v 2 EW /v 2 R : Thus the cross section is suppressed by the right- which results in a Breit-Wigner enhancement. Even in this case, we must ensure that the maximum value of the cross section at the resonance obeys the partial wave unitarity limit [92]. For the resonance R(= ∆ 0 R , Z R ) just above the threshold 2M N M R , Eq. (3.3) can be approximated by where Γ N N and Γ SM are the partial decay widths of R to the DM pair and SM particles respectively, and B N N and B SM are the corresponding branching ratios. From Eq. (3.4), we find that the annihilation rate decreases with increasing DM mass, which leads to the unitarity bound of ∼ 20 TeV in our model [93]. Our goal in this subsection is to point out that in our model, there is a way to relax this generic bound without resorting to fine tuning of parameters (as required, e.g. in using an s-channel resonance below the threshold when 2M N M ∆ 0 R [51]). The basic idea is that in our model, there are naturally occurring long lived colored fermions (the SU (2) R quarks) which are heavy e.g. the lightest heavy SU (2) R sector quark, or the next to lightest RH neutrino, N 2 . We will show that there is a range of parameters in the model where, these long lived particles (denoted by X = Q, N 2 ) can decay to relativistic lighter SM fermions below the freeze-out temperature of the RHN DM. In that decay process sufficient entropy release can occur leading to dilution of the DM density to the right level. 6 There is a range of parameters of the model, where the heavy particle X decays at temperature T X to relativistic species after N 1 relic density is frozen (roughly around T F ∼ M N 1 /20) i.e. T X < T F . It will then generate entropy which can dilute the relic density of N 1 . There is also some dilution each time heavy species annihilate and disappear from the cosmic soup. The X decay temperature is given by We calculate Γ X for each choice of X = Q, N 2 and verify that X decays after T F by adjusting the mixing of heavy quarks to light quarks δ U,D and the mass M X . The dilution factor is then calculated as follows: equating the energy density before and after the decay we get where Y X = n X /s is the number density of X normalized to the entropy density. The dilution factor is given by As an illustration, we show below typical parameter values for the case of dilution by N 2 . We choose M N 2 ∼ 10 PeV and heavier quarks coupled to W R to have masses such that N 2 decay to them is kinematically suppressed/forbidden. The decay of N 2 then takes place via the mixings δ U,D of the heavy quarks to the lighter ones. We find that if δ U δ D ∼ 10 −3 GeV 2 we can get T N 2 ∼ 10 GeV (i.e. above the QCD phase transition scale). Taking the annihilation cross section for N 2 N 2 → SM fermions we estimate Y N 2 M N 2 at the time of N 2 decay which gets converted to entropy at T N 2 leading to a dilution factor d 100 for α Z R 10 −2 .
For the heavy quarks, due to QCD couplings being much larger than Z R coupling to N 2 , we find that the dilution factor from heavy quark decay is not very strong (roughly of order 2-3). There is also some dilution of DM relic density due to increase in entropy at the QCD phase transition point. Overall, it is possible to get a total dilution factor as big as d ∼ 10 6 , which can relax the unitarity bound for M N up to a PeV or so [93]. experiments. As a Majorana DM candidate, the RHN has both spin-independent (SI) and spin-dependent (SD) interactions with nuclei, which are mediated by the scalars and gauge bosons respectively. The SI scattering cross section reads, at the zero-momentum transfer limit,  Fig. 4, where we set the RH VEVs v R = w R = 1, 3 and 10 TeV, the gauge coupling g R = g L and the quartic coupling λ = 1. The colorful bands are due to the sizable uncertainties of effective DM-nucleon couplings f p, n [cf. Eq. (4.2)]. The SI DMnucleon scattering cross section is severely constrained, and currently the most stringent limit comes from the LUX experiment [101]. (with comparable limits from PandaX-II [102] up to 1 TeV DM mass). The future projected limits from XENON1T (with two exposure values of 2 and 20 ton·year) [103] and LZ [104] are also shown in Fig. 4 for comparison with the model predictions. The gray curves indicate the parameter space with the right relic density Ω N h 2 = 0.12, whereas the dashed gray part is excluded by direct LHC searches of Z R in the dilepton channel (cf. Section 3.1). Note that for a given DM mass, it is possible to have multiple solutions for the correct relic density depending on the Z R and ∆ 0 R masses (as in Fig. 3, where there are two solutions on either side of each of the ∆ 0 R and Z R resonances). Although the current direct detection limits are not stringent enough to probe the allowed parameter space of the model, the upcoming ton-scale experiments will have the sensitivity to probe part of this region, with the RH scale w R up to 10 TeV, far beyond the capability of direct production at LHC.
It is instructive to recast the SI limit shown in Fig. 4 onto constraints on the RH scale w R in Eq. (4.1) as a function of the DM mass M N , which is presented in Fig. 5. The parameter space below the colored blue curve in Fig. 5 is excluded by LUX, which sets a lower limit on the RH scale w R for any given value of the DM mass M N . When combined with the relic density constraint, this gives a lower limit on the DM mass of order 1 TeV, irrespective of the collider constraint on Z R . Note that for the relic density curves 7 These values are in good agreement with those extracted from an effective field theory approach [99,100]. We show also the current upper limits from LUX [101], as well as the future reaches of XENON1T [103] and LZ [104]. The gray curves correspond to the parameter space producing the observed relic density Ω N h 2 = 0.12 (with the dashed part excluded by direct collider searches of Z R in the dilepton channel, while the solid part is consistent with all constraints).
in Fig. 5, the vertical peaks (from left to right) correspond to the SM Z and Higgs bosons and the ∆ R resonances, whereas the slanted peak corresponds to the Z R resonance, whose location depends on the RH scale w R . As expected from Eq. (4.1), the direct detection constraints become more stringent when the DM becomes heavier. It is interesting to note that the future direct detection experiments are sensitive to the multi-TeV scale DM in our model, even up to 10 TeV, thus complementing the direct LHC searches. In this plot we set explicitly the quartic coupling λ = 1. For a coupling λ = 1, the constraints in Fig. 5 should be rescaled by a factor of √ λ, as indicated in Eq. (4.1). As for the SD cross sections, the leading order SD scattering of RHN DM from nuclei is from the axial-vector coupling of SM Z boson and Z R boson to partons in the detector nuclei. The Z portal is suppressed by the Z − Z R mixing M 2 Z /M 2 Z R while the Z R channel by its mass M Z R . It turns out that the two contributions are eventually of the same order, given by where J N is the total angular momentum quantum number of the nucleus, which equals to 1/2 for free nucleons, and λ q depends on the nucleus and reduces to the light quark contributions ∆ p q (∆ n q ) for scattering off free proton (neutron) [105]:  [neutron] [neutron] Ω N h 2 = 0 . 1 2 g R /g L = 1.0 Figure 6. Predictions for the SD scattering cross sections off neutron (blue horizontal line) and proton (orange horizontal line) assuming the gauge coupling g R = g L and w R = 1 and 3 TeV. We show also the current limit from LUX [106] and future limit from LZ [104]. See text and the caption of Fig. 4 for more details.
where we neglect the uncertainties on the effective DM-nucleon couplings. The predictions of SD cross sections are presented in Fig. 6, where we set g R = g L and w R = 1, 3 TeV. When the collider constraints on M Z R are taken into consideration, the Z − Z R mixing is so small that even with the future LZ limits it is rather challenging to test the model in terms of SD scattering, unlike the case of SI scattering, where large quartic couplings can enhance the testability of the model. So far in this section we have used the gauge coupling g R = g L . A smaller or larger  g R is also phenomenologically viable in different classes of LR models.
As demonstrative examples, we show the corresponding SI and SD scattering cross section plots with g R /g L = 0.6 and 1.5 in Figs. 7 and 9, and the lower limits on w R in Fig. 8. For the SI scattering, which is from the scalar channel, only the relic density lines Ω N h 2 = 0.12 are changed and shifted. For a smaller g R , the collider constraints of M Z R become more stringent (cf. Table 2], and the allowed RH scale w R is shifted to higher values, as presented in the left panels of Fig. 7 and 8. On the contrary, with a large g R , w R could be lowered and more parameter space can be probed by the upcoming DM direct detection experiments, as shown by the right panels in Fig. 7 and 8. The SD DM-nucleon scattering depends directly on the gauge coupling g R [cf. Eq. (4.3], and the most significant effect of changing g R is on the predictions of σ SD . It is readily understood that for a smaller (larger) g R , the effective coupling of Z and Z R to the RHN DM N is larger (smaller), and thus the cross section σ SD becomes larger (smaller). However, in both the cases shown in Fig. 9, the cross section σ SD for the allowed range of model parameters lies below the future projected limits. [neutron] [neutron] Ω N h 2 = 0 . 1 2 g R /g L = 1.5 Figure 9. Same as in Fig. 6, but with g R /g L = 0.6 (left panel) and 1.5 (right panel).

Indirect constraints
The RHN DM pairs in our model annihilate in the Universe into electrically charged SM particles, e.g. + − , qq and W + W − , mediated by the scalar and neutral gauge boson channels [cf. Eq. (3.1)] which could lead to energetic gamma-rays and thus be constrained by current and future gamma-ray observations. However, due to suppression caused by either small mixing angles ζ S, Z or large masses M Z R , ∆ 0 R , the resultant gamma-ray signals are below the current Fermi-LAT sensitivity [107], if the RHN DM is to have the observed relic density. Interpreting the Fermi-LAT constraints onto the limits on the model parameter space, we have the orange curve in the (M N , w R ) plot in Fig. 10 where we have set the gauge coupling g R = g L . In this plot, we show the combined Fermi-LAT constraint taking into account all allowed SM final states relevant to our model. The area below (or inside part of) the curve in excluded.
On the other hand, DM annihilation injects energy into the thermal bath in the early Universe, thus potentially altering the recombination history and thus changing the temperature and polarization power spectra of the cosmic microwave background (CMB). The Planck observations [1] of CMB anisotropies could therefore constrain the DM annihilation rate, which is complementary to other indirect detections of DM, such as the gamma-rays. The calculation procedure is quite similar to that for the Fermi-LAT limits, with one significant difference here being that we need to rescale the Planck limits given in Ref.  by the red curve in Fig. 10, where the area below (or inside part of) the curve is excluded. We find that the Planck limit is comparable to the Fermi-LAT limit, and in some region of the parameter space, even more stringent, in particular for a heavy DM of interest with M N 1.5 TeV. However, both the Fermi-LAT and Planck limits still cannot constrain any parameter space with the correct relic density. The same conclusion holds even for a different gauge coupling, as illustrated in Fig. 11 for g R /g L = 0.6 and 1.5. Other existing indirect constraints, e.g. from IceCube [109,110], are weaker than the Planck limits, and therefore, not shown here.

Collider phenomenology
As the time-reversed process of annihilation in the Universe, the RHN DM can be pairproduced at high energy colliders, such as the LHC, through the scalar (h and ∆ 0 R ) and neutral gauge boson (Z and Z R ) channels. As the heavy neutrinos are doublets under the gauge symmetry SU (2) R , it can also be produced through the W (R) portal, which is accompanied by a (off-shell) heavy charged lepton E in the final state. The heavy lepton E decays fast into the DM N plus a charged W (R) boson, which is different from the smokinggun signal of W R bosons at colliders in conventional LR models. However, the W − W R mixing in our model is heavily suppressed by both the loop factor and the soft breaking terms δ U, D (see Section 2.4), 8 thus it is rather challenging to search for the W R boson in the present model, and we will not consider such W (R) mediated processes in the discussions of DM searches at colliders below.

DM searches
As in case of typical WIMP DM searches at colliders, emission of one or more SM particles from the initial state partons is the smoking-gun collider signal for the RHN DM in our model. Emission of a hard gluon jet from one of the initial partons leads to a monojet plus large missing E T signal at the LHC, which has been searched for in the √ s = 13 TeV data [112,113]. There are also DM searches at the LHC in the mono-Higgs [114], mono-W/Z [115-117] and mono-photon [118] channels. Among these, the monojet channel turns out to be the dominant one for our RHN DM production at hadron colliders.
For pair production of the RHN DM with a high transverse momentum jet, the DM N could interact with the SM fermions through the SM Higgs or Z boson portal. However, both these two channels are from the heavy-light mixing angles ζ S or ζ Z , and thus the dominant channel is through the heavy Z R boson. Taking into consideration of the current Z R mass limit discussed in Section 3.1, our parton-level simulations reveal that the monojet cross sections in our present model are about three orders of magnitude smaller than the current LHC constraints [112]. Thus it is almost hopeless to see the RHN DM directly at the LHC in the monojet channel. Even at future 100 TeV collider [119], due to the rather low signal to background ratio, it is rather challenging to directly observe the RHN DM with the designed luminosity. The future lepton colliders, e.g. ILC, FCC-ee or CEPC, although much cleaner, have too low colliding energy to set any limits on the RHN DM in our model. So we have to look for alternative collider signatures of our model, as discussed below.

Searches for long lived particles at the LHC
As is clear from the discussion in Section 2.5, there is a viable range of the heavy-light quark mixing parameters δ U,D in our model for which the heavy quarks and resulting heavy baryons are long lived. Just to give a feeling for the numbers, we recall the expression for lifetime of a generic heavy quark decay Q → qZ in terms of the mixing δ Q as given in Eq. (2.18). For M Q ∼ 1 TeV, v R ∼ 10 TeV and δ Q ∼ 10 −6 GeV, this roughly gives a lifetime τ Q ∼ 10 −7 sec or the decay length L 0 ≡ cτ ∼ 30 m. For higher δ Q , this distance goes down like δ Q /10 −6 GeV −2 and could give rise to displaced vertices at the LHC. We expect the number of such events to be ∼ N Q (1 − e −L/L 0 ) /b where b = | p|/M Q is the boost factor for the produced heavy quark, N Q is the number of heavy quarks produced, L is the distance of the detector from the production point and is the detector efficiency. The current searches for displaced vertices [120] are roughly sensitive to δ Q ∼ 10 −5 GeV. The dominant production for the heavy quarks in our model at hadron colliders is through the processes gg, qq → QQ, with the cross sections with √ s = 14 and 100 TeV presented in Fig. 12 [121]. Note that the single production of heavy quarks in our model is suppressed by the mixing parameter δ Q . For a 1 TeV heavy quark, we expect a signal number of N Q = 3.6 × 10 5 at 14 TeV with an integrated luminosity of 3000 fb −1 , and a much larger number of N Q = 1.2 × 10 9 at 100 TeV with a luminosity of 30 ab −1 .
For neutral long-lived particles such as N 2 in our model, one may consider a different set-up as has been recently proposed [122]. More details of this will be given elsewhere.

Discussion and Conclusion
Before concluding, we briefly comment on some other phenomenological implication of the model: • Even though there are RH neutrinos with coupling to RH gauge boson W R , there are no new dominant contributions to neutrinoless double beta decay and the only observable channels are due to the canonical light neutrino exchange.
• Since the neutrino mass matrix fixes the Yukawa couplings of the left-handed triplet Higgs field ∆ L in our model with type-II seesaw, there will be contributions to lepton flavor changing processes such as µ → eγ, µ → 3e etc coming from the f -couplings. The severe experimental constraints on these processes impose constraints on the mass of the members of this triplet to be correlated to the value of v L ≡ ∆ 0 L . For example, for f ∼ 10 −2 , M ∆ 3 TeV.
• If needed, our model can naturally accommodate very long lived DM by soft breaking of the Z 2, symmetry stabilizing the DM. Depending on the DM mass, this could have additional phenomenological implications, e.g. at IceCube for a PeV-scale decaying DM.
In summary, we have discussed the phenomenology and cosmology of a TeV scale DM in the context of an SU (2) L × SU (2) R × U (1) Y L × U (1) Y R model where there is an automatic Z 2, symmetry that guarantees the stability of the DM. As a result, it provides an interesting possibility where the lightest right-handed neutrino N plays the role of cold DM. The DM relic density and direct detection constraints imply a lower limit of order 1 TeV for the RHN DM. The model also predicts the existence of long-lived heavy quarks which can be searched for at the LHC.
The subscripts S, f and V stand for the SM scalar, fermion and vector final states (all the SM fermion states are assumed to be summed over with the appropriate color factor N C ), and the superscripts S and V for the scalar and vector mediators, with the second superscripts V and A in Eqs. (A.2), (A.6) and (A.7) denoting the "vector" and "axialvector" parts of the gauge couplings to SM fermions. Those channels missing for the coefficient a are all vanishing in our model. δ S, V is a symmetry factor, with the value of 1 for identical final states and 2 for different states. C S i XX and C V i XY are the Yukawa and gauge couplings of the mediators to the DM particle N and the SM final states. In the present model the couplings of SM h and Z bosons to all the SM final state are the same as in the SM at the leading order, while the Yukawa couplings C ∆ 0 R N N = 1 2 f and axial-vector current coupling C Z R N N = g R /2 cos φ. The Z R boson couples to the SM fermions via the couplings C Z R f f = − 1 4 g R sin φ tan φ(Y L,f SM + Y R,f SM ) with Y L, R SM the SM hypercharge for the left-and right-handed fermions. The couplings of heavy bosons ∆ 0 R to the SM particles and the h and Z couplings to the DM N are essentially from the scalar and vector mixings and thus rescaled respectively by the small mixing angles of ζ S and ζ Z in Eqs. (2.10) and (2.14). P X = (4M 2 N − M 2 X + iM X Γ X ) −1 is the standard propagator for the mediator X. The velocities of final states are, at the leading order,