Vacuum Stability in Inert Higgs Doublet Model with Right-handed Neutrinos

We analyze the vacuum stability in the inert Higgs doublet extension of the Standard Model (SM), augmented by right-handed neutrinos (RHNs) to explain neutrino masses at tree level by the seesaw mechanism. We make a comparative study of the high- and low-scale seesaw scenarios and the effect of the Dirac neutrino Yukawa couplings on the stability of the Higgs potential. Bounds on the scalar quartic couplings and Dirac Yukawa couplings are obtained from vacuum stability and perturbativity considerations. The regions corresponding to stability, metastability and instability of the electroweak vacuum are identified. These theoretical constraints give a very predictive parameter space for the couplings and masses of the new scalars and RHNs which can be tested at the LHC and future colliders. The lightest non-SM neutral CP-even/odd scalar can be a good dark matter candidate and the corresponding collider signatures are also predicted for the model.


Introduction
The last missing piece of the Standard Model (SM) particle spectrum was found in 2012 with the discovery of a SM-like Higgs boson with a mass of about 125 GeV at the Large Hadron Collider (LHC) [1,2], followed by increasingly-precise measurements [3][4][5][6] on its spin, parity, and couplings to SM particles, all of which are consistent within the uncertainties with those expected in the SM [7]. On the other hand, there are ample experimental evidences, ranging from observed dark matter (DM) relic density and matter-antimatter asymmetry in the universe to nonzero neutrino masses, that necessitate an extension of the SM, often involving the scalar sector. Moreover, from the theoretical viewpoint, it is known that the SM by itself cannot ensure the absolute stability of the electroweak (EW) vacuum up to the The tree-level Higgs potential symmetric under the SM gauge group SU (2) L × U (1) Y is given by [110] V scalar = m 2 11 where the mass terms m 2 11 , m 2 22 and the quartic couplings λ 1,2,3,4 are all real, whereas m 2 12 and the λ 5,6,7 couplings are in general complex. To avoid the dangerous flavor changing neutral currents at tree-level and to make Φ 2 inert for getting a DM candidate, we impose an additional Z 2 symmetry under which Φ 2 is odd and Φ 1 is even. This removes the m 12 , λ 6 and λ 7 terms from the potential and Eq. (2.2) reduces to 3) The EW symmetry breaking is achieved by giving real vacuum expectation value (VEV) to the first Higgs doublet, i.e with v 246 GeV, whereas the second Higgs doublet, being Z 2 -odd, does not take part in symmetry breaking (hence the name 'inert 2HDM').
Using minimization conditions, we express the mass parameter m 11 in terms of other parameters as follows: whereas the physical scalar masses are given by Since Φ 2 is inert, there is no mixing between Φ 1 and Φ 2 and the gauge eigenstates are same as the mass eigenstates for the Higgs bosons. The Z 2 -symmetry prevents any such mixing through the Higgs portal. In this scenario, the second Higgs doublet does not couple to fermions. Moreover, we get two CP even neutral Higgs bosons h and H, where h is identified as the SM-like Higgs boson of mass 125 GeV discovered at the LHC. We also get one pseudoscalar Higgs boson A and a pair of charged Higgs bosons H ± . Notice from Eq. (2.6) that the heavy Higgs bosons H, A and H ± are nearly degenerate. Depending upon the sign of λ 5 one of scalars between H and A can be a cold DM candidate. Since all the physical Higgs bosons except h are Φ 2 -type, i.e., Z 2 -odd, this also restricts their decay modes.

The Fermion Sector
In the fermion sector, we just add SM gauge-singlet RHNs to the SM particle content to generate tree-level neutrino mass via seesaw mechanism. In the canonical type-I seesaw, we just add three RHNs N R i , where i = 1, 2, 3 and the relevant part of the Yukawa Lagrangian is given by where L ≡ (ν, ) L is the SM lepton doublet, Φ 1 = iσ 2 Φ 1 (with σ 2 being the second Pauli matrix), N c R ≡ N R C −1 (with C being the charge conjugation matrix), Y N is the 3×3 Yukawa matrix and M R is the 3×3 diagonal mass matrix for RHNs.
After EW symmetry breaking by the VEV of Φ 1 , the Y N couplings generate the Dirac mass terms for the neutrinos: which mix the left-and right-handed neutrinos. This leads to the full neutrino mass matrix After block diagonalization and in the seesaw limit ||M D || ||M R ||, we obtain the mass eigenvalues for the light neutrinos as whereas the RHN mass eigenstates have masses of order M R . From Eq. (2.10), it is clear that in order to have the correct order of magnitude of light neutrino mass m ν 0.1 eV, as required by oscillation data as well as cosmological constraints, the Yukawa couplings in the canonical seesaw have to be very small, unless the RHNs are super heavy. For instance, for M R ∼ O(100 GeV), we require Y N O(10 −6 ). We will see later that these coupling values are too small to have any impact in the RG evolution of other couplings, and thus, the RHNs in the canonical seesaw have effectively no contribution to the vacuum stability in this model.
However, most of the experimental tests of RHNs in the minimal seesaw rely upon larger Yukawa couplings [111,112]. There are various ways to achieve this theoretically, even for a O(100 GeV)-scale RHN mass. One possibility is to arrange special textures of M D and M R matrices and invoke cancellations among the different elements in Eq. (2.10) to obtain a light neutrino mass [113][114][115][116][117][118][119][120]. Another possibility is the so-called inverse seesaw mechanism [121,122], where one introduces another set of fermion singlets S i (with i = 1, 2, 3), along with the RHNs N R i . The corresponding Yukawa Lagrangian is given by where M R is a 3×3 Dirac mass matrix in the singlet sector and µ S is the small lepton number breaking mass term for the S-fields. In the basis of {ν c L , N R , S}, the full 9×9 neutrino mass matrix takes the form (2.12) After diagonalization of the mass matrix Eq. (2.12) we get the three light neutrino masses whereas the remaining six mass eigenstates are mostly sterile states with masses given by M R ± µ S /2. The key point here is that the presence of additional fermionic singlet and the extra mass term µ S give us the freedom to accommodate any M R values while having sizable Yukawa couplings.
Irrespective of the underlying model framework, if we take large Y N ∼ O (1), it will have a significant negative contribution to the running of quartic couplings via the RHN loop at scales µ > M R . This must be taken into account in the study of vacuum stability in low-scale seesaw scenarios, as we show below.

RG Evolution of the Scalar Quartic Couplings
To study the RG evolution of the couplings, the inert 2HDM+RHN scenario was implemented in SARAH 4.13.0 [124] and the β-functions for various gauge, quartic and Yukawa couplings in the model are evaluated up to two-loop level. The explicit expressions for the two-loop β-functions can be found in Appendix A, and are used in our numerical analysis of vacuum stability in the next section. To illustrate the effect of the Yukawa and additional scalar quartic couplings on the RG evolution of the SM Higgs quartic coupling λ 1 in the scalar potential (2.3), let us first look at the one-loop β-functions. At the one-loop level, the β-function for the SM Higgs quartic coupling in this model receives three different contributions: one from the SM gauge, Yukawa and quartic interactions, the second from the RHN Yukawa couplings and the third from the inert scalar sector: Here g 1 , g 2 are respectively the U (1) Y , SU (2) L gauge couplings, and Y u , Y d , Y e are respectively the up, down and electron-type Yukawa coupling matrices in the SM. We use the SM input values for these parameters at the EW scale [12]: λ 1 = 0.1264, g 1 = 0.3583, g 2 = 0.6478, y t = 0.9369 and other Yukawa couplings are neglected [11].
It is important to note that the RHN contribution to the RG evolution of λ 1 is applicable only above the threshold of M R . For illustration, we assume M R = 100 GeV and fix all other quartic coupling values to λ i = 0.1 (with i = 2, 3,4,5) at the EW scale. The added effects of these new contributions in Eq. (3.1) on the RG evolution of the SM Higgs quartic coupling λ 1 ≡ λ h as a function of the energy scale µ are shown in Figure 1. Here the red curve shows the RG evolution of λ h using β SM λ 1 only [cf. Eq. (3. 2)], while the blue curve shows the evolution using β SM λ 1 + β RHN λ 1 , and finally the green curve shows the full evolution using , that pushes the stability scale back to 10 8.5 GeV and makes λ h > 0 again near the Planck scale. As shown in the middle and right panels, for smaller Y N values, the RHN contribution to the running of λ h is negligible, and therefore, the red and blue curves almost coincide. In these cases, the addition of inert scalar contribution pushes the stability scale up to 10 10 GeV, and then λ h again becomes positive at ∼ 10 15 GeV.
For completeness, we show the full two-loop evolution using the β-functions given in Appendix A in Figure 2. In this case, the stability scale in the SM is 10 9.5 GeV, whereas including the inert scalar contribution always leads to a stable vacuum all the way up to the Planck scale, even for the case when the Yukawa coupling is chosen to be large, Y N = 0.4 (left panel). From this illustration, we conclude that although large Yukawa couplings involving RHNs in low-scale seesaw models tend to destabilize the vacuum at energy scales lower than that in the SM, the additional scalar contributions in the inert 2HDM extension under consideration here have the neutralizing effect of bringing back (or even enhancing) the stability up to higher scales, and in the particular example shown above, all the way up to the Planck scale.

Stability Bound
The variation of the stability scale with the size of Y N and λ i is depicted in Figure 3. For smaller values of λ i , say 0.1 (red curve), the stability can be ensured up to the Planck scale only for Y N ≤ 0.30, beyond which the negative contribution from the RHNs take over and pull λ h to negative values at scales below the Planck scale. As we increase the λ i values, the compensating effect from the scalar sector gets enhanced and stability can be ensured up to the Planck scale for higher values of Y N . This is illustrated by the blue curve corresponding to λ i = 0.2, for which Y N ≤ 0.50 is allowed. However, arbitrarily increasing λ i does not help, as the theory encounters

Perturbativity Bound
Apart from the stability constraints on the model parameter space, we also need to consider the perturbativity behavior of the dimensionless couplings as we increase the validity scale of the theory. We impose the condition that all dimensionless couplings of the model must remain perturbative for a given value of the energy scale µ, i.e. † † the couplings must satisfy the following constraints: where λ i with i = 1, 2, 3, 4, 5 are all scalar quartic couplings, g j with j = 1, 2 are EW gauge couplings, 3 and Y k with k = u, d, e, N are all Yukawa couplings. Figure 4 describes the variations of different dimensionless couplings with the energy scale µ. Here we have shown the two-loop RG evolution of g 1 (yellow), g 2 (dotted blue), λ h (green), λ 3 (red), λ 4 (purple) and λ 5 (blue) as a function of the energy scale µ for benchmark values of Y N = 0.53 and M R = 100 GeV and with the initial conditions g 1 =0.3583, g 2 =0.6478, λ h =0.1264, and λ i = 0.4 (for i = 3, 4, 5) at the EW scale. Three important features are to be noted from this plot: (i) the λ h coupling becomes non-perturbative at around the scale µ 10 15 GeV, driven by the large Y N value; (ii) the λ 3 coupling becomes non-perturbative around 10 8.5 GeV until about 10 16 GeV, again driven by the large Yukawa coupling; and (iii) the gauge couplings g 1 and g 2 hit Landau pole at around the scale µ 10 17.3 GeV. Together, these features imply that the model becomes non-perturbative below the Planck scale for the choice of parameters shown here, in particular for the large Yukawa coupling Y N chosen in Figure 4. Thus, perturbativity of the couplings up to the Planck scale is an additional constraint we have to take into account along with the vacuum stability constraint. The perturbativity behavior of the scalar quartic couplings λ 3,4,5 is studied in Figures 5-7 respectively. In each case, we consider three benchmark values for the Yukawa coupling Y N = 0.1 (left), 0.4 (middle) and 0.9 (right). In each subplot, the various curves correspond to different benchmark initial values for the remaining unknown quartic couplings at the EW scale: red, green, blue and purple respectively for very weak coupling (λ i = 0.01), weak coupling (λ i = 0.1), moderate coupling (λ i = 0.4) and strong coupling (λ i = 0.8), while the SM Higgs quartic coupling is fixed at λ h = 0.126 and one of the quartic coupling value is varied (as shown along the x-axis) at the EW scale. From Figure 5, we see that for a given Y N value, the scale at which λ 3 hits the perturbative limit decreases as the scalar effect is increased. For example, in the strong coupling limit (with λ 2,4,5 = 0.8 at the EW scale), λ 3 hits the Landau pole at µ ∼ 10 6 GeV making the theory non-perturbative much below the Planck scale. As we increase the Y N value (going from left to right panel), the perturbative limit is reached even for smaller values of λ i . For instance, for Y N = 0.9 (right panel of Figure 5), λ 3 hits the Landau pole even in the very weak coupling limit (with λ i = 0.01) at µ ∼ 10 12 GeV. The results for λ 4 (cf. Figure 6) and λ 5 (cf. Figure 7) are very similar to those of λ 3 discussed above. Figure 8 shows the bounds on Yukawa coupling Y N from perturbativity of λ i for different initial λ i values. Here the color coding refers to the size of the Yukawa coupling. For small Y N ∼ 10 −7 corresponding to the canonical type-I seesaw limit (sky-blue region), no significant effect of RHN is noticed on the perturbativity bound.  values, i.e. λ i ≤ 0.15 and Y N ≤ 0.3 for the given theory to remain perturbative till the Planck scale. For comparison, it is worth noting that the perturbativity limit on Y N derived here is a factor of few weaker than those coming from EW precision data, which vary between 0.02 to 0.07, depending on the lepton flavor, for the minimal seesaw case (i.e. without the inert doublet) [125][126][127][128][129].
approach by Coleman and Weinberg [130], and calculate the effective potential at one-loop for our model. The parameter space of the model is then scanned for the stability, metastability and instability of the potential by calculating the effective Higgs quartic coupling and demanding appropriate limits. We then translate it into constraints on the model parameter space. The tree-level potential of our inert 2HDM is given in Eq. (2.3). To ensure that the potential is bounded from below in all the directions the tree-level stability conditions are given by [110] Considering the running of couplings with the energy scale in the SM, we know that the Higgs quartic coupling λ h gets a negative contribution from top Yukawa coupling y t , which makes it negative around 10 10−11 GeV and we expect a second deeper minimum for the high field values. Since the other minimum exists at much higher scale than the EW minimum, we can safely consider the effective potential in the h-direction to be where λ eff (h, µ) is the effective quartic coupling which can be calculated from the RG-improved potential. The stability of the vacuum can then be guaranteed at a given scale µ by demanding that λ eff (h, µ) ≥ 0. We follow the same strategy as in the SM in order to calculate λ eff (h, µ) in our model, as described below.

Effective Potential
The one-loop RG-improved effective potential in our model can be written as where V 0 is the tree-level potential given by Eq. are the corresponding one-loop effective potential terms from the inert scalar doublet and the RHN loops in the model. In general, V 1 can be written as where the sum runs over all the particles that couple to the h-field, F = 1 for fermions in the loop and 0 for bosons, n i is the number of degrees of freedom of each particle, M 2 i are the tree-level field-dependent masses given by with the coefficients given in Table 1. In the last column, m 2 corresponds to the tree-level Higgs mass parameter. Note that the massless particles do not contribute to Eq. (4.5), and hence, neither to Eq. (4.4). Therefore, for the SM fermions, we only include the dominant contribution from top quarks, and neglect the other quarks. It is also important to note that the RHN contributions come after each threshold value of M R i . Using Eq. (4.4) for the one-loop potentials, the full effective potential in Eq. (4.3) can be written in terms of an effective quartic coupling as in Eq. (4.2). This effective coupling can be written as follows: Note that in the inverse seesaw case and in the limit µ S → 0, each of the RHN mass eigenvalue is double-degenerate, and therefore, we have an extra factor of two for each RHN contribution in Eq. (4.6). The nature of λ eff (h, µ) in our model thus guides us to identify the possible instability and metastability regions, as discussed below. We take the field value h = µ for the numerical analysis as at that scale the potential remains scale-invariant [131].

Stable, Metastable and Unstable Regions
The parameter space where λ eff ≥ 0 is termed as the stable region, since the EW vacuum is the global minimum in this region. For λ eff < 0, there exists a second minimum deeper than the EW vacuum. In this case, the EW vacuum could be either unstable or metastable, depending on the tunneling probability from the EW vacuum to the true vacuum. The parameter space with λ eff < 0, but with the tunneling lifetime longer than the age of the universe is termed as the metastable region. The expression for the tunneling probability to the deeper vacuum at zero temperature is given by where T 0 is the age of the universe and µ denotes the scale where the probability is maximized, i.e. ∂P ∂µ = 0. This gives us a relation between the λ values at different scales: where v 246 GeV is the EW VEV. Setting P = 1, T = 10 10 years and µ = v in Eq. (4.7), we find λ eff (v) =0.0623. The condition P < 1, for a universe about T = 10 10 years old is equivalent to the requirement that the tunneling lifetime from the EW vacuum to the deeper one is larger than T 0 and we obtain the following condition for metastability [8]: (4.9) The remaining parameter space with λ eff < 0, where the condition (4.9) is not satisfied is termed as the unstable region. As can be seen from Eq. (4.6), these regions depend on the energy scale µ, as well as the model parameters, including the RHN mass and the gauge, scalar quartic and Yukawa couplings (see also Ref. [132]). Figure 9 shows the variation of λ eff in our model with the energy scale for different values of λ i (with i = 2, 3, 4, 5) and M R values with a fixed Y N = 0.4. The three different lines correspond to different values of the top Yukawa coupling by varying the top mass from 170 to 176 GeV with median value at 173 GeV [10]. The red region in Figure 9 corresponds to the instability region and the yellow region below the horizontal line λ eff = 0 corresponds to the metastable region, whereas the green region above λ eff = 0 is the stability region. Figure 9(a) and Figure 9 the top mass). Figure 9(a), Figure 9(c) and Figure 9(e) [or Figure 9(b), Figure 9(d) and Figure 9(f)] show that for fixed λ i and Y N , the stability scale also gets enhanced as we increase RHN mass M R , because the RHNs contribute to the β-function only at scales µ ≥ M R . This is the reason for the discontinuity at M R value, which is In all the subplots, we have fixed λ 2 = λ 5 = 0.01. The red, yellow and green regions correspond to the unstable, metastable and stable regions, respectively. obvious in Figure 9(e) and Figure 9(f).
To see the individual effects of the scalar quartic couplings λ 2,3,4,5 on the stability scale, we show in Figure 10 the three-dimensional correlation plots for λ 3 versus λ 4 with energy scale µ for different values of Y N and M R with a fixed λ 2 = λ 5 = 0.01. As in Figure 9, the red, yellow and green regions correspond to the unstable, metastable and stable regions respectively. Figure 10 As can be seen from Figure 9, the stability scale crucially depends on the top Yukawa coupling. The running of λ eff also depends on the initial value of λ h , which comes from the experimental value of the SM Higgs mass. Figure 11  contours show the current experimental 1σ, 2σ, 3σ regions in the (M h , M t ) plane, while the dot represents the central value [12]. Figure 11(a) describes that for small Y N = 10 −7 , the current 3σ values for the Higgs boson mass and top mass mostly lie in the stable region. However, as Y N is increased to a large value of 0.38 in Figure 11(b), the Higgs boson mass value lies in the stable region but the top mass value lies in the unstable/metastable region. The bound that comes on Y N from stability for which both Higgs boson mass and the top mass lie in the stability region is Y N 0.32 for M R = 100 GeV and λ i = 0.1.

LHC Phenomenology
The collider phenomenology of inter Higgs doublet with RHN is quite interesting as some decay modes involving RHNs are not allowed due to the Z 2 symmetry and this feature can be used to distinguish it from other scenarios. The pseudoscalar boson, the heavy CP-even Higgs boson and the charged Higgs boson (A, H, H ± ) are all from the inert doublet Φ 2 , which is Z 2 odd and their mass splittings are mostly M W [cf. Eq. (2.6)]. However, mass splittings around > ∼ M W ± ,Z are also possible some parameter space. The Z 2 symmetry prohibits any kind of mass-mixing of these inert Higgs bosons with the SM-like Higgs boson, which is coming from Z 2 -even Φ 1 . The couplings of Φ 2 with fermions are also prohibited, leaving only the gauge and scalar couplings. Nevertheless, as shown above, the inert Higgs doublet Φ 2 plays a crucial role in determining the stability and perturbativity conditions, and therefore, it is important to study their potential signatures at colliders. In Table 2     are with BR(N i → HH ± ∓ ) ∼ 2.4 × 10 −4 % and BR(N i → AH ± ∓ ) ∼ 5.2 × 10 −5 % respectively, as given in Table 3 for Y N = 0.01 and M R = 1 TeV. As for the RHN production at the LHC, being SM gauge-singlets, they can only be produced via their mixing with active neutrinos in the minimal seesaw model. The dominant production modes are shown in Figure 13. There are two types of processes: (a)-(d) involve RHN production via off-shell Higgs boson from gluon-gluon fusion, whereas (e)-(f) involve production via off-shell W ± /Z from Drell-Yan processes. The next-to-leading order (NLO) cross-sections for Y N = 0.1, 0.4 and M R = 500 GeV, 1 TeV are given in Table 4 where other parameters are kept as in BP3 of Table 2. For the process N ν [cf. Figure 13(a)], the production cross-section at NLO for Y N = 0.1 and M R = 500 GeV is: σ(gg → i=1,2,3 N i ν i ) is ∼ 0.15 and 9.7 fb respectively at the LHC with 14 TeV and 100 TeV center of mass energy [133]. For pair production the cross-sections are 1.8 × 10 −4 and 1.2 × 10 −3 respectively at the LHC with 14 TeV and 100 TeV center of mass energy. Here we have used CalcHEP 3.7.5 [134] for calculating the tree-level cross sections and decay branching fraction and have chosen NNPDF 3.0 QED NLO [135] and √ŝ (parton-level center of mass energy) as the energy scale for the cross-section calculations. The third column of Table 4 also give NLO Drell-Yan cross-sections for the same scale and PDF. We can see that for √ s = 14 TeV at the LHC Drell-Yan processes are more dominant than gluon gluon fusion, whereas at √ s = 100 TeV gluon gluon fusion processes surpass Drell-Yan ones. Though the overall cross-sections are small, but higher luminosity LHC can probe these three-body decays. The maximum cross-section comes for Y N = 0.4 and M R = 500 GeV and for √ s = 100 TeV and these are 155.40 fb, 95.60 fb, 0.50 fb respectively for (gg Note that although such large values of Y N might have been excluded from indirect constraints such as EW precision data, it is still useful to get an independent direct constraint from the collider searches. Coming to the inert Higgs boson signatures we have to rely on the mass spectrum of the Higgs bosons which depend on the couplings λ 3,4,5 as shown in Eq. (2.6). Table 2 shows benchmark points with the λ 3,4,5 that are allowed by the vacuum stability and perturbativity conditions. Depending on the phase space available, the charged Higgs boson in this model can decay into AW ± and/or HW ± mostly via off-shell W boson as the heavy Higgs bosons stay degenerate. The lighter of A and H is the DM candidate and thus can give rise to the signature of monolepton plus missing energy or dijet plus missing energy. However, because of the Z 2 -odd nature of H, A, H ± we can only produce the charged Higgs bosons as pair or in association with H/A. The heavier of A/H in that case decays to dilepton plus missing energy via off-shell Z boson. The production of H ± pair gives rise to dilepton plus missing energy and H ± A/H give rise to trilepton or mono-lepton plus missing energy signatures, which can be searched for at the LHC and FCC-hh [136]. The inert Higgs boson productions in association with the DM candidate leaving to jet plus lepton and missing energy signatures are studied in Ref. [137]. The inert doublet signatures along with the three-body decays of RHNs with Higgs boson in the final state can shed light on this model at the LHC with higher luminosity.
The LHC phenomenology discussed here is different from U (1) extensions where the RHNs can be pair-produced at the LHC via the U (1) gauge boson [138][139][140][141][142]. Phenomenological signatures of such RHN decays in the type-I seesaw in presence of extra scalars have been studied in the literature [143][144][145][146][147]. Similarly, in the case of type-III seesaw, the RHNs have charged partner and couple to W ± bosons [148]. The LHC phenomenology of such extensions with and without additional Higgs doublet has also been looked into [149][150][151]. The inverse-seesaw phenomenologies probing the RHNs at the LHC along with heavier Higgs bosons were also examined [152,153].

Conclusion
We have considered a simple extension of the SM with a Z 2 -odd inert Higgs doublet, supplemented by right-handed neutrinos with potentially large Dirac Yukawa couplings. The neutral part of the inert-Higgs doublet is a suitable DM candidate, while the RHNs are responsible for the correct light neutrino masses via seesaw mechanism. We have studied the effect of these new scalars and fermions on the stability of the EW vacuum by performing an RG analysis for the scalar quartic couplings.
We find that the additional scalars enhance the EW stability bound with respect to the SM case, as expected. Although the introduction of RHNs with relatively larger Yukawa couplings can be a spoiler for vacuum stability, the inert doublet comes to a rescue by contributing positively to the β-functions. On the other hand, the scalar quartic couplings cannot take arbitrarily large values at the EW scale due to perturbativity considerations at higher scales. In particular, we find upper bounds on the scalar quartic couplings λ i (with i = 2, 3, 4, 5) and the Dirac Yukawa couplings Y N , depending on the RHN mass scale M R , to satisfy both stability and perturbativity constraints.
We also analyzed the RG-improved effective potential to identify the regions of parameter space giving rise to stable, metastable and unstable vacua. For fixed values of λ i , increasing Y N enlarges the unstable vacuum region, whereas decreasing Y N and/or increasing the RHN mass scale M R enhances the stability prospects.
We also studied the phenomenological signatures of the heavy Higgs bosons along with RHNs at the LHC and future 100 TeV collider. Since the heavy Higgs bosons in this model come from the Z 2 -odd doublet, they are relatively non-interacting with the SM particles and are almost mass-degenerate, thus making their collider searches rather difficult. We have identified some new three-body decay modes of the RHNs to heavy Higgs bosons (assuming that the RHNs are heavier than the Higgs bosons) which can be used to distinguish this model from other vanilla RHN models.