Study of Electroweak Vacuum Stability from Extended Higgs Portal of Dark Matter and Neutrinos

We investigate the electroweak vacuum stability in an extended version of the Standard Model which incorporates two additional singlet scalar fields and three right handed neutrinos. One of these extra scalars plays the role of dark matter while the other scalar not only helps in making the electroweak vacuum stable but also opens up the low mass window of the scalar singlet dark matter (<500 GeV). We consider the effect of large neutrino Yukawa coupling on the running of Higgs quartic coupling. We have analyzed the constraints on the model and identify the range of parameter space which is consistent with neutrino mass, appropriate relic density and direct search limits from the latest XENON 1T preliminary result as well as in realizing the stability of the electroweak vacuum upto the Planck scale.


Introduction
The discovery of the Higgs boson [1][2][3] has been considered as the greatest triumph in present day particle physics. Although the experimental search is still on in order to investigate the Higgs boson's properties, several theoretical and phenomenological reasons are there to push us toward hunting for an enlarged Higgs sector compared to the one present in Standard model (SM). For example the Higgs quartic coupling λ H in SM becomes negative at a high energy scale (Λ SM I ∼ 10 10 GeV) leading to a possible instability of Higgs vacuum. The present measured values of Higgs mass∼ 125.09 GeV [4] and top mass ∼ 173.2 GeV [4], suggest that the electroweak (EW) vacuum can be metastable [5][6][7][8][9][10][11].
However the conclusion exclusively depends on the precise measurement of the top mass [12,13].
Also the metastability of the Universe is not a very robust situation in the context of cosmological inflation [14]. One of the possible solutions to this is to introduce new physics between EW scale and Λ SM I . In view of SM's incompetence in resolving some of the other issues like dark matter, neutrino mass, matter antimatter symmetry, inflation etc, the introduction of new physics is of course a welcome feature. 1 pghoshiitg@gmail.com 2 abhijit.saha@iitg.ernet.in 3 asil@iitg.ernet.in In particular, SM fails to accommodate a significant share in terms of its content called dark matter(DM). The most economical and popular scenario is the singlet scalar extension of the SM [15][16][17][18][19][20][21][22][23] having Higgs portal interaction. The stability of the dark matter is ensured by imposing a Z 2 symmetry on it. The relic abundance and corresponding direct detection cross section are solely determined by the DM (the scalar singlet) mass and its coupling with Higgs ( portal coupling). However present experiments, LUX [24] and XENON1T [25], strongly disfavor the model below m DM < 500 GeV except the resonance region. The bound on Higgs invisible decay width further constrain the model for m DM < 62.5 GeV [26]. Hence a large range of DM mass seems to be excluded within this simplest framework which otherwise would be an interesting region of search for several ongoing and future direct [24,25,27] and indirect experiments [28]. It is interesting to note that the presence of extra scalar in the form of DM can shift the instability scale (Λ I ) toward larger value compared to the SM one (Λ I > Λ SM I ) [29][30][31][32][33][34][35]. On the other hand to accommodate non-zero neutrino mass via type-I seesaw mechanism, one can extend SM with three right handed(RH) neutrinos. The RH neutrinos, being SM singlet, will have standard Yukawa like coupling involving Higgs and lepton doublets. The presence of the neutrino Yukawa coupling affects the running of the Higgs quartic coupling similar to the top Yukawa coupling.
With an endeavour to make the EW vacuum absolutely stable upto the planck scale M P , in a scenario that can accommodate both the DM and massive neutrinos with large Y ν (in type-I seesaw) and simultaneously to reopen the window for lighter scalar singlet DM mass (< 500 GeV), we incorporate two SM real singlet scalars and three SM singlet RH neutrinos in this work. Similar models to address DM phenomenology involving additional scalars (without involving RH neutrinos) have been studied [35,[51][52][53][54][55], however with different motivations. Our set up also differs from them in terms of inclusion of light neutrino mass through type-I seesaw. The proposed model has several important ingredients which are mentioned below along with their importance.
• One of the additional SM singlet scalars is our DM candidate whose stability is achieved with an unbroken Z 2 symmetry.
• The other scalar would acquire a nonzero vacuum expectation value (vev). This field has two fold contributions in our analysis: (i) it affects the running of the SM Higgs quartic coupling and (ii) the dark matter phenomenology becomes more involved due to its mixing with the SM Higgs and the DM.
• The set up also contains three RH neutrinos in order to generate non-zero light neutrino mass through type-I seesaw mechanism. Therefore, along with the contributions from the additional scalar fields, neutrino Yukawa coupling, Y ν , is also involved in studying the running of the Higgs quartic coupling.
We observe that the presence of the scalar 4 with non-zero vev affects the DM phenomenology in such a way that m DM less than 500 GeV becomes perfectly allowed mass range considering the recent XENON-1T result [25], which otherwise was excluded from the DM direct search results [56]. We also include XENON-nT [25] prediction to further constrain our model. On the other hand, we find that the SM Higgs quartic coupling may remain positive till M P (or upto some other scale higher than Λ SM I ) even in presence of large Y ν , thanks to the involvement of the scalar with non-zero vev. We therefore identify the relevant parameter space (in terms of stability, metastability and instability regions) of the model which can allow large Y ν (with different mass scales of RH neutrinos) and scalar DM below 500 GeV. Bounds from other related aspects, e.g. lepton flavor violating decays, neutrinoless double beta decay etc., are also considered. The set-up therefore demands rich phenomenology what we present in the following sections.
The paper is organized as follows. In section 2, we discuss the set-up of our model and in section 3, we include the constraints on our model parameters. Then in the subsequent sections 4 and 5, we discuss the DM phenomenology and vacuum stability respectively in the context of our model. In section 6, we discuss connection of the model with other observables. Finally we conclude in section 7.

The Model
As mentioned in the introduction, we aim to study how the EW vacuum can be made stable in a model that would successfully accommodate a scalar DM and neutrino mass. For this purpose, we extend the SM by introducing two SM singlet scalar fields, φ and χ, and three right-handed neutrinos, N i=1,2,3 . We have also imposed a discrete symmetry, Z 2 × Z 2 . The field φ is odd (even) under Z 2 (Z 2 ) and χ is even (odd) under Z 2 (Z 2 ) while all other fields are even under both Z 2 and Z 2 . There exists a non-zero vev associated with the χ field. The unbroken Z 2 ensures the stability of our dark matter candidate φ. On the other hand, the inclusion of Z 2 simplifies the scalar potential in the set-up 5 .
The RH neutrinos are included in order to incorporate the light neutrino mass through type-I seesaw mechanism.
The scalar potential involving φ, χ and the SM Higgs doublet (H) is given by where The relevant part of the Lagrangian responsible for neutrino mass is given by where l Li are the left-handed lepton doublets, M N is the Majorona mass matrix of the RH neutrinos.
This leads to the light neutrino mass, with v = 246 GeV as the vacuum expectation value of the SM Higgs. Minimization of the potential V leads to the following vevs of χ and H 0 (the neutral component of H), as given by ( So after χ gets the vev and electroweak symmetry is broken, the mixing between H 0 and χ will take place and new mass or physical eigenstates, H 1 and H 2 , will be formed. The two physical eigenstates are related with H 0 and χ by where the mixing angle θ is defined by Similarly the mass eigenvalues of these physical Higgses are found to be 6 Note that due to the absence of any Z 2 breaking term in the Lagrangian of our model, panic vacua [55,64,65] do not appear here. Using Eqs. (5,6,7), the couplings λ H , λ χ and λ χH can be expressed in terms of the masses of the physical eigenstates H 1 and H 2 , the vevs (v, v χ ) and the mixing angle θ as Among H 1 and H 2 , one of them would be the Higgs discovered at LHC. The other Higgs can be heavier or lighter than the SM Higgs. Below we proceed to discuss the constraints to be imposed on the couplings and mass parameters of the model before studying the DM phenomenology and vacuum stability in the subsequent sections.

Constraints
Here we put together the constraints (both theoretical and experimental) that we will take into account to find the parameter space of our model.
• In order to keep the entire potential stable, one needs to maintain the following conditions involving the couplings present in V (considering all couplings as real) which followed from the co-positivity of the mass-squared matrix involving H, χ and φ [66,67].
• In addition, the perturbative unitarity associated with the S matrix corresponding to 2 → 2 scattering processes involving all two particles initial and final states [68,69] are considered. In the specific model under study, there are eleven neutral and four singly charged combinations of two-particle initial/final states. The details are provided in Appendix A. It turns out that the some of the scalar couplings of Eq.(1) are bounded by The other scalar couplings are restricted (in form of combinations among them) from the condition that the roots of a polynomial equation should be less than 16π (see Eq.(A.9) of Appendix A).
• To maintain the perturbativity of all the couplings, we impose the condition that the scalar couplings should remain below 4π while Yukawa couplings are less than √ 4π till M P . An upper bound on tan β(= v/v χ ) follows from the perturbativity of λ χ [70] with a specific choice of m H2 .
•  [70]. For our purpose, we consider sin θ 0.3 for the analysis.
Apart from these, we impose the constraints on Y ν from lepton flavor violating decays. Also phenomenological limits obtained on the scalar couplings involved in order to satisfy the relic density (0.1175 ≤ Ωh 2 ≤ 0.1219) [79] and direct search limits [25] by our dark matter candidate φ are considered when stability of the EW minimum is investigated.

Dark matter phenomenology
The scalar field φ playing the role of dark matter has a mass given by m 2 DM = (µ 2 φ + 1 2 λ φH v 2 ) as followed from Eq.(1). Before moving toward the relic density calculation in our model, we would like to comment on the simplest Z 2 odd scalar dark matter scenario in view of recent XENON 1T [25] result. Note that for the purpose of DM phenomenology in this case, the only relevant parameters are given by m DM and the Higgs portal coupling λ φH (or µ φ and λ φH ). In Fig.1 (left panel), we provide a contour plot for relic density consistent with the Planck result [79] in the λ φH − m DM plane indicated by the blue solid line. In the right panel of Fig.1, we provide the DM-nucleon cross section evaluated with the value of λ φH corresponding to the m DM value as obtained from the left panel plot. We then incorporate the direct search limits on the DM-nucleon cross section as obtained from LUX 2016 [24], and the most recent XENON 1T [25] result [25] in the same plot denoted by blue and red dashed lines respectively. We conclude that the dark matter mass below 500 GeV is excluded from the present XENON 1T [25] result. This result is indicated by the red portion of the contour line in the left panel, Following Eq.(13) we draw the Feynman diagrams for DM annihilation channels into SM particles and to the second Higgs in Fig.2.
It is expected that the DM candidate is in thermal equilibrium with the SM degrees of freedom in the early universe. We therefore proceed to evaluate their abundance through the standard freeze-out mechanism. The Boltzmann equation, is employed for this purpose, where n φ is the number density of the dark matter φ, H is the Hubble parameter, σv φφ represents the total annihilation cross-section as given by σv φφ = σv φφ→SM,SM + σv φφ→H1H2 + σv φφ→H2H2 . We consider here the RH neutrinos to be massive enough compared to the DM mass. So RH neutrinos do not participate in DM phenomenology. We have then used the MicrOmega package [80] to evaluate the final relic abundance of DM.
We first assume H 1 as the Higgs discovered at LHC i.e. m H1 ∼125.09 GeV [4] and the other Higgs is the heavier one (m H2 > m H1 ). It would be appealing in view of LHC accessibility to keep m H2 below 1 TeV. In this case limits on sin θ, tan β are applicable as discussed in Section 3 depending on specific value of m H2 [72]. Now in this regime (where m H2 is not too heavy, in particular m H2 < 1 TeV), sin θ is bounded by sin θ � 0.3 [72] and we have taken here a conservative choice by fixing sin θ = 0.2. Note that in the small sin θ approximation, H 1 is mostly dominated by the SM Higgs doublet H. In this limit the second term in Eq.(8) effectively provides the threshold correction to λ H [57, 81, 82] which helps in achieving vacuum stability as we will see later. Furthermore considering this threshold effect to be equal or less than the first term in Eq.(8) (i.e. approximately the SM value of λ H ), we obtain an upper bound on m H2 as m H2 < mH 1 tan θ . Therefore in case with m H2 > m H1 , our working regime of m H2 can be considered within mH 1 tan θ > m H2 > m H1 . We take m H2 to be 300 GeV for our analysis. Note that with small θ, λ χ almost coincides with the second term in Eq. (9). It is quite natural to keep the magnitude of a coupling below unity to maintain the perturbativity limit for all energy scales including its running. Hence with the demand λ χ < 1, one finds v χ > √ 3m H2 . To show it numerically, let us choose sin θ = 0.2, then we obtain 125 GeV < m H2 < 620 GeV. Therefore with 9 Figure 2: Diagrams contributing to φφ annihilation to SM particles and the other Higgs.
The parameters v χ is involved in the definition of tan β = v/v χ . Parameters (λ H , λ χ , λ χH ) can be written in terms of other parameters as shown in Eqs. (8,9,10). Among all the parameters in Eq.(15), λ φ does not play any significant role in DM analysis.
We first assume H 1 as the Higgs discovered at LHC i.e. m H1 ∼125.09 GeV [4] and the other Higgs is the heavier one (m H2 > m H1 ). It would be appealing in view of LHC accessibility to keep m H2 below 1 TeV. In this case limits on sin θ, tan β are applicable as discussed in Section 3 depending on specific value of m H2 [72]. Now in this regime (where m H2 is not too heavy, in particular m H2 < 1 TeV), sin θ is bounded by sin θ 0.3 [72] and we have taken here a conservative choice by fixing sin θ = 0.2. Note that in the small sin θ approximation, H 1 is mostly dominated by the SM Higgs doublet H. In this limit the second term in Eq. (8)  tan θ > m H2 > m H1 . We take m H2 to be 300 GeV for our analysis. Note that with small θ, λ χ almost coincides with the second term in Eq. (9). It is quite natural to keep the magnitude of a coupling below unity to maintain the perturbativity limit for all energy scales including its running. Hence with the demand λ χ < 1, one finds v χ > √ 3m H2 . To show it numerically, let us choose sin θ = 0.2, then we obtain 125 GeV < m H2 < 620 GeV. Therefore with m H2 = 300 GeV, a lower limit on v χ ≥ 520 GeV can be set. We consider v χ to be 800 GeV so that tan β turns out to be 0.307.
On the other hand, if we consider the other Higgs to be lighter than the one discovered at LHC, we identify m H2 to be the one found at LHC and hence m H1 ≤ 125 GeV. Then Eq.(4) suggests sin θ → 1 as the complete decoupling limit of the second Higgs. Following the analysis in [72,[83][84][85][86][87], we infer that most of the parameter space except for a very narrow region both in terms of mixing angle (sin θ ∼ 0.9) and mass of the lighter Higgs (m H1 ∼ 85 − 100) GeV, is excluded from LEP and LHC searches. Such a range is not suitable for our purpose as can bee seen from Eq. (8). In this In a direct detection experiment, the DM scatters with the nucleon through the exchange of H 1 and H 2 as shown schematically in Fig.3 . The resulting spin-independent cross-section of DM-nucleon elastic scattering is given by [35] : where µ n = mnmDM mn+mDM , f n = 0.284 [88,89]. The couplings appeared as λ 1 , λ 2 are specified in the list of vertices in Eq. (13). Below we discuss how we can estimate the relevant parameters (λ φ H ,λ χφ and m DM ) from relic density and direct search limits. For this purpose, we consider m H2 = 300 GeV and v χ = 800 GeV as reference values, unless otherwise mentioned.

DM mass in region R1: [150 GeV < m DM ≤ 500 GeV]
In this region any decay mode of H 1 and H 2 into DM is kinematically forbidden following our consideration for m H2 = 300 GeV. As stated before, we consider m H1 to be the SM like Higgs discovered at LHC, with v χ = 800 GeV and tan β is fixed at 0.307. Therefore in order to satisfy the relic density Ωh 2 = 0.1161 ± 0.0028 [79], we first scan over λ φH and λ χφ for different ranges of dark matter mass where sin θ is kept fixed at 0.2. The allowed range of parameter space contributing to the relic abundance satisfying the correct relic density is indicated on λ φH − λ χφ plane in Fig From the top left panel of Fig.4, the relic density contour plot (with a particular m DM ) in λ χφλ φH plane shows that there exists a range of λ φH for which the plot is (almost) insensitive to the change in λ χφ . This becomes more prominent for plots associated with higher dark matter mass. In particular, the contour line satisfying the correct relic density with m DM = 500 GeV depicts a sharp variation in λ χφ (below 0.4) with almost no variation of λ φH around 0.13. We now discuss the reason behind such a behaviour. We note that for λ φH > 0.13, the total annihilation cross section satisfying the relic density is mostly dominated by the φφ → SM, SM process, specifically φφ → W W, ZZ dominate. In our scenario, φφ → W W, ZZ processes are mediated by both the Higgses, H 1 and H 2 .
Although λ χφ is involved in the vertices characterizing these processes, it turns out that once both the H 1 , H 2 contributions are taken into account, the λ χφ dependence is effectively canceled leaving the φφ → W W, ZZ annihilation almost independent of λ χφ . Hence φφ → SM, SM depends mostly on λ φH . The other processes like φφ → H 1 H 2 (H 2 H 2 ) are subdominant (these are allowed provided m DM > 212.5(300) GeV) in this region with large λ φH . Then the total cross section σv φφ and hence the relic density contour line becomes insensitive to the change in λ χφ as long as it remains below 0.4 while λ φH > 0.13. This is evident in the top left panel of Fig.4. Similar effects are seen in case of lower m DM (< 500 GeV) as well.
Once we keep on decreasing λ φH below 0.13, it turns out that φφ → SM, SM becomes less important compared to the φφ → H 2 H 2 (in particular the t channel) with λ χφ beyond 0.4 (in case of m DM = 500 GeV). Note that the plot shows the insensitiveness related to λ φH in this low λ φH region for obvious reason. Similar results follow with m DM < 300 GeV also, where φφ → H 1 H 2 provides the dominant contribution in σv φφ . Based on our discussion so far we note that for λ χφ λ φH the channels with Higgses in the final states contribute more to total σv φφ . On the other hand for low values of λ χφ (although comparable to λ φH ), the model resembles the usual Higgs portal dark matter scenario where W bosons in the final state dominate. To summarize, • 150 GeV < m DM < 212.5 GeV: For low λ χφ , φφ → W + W − dominates. However for large λ χφ , φφ → H 1 H 1 becomes the main annihilation channel.
• 212.5 GeV < m DM < 300 GeV: New annihilation process φφ → H 1 H 2 opens up. This with φφ → H 1 H 1 contribute dominantly for large λ χφ . Otherwise the channels with SM particles in final states dominate.
• 300 GeV < m DM < 500 GeV: The annihilation channel φφ → H 2 H 2 opens up in addition to Then total annihilation cross section will be enhanced for m DM > 300 GeV case, i.e σv φφ = σv φφ→SM,SM + σv φφ→H1,H2 + σv φφ→H2,H2 becomes large compared to the 280 GeV < m DM < 300 GeV mass range where σv φφ→H2H2 is not present. This enhancement has to be nullified in order to realize the correct relic density and this is achieved by reducing λ χφ compared to its required value for a fixed λ φH and m DM in 280 GeV ≤ m DM < 300 region. Note that in view of our previous discussion, we already understand that φφ → H 2 H 2 becomes important compared to φφ → SM, SM process in the region with λ χφ λ φH . Hence the two mass regions (below and above 300 GeV) overlap in λ φH − λ χφ plane as seen in the top left panel of Fig.4 as well in Fig.5. The total annihilation cross section of DM depends on its mass also. However the small mass differences between the two overlapped regions have very mild effect on σv Tot . Similar effect should be observed below and above m DM ∼ (m H1 + m H2 )/2 = 212.5 GeV as φφ → H 1 H 2 opens up there. However we find that around the m DM = 212.5 GeV, even with λ χφ λ φH , the contribution from this particular channel to σv Tot is negligible as compared to φφ → SM SM contribution and hence we do not observe any such overlapped region there.
In the top right panel of Fig.4 we provide the spin-independent (SI) direct detection (DD) cross sections corresponding to the points in the left panel satisfying relic density data having different range of dark matter masses as indicated by the colored patches. We further put the LUX 2016 [24], XENON 1T [25] and nT (expected) lines on it. As known, for a lowerer cross section, it reaches the neutrino floor where signals from DM can not be distinguished from that of neutrino. We find that the scenario works with reasonable values of the parameters, i.e. not with any un-naturally small or large values of couplings. Note that once we use the XENON 1T [25] and projected XENON nT [27] limits on the scattering cross section, we would obtain more restricted region of parameter space for λ φH − λ χφ as shown in left (with XENON 1T [25]) and right (with XENON nT [27]) figures of the bottom panel.
From the plot with XENON-nT prediction, we find that the scenario works even with reasonably large values of λ φH , λ χφ required for satisfying the relic density, although they are comparable to each other. This is because of the fact that to keep the direct detection cross section relatively small (even smaller than the XENON nT), it requires a cancellation between λ φH and λ χφ as can be seen from Eq. (16) in conjugation with definition of λ 1 and λ 2 for a specific sin θ = 0.2 value. Such a cancellation is not that important for plots with LUX 2016 [24] or XENON 1T [25] results and hence showing a wider region of parameter space for λ χφ and λ φH . It can be concluded from upper panel of Fig.4 that the presence of additional singlet scalar field χ helps in reducing the magnitude of λ φH that was required (say λ 0 φH ) to produce correct relic density in minimal form of singlet scalar DM or in other words it dilutes the pressure on λ φH to produce correct relic density and to satisfy DD cross section simultaneously. For illustrative purpose, let us choose a dark matter mass with 500 GeV. From Fig.1, we found that in order to satisfy the relic density, we need to have a λ 0 φH ∼ 0.15 which can even be 0.02 in case with large λ χφ ∼ 0.6. Similarly we notice that for m DM = 300 GeV, λ 0 φH was 0.086 in order to produce correct relic density which however was excluded from direct search point of view. This conclusion changes in presence of λ χφ as we can see from Fig.4, (left panel) that m DM = 300 GeV can produce correct relic density and evade the direct search limit with smaller λ φH : 0.065 − 0.086. This is possible in presence of nonzero λ χH and small sin θ(∼ 0.2 here) which redistribute the previously obtained value of λ 0 φH into λ φH and λ χφ while simultaneously brings the direct search cross section less than the experimental limit due to its association with sin θ (see the definition of λ 1 and λ 1 ). for DM heavier than 150 GeV, we find m DM ∼ 300 GeV can correctly produce the relic density in the observed range and simultaneously evade the DD limit set by LUX 2016 [24]. This result is consistent with the plot in Fig.4. Similarly m DM ∼ 500 GeV is in the acceptable range, which is in line with observation in Fig.4. In the left panel of Fig.7 we also have another region of DM mass∼ 75 GeV having correct relic abundance however discarded by LUX 2016. The region was not incorporated in top left panel of Fig.4 as we have started with m DM bigger than 150 GeV only. The possibility of having dark matter lighter than 150 GeV in the present scenario will be discussed in the next 7 As expected, it would be always possible to satisfy the relic density and DD limits within this region. subsection. Since in obtaining the Fig.4, we have fixed sin θ, tan β and m H2 , below in Fig.6 and 8, we provide the expected range of two couplings λ χH and λ φH when sin θ, tan β are varied for dark matter mass m DM = 300 GeV . We find the variation is little sensitive with the change of both v χ and sin θ.
As v χ or sin θ increases for m DM = 300 GeV, it requires less λ χφ for a particular λ φH to satisfy the relic density. We have also applied the LUX 2016 [24] DD cross section limit in those plots and are indicated by solid black patches. In Fig.8, one dark blue dot has been put on the sin θ = 0.2 contour which will be used in study of Higgs vacuum stability as a reference point.

DM mass in region R2: (m DM < 150 GeV)
Here we briefly discuss the DM phenomenology in the low mass region m DM < m H 2 2 = 150 GeV.
In this region, the decay process of heavy higgs to DM (H 2 → φφ) will be active. For further low m DM < m H1 /2 62.5, both H 2 → φφ and H 1 → φφ decay modes will be present.
We perform a scan over the λ φH − λ χφ region to find the correct relic density satisfied parameter space with allowed direct detection cross section from LUX 2016 [24] and XENON 1T experiments [25].
The results are shown in Fig.9  as obtained in Fig.4. We also note that there exists a resonance region through H 1 near m DM ∼ 63 GeV, indicated by the blue patch. In this resonance region, the relic density becomes insensitive to the coupling and hence the blue patch is extended over the entire region of λ χφ , λ φH in the Fig.9.
Finally we attempt to estimate the sin θ required to provide the correct amount of modification over the minimal version of a real singlet DM having interaction with SM Higgs only in order to revive the 'below 500 GeV' DM into picture. In other words, the amount of sin θ should be enough to satisfy correct relic abundance and DD cross section limits of LUX 2016 [24] and XENON 1T [25] for this particular mass range. To do the analysis, we fix λ χφ ∼ 0.2 while three different values of λ φH at 0.04, 0.08 and 0.10 are considered for the study. We then provide the sin θ versus m DM plot in Fig.10 which is consistent with relic density and LUX 2016 limits. We infer that a sizable value of sin θ is required for this. With λ φH = 0.1, we have noted earlier from Fig. 1

Vacuum stability
In this section, we will discuss how the EW vacuum stability can be achieved in our model. For clarification purpose and a comparative study of it, we first discuss how the presence of different ingredients (three RH neutrinos, DM and extra scalar χ) can affect the running of the Higgs quartic coupling when added one after other. We first comment on the inclusion of the RH neutrinos and investigate the running of λ H . Then we study how the involvement of the scalar singlet DM field φ can alter the conclusion. Finally we discuss the result corresponding to our set-up, i.e. including the χ field as well.
In doing this analysis, the absolute stability of the Higgs vacuum is ensured by λ H (µ) > 0 for any energy scale µ where the EW minimum of the scalar potential is the global minimum. However there may exist another minimum which is deeper than the EW one. In that case we need to calculate the tunneling probability of the EW vacuum to the second minimum. The Universe will be in metastable state provided the decay time of EW vacuum is longer than the age of the universe. The tunneling probability is given by [5,6], where T U is the age of the universe. µ B is the scale at which probability is maximized, determined from β λ H (µ B ) = 0. Hence for metastable Universe requires [5] λ H (µ B ) > −0.065 where T U 10 14 yr is used. As noted in [6], for µ B > M P , one can safely consider λ H (µ B ) = λ H (M P ).
Before proceeding further, some discussion on the involvement of light neutrino mass in the context of vacuum stbaility is pertinent here. As stated before, the light neutrino mass is generated through type-I seesaw for which three RH neutrinos are included in the set up. We now describe the strategy that we adopt here in order to study their impact on RG evolution. For simplicity, the RH neutrino mass matrix M N is considered to be diagonal with degenerate entries, i.e. M i=1,2,3 = M R . As we will see, it is Tr[Y † ν Y ν ] which enters in the β function of the relevant couplings. In order to extract the information on Y ν , we employ the type-I mass formula 2M R . Naively one would expect that large Yukawas are possible only with very large RH neutrino masses. For example with M R ∼ 10 14 GeV, Y ν comes out to be 0.3 in order to obtain m ν 0.05 eV. Contrary to our naive expectation, it can be shown that even with smaller M R one can achieve large values of Tr[Y † ν Y ν ] once a special flavor structure of Y ν is considered [38]. Note that we aim to study the EW vacuum stability in presence of large value of Tr[Y † ν Y ν ]. For this purpose, we use the parametrization by [90] and write Y ν as where m d ν is the diagonal light neutrino mass matrix and U PMNS is the unitary matrix diagonalizing the neutrino mass matrix m ν such that m ν = U * PMNS m d ν U † PMNS . Here R represents a complex orthogonal matrix which can be written as R = Oexp(iA) with O as real orthogonal and A as real antisymmetric matrices respectively. Hence one gets Note that the real antisymmetric matrix A does not appear in the seesaw expression for m ν = Y T ν Yν v 2 2M R . Therefore with any suitable choice of A, it would actually be possible to have sizeable Yukawas even with light M R and hence this can affect the RG evolution of λ H significantly. As an example, let us consider magnitude of all the entries of A to be equal, say a with all diagonal entries as zero. Then can be as large as 1 with a = 8.1 [90,91]. Below we specify the details of Higgs vacuum stability in presence of RH neutrinos only.

Higgs vacuum stability with right-handed neutrinos
In presence of the RH neutrino Yukawa coupling Y ν , the renormalization group (RG) equation of SM couplings will be modified [92].
Below we present the one loop beta functions of Higgs quartic coupling λ H , top quark Yukawa coupling y t and neutrino Yukawa coupling Y ν , where β SM λ H and β SM yt represent the β functions of λ H and y t respectively in SM. The Y ν dependence is to be evaluated in accordance with the type-I seesaw expression,  (elements of A), it is found [38] that Tr[(Y † ν Y ν ) 2 ] Tr[Y † ν Y ν ] 2 and we will be using this approximated relation in obtaining the running of the couplings through Eqs. (21,22,23). Here we have used the best fit values of neutrino oscillation parameters for normal hierarchy [93,94]. We have also considered the mass of lightest neutrino to be zero.
Note that just like the top quark Yukawa coupling, the neutrino Dirac Yukawa is having a similar impact on the Higgs quartic coupling, in particular with large Y ν . Also the top quark Yukawa would have a contribution dependent on Y ν . This has been studied in several works [36,[38][39][40][41][42][43][44][45][46]. We summarize here the results with some benchmark values of RH neutrino masses. These will be useful for a comparative study with the results specific to our model. In Fig.11 Table 1 at an energy scale µ = m t . Here we consider m h = 125.09 GeV, m t = 173.2 GeV and α s = 0.1184. In Fig.12, we have shown a region plot for and m t with fixed M R at 10 8 GeV in terms of stability (λ H remains positive all the way upto M P ), metastability and instability of the EW vacuum of the SM. The top quark mass is varied between 168 GeV to 178 GeV. The region in which EW vacuum is stable is indicated by green and the metastable region is indicated by white patches. The instability region is denoted with pink shaded part. It can be noted that the result coincides with the one obtained in [41]. We aim to discuss the change obtained over this diagram in the context of our model.

Higgs vacuum stability from Higgs Portal DM and RH neutrinos
Here we discuss the vacuum stability scenario in presence of both the scalar DM (φ) and three RH neutrinos (N ). In that case, effective scalar potential becomes V I + V H only. Note that the DM phenomenology is essentially unaffected from the inclusion of the heavy RH neutrinos with the assumption M R m DM . On the other hand combining Eq. (21,22,23), we obtain the corresponding beta functions for the couplings as provided below; From the additional term β II λ H , we expect that the involvement of DM would affect the EW vacuum stability in a positive way (i.e. pushing the vacuum more toward the stability) as shown in [29][30][31][32][33][34] whereas we noted in the previous subsection that the Yukawa coupling (if sizable) has a negative impact on it.
The interplay between the neutrino Yukawa coupling and Higgs portal coupling with DM is shown in Fig. 13, left and right panels (top and bottom). For the purpose of comparison, we have kept the same set of choices of parameters as in Fig.11, (left and right panels ). For the top panels, we consider mass of the dark matter to be m DM = 300 GeV and for the bottom set, m DM = 920 GeV is taken. The choice of m DM could in turn fix the λ φH coupling from the relic density plot of Fig.1. For example with m DM = 300 GeV λ φH is 0.075 and for m DM = 920 GeV, λ φH is given by 0.286 value. It is evident that the presence of Higgs portal coupling only has a mild effect as compared to the impact created by the neutrino Yukawa coupling. Finally in Fig.14 we provide the region plot in Tr[Y † ν Y ν ] -m t plane where the stable and instable regions are indicated by green and pink patches. This plot while compared with Fig.12, indicates that there is no such noticeable improvement except the mild enhancement of the metastable region due to the involvement of singlet scalar (DM) with Higgs portal coupling. With an aim to accommodate both the massive neutrinos and a relatively light dark matter (< 500 GeV), we move on to the next section where the χ field is included.

EW vacuum stability in extended Higgs portal DM and RH neutrinos
Turning into the discussion on vacuum stability in our framework of extended Higgs portal having three RH neutrinos, DM and the χ fields, we first put together the relevant RG equations (for µ > mt =173.2GeV    m DM , m H2 ) as given by, We note that the couplings λ χφ , λ φH and λ χH which played important role in DM phenomenology, are involved in the running of couplings as well. From the discussion of the DM section, we have estimated these parameters in a range so as to satisfy the appropriate relic density and be within the direct search limits for a specific choice of other parameters at their reference values: m H2 = 300 GeV and v χ = 800 GeV, sin θ = 0.2 (henceforth we describe this set as A). In particular an estimate for λ χφ , λ φH are obtained from Fig.4 (for 150 GeV< m DM < 500 GeV) and from Fig. 9 (for m DM < 150 GeV) having different choices of m DM and sin θ. The parameter λ χH dependence is mostly realized through sin θ following Eq.(10), where m H2 , tan β are fixed from set A. This sin θ is the most crucial parameter which control both the DM phenomenology and the vacuum stability. We have already seen that it allows the scalar singlet DM to be viable for the low mass window by relaxing λ φH from its sole role in case of single scalar singlet DM. On the other hand, a non-zero sin θ provides a positive shift (it is effectively the threshold effect in the small θ limit as seen from Eq. (8)) to the Higgs quartic coupling and hence guides the λ H toward more stability. Hence sin θ would be a crucial parameter in this study. Note that the RH neutrinos being relatively heavy as compared to the DM, neutrino Yukawa coupling does not play much role in DM phenomenology.
Assuming the validity of this extended SM (with three RH neutrinos and two singlets, φ, χ) upto the Planck scale, we study the running of the Higgs quartic coupling λ H from EW scale to M P as shown in Fig.15. In obtaining the running, we have considered m H2 = 300 GeV, sin θ = 0.2 and m DM is considered to be 300 GeV. The values of λ χφ and λ φH are fixed at 0.135 and 0.06 respectively (this particular point is denoted by a blue dot, named X, on Fig.8 ). It turns out that any other set of λ χφ and λ φH other than this blue dot from Fig. (while m DM = 300 GeV is fixed) would not change our conclusion significantly as long as sin θ is considered at 0.2. In order to compare the effect of the extra M R is fixed at 10 8 GeV. Contrary to our previous finding in section (see Fig.11, 13 ), we clearly see here that with M R = 10 14 GeV and m t = 171 GeV, λ H remains positive upto M P even in presence of Hence EW vacuum turns out to be absolutely stable. Although there exists other values of M R and/or m t , for which EW vacuum still remains unstable, the scale at which λ H enters into the instable region is getting delayed with a noticeable change from earlier cases (Figs.11,   13). This becomes possible due to the introduction of the χ field having contribution mostly from the sin θ parameter. In order to show its impact on stability, in Fig.16  The value of λ χ is then followed from Eq.(9) and λ φ is chosen to be at 0.7. Values of Tr[Y † ν Y ν ] = 0.24 and m t = 173.2 GeV are chosen for this purpose from Fig.17 (here the benchmark values are denoted by a black dot Y ). We conclude that all the stability criteria are fulfilled within the framework. Lastly we comment that instead of picking up the point X from relic density contour with sin θ ∼ 0.2 in Fig.8 to study vacuum stability in our model, we could have chosen any other point from that curve. As the stability of Higgs vacuum primarily depends on the value of θ, our conclusion would not change much.
However choice of any point having large λ χφ could make it reaching Landau pole well before M P in its RG running through Eq. (30). To avoid that one can reduce the value of λ φ ∼ O(10 −2 ) or less (earlier it was 0.7) which has no direct connection or impact on DM phenomenology and vacuum stability analysis in the proposed set up. In Fig.20, we have shown the running of all parameters from M R to M P involved in perturbative unitarity bound for the benchmark point: m H2 = 300 GeV, tan β = 0.30, sin θ = 0.2, m DM = 300 GeV, λ φH =0.06, λ χφ = 0.135, M R = 10 8 GeV and Tr[Y † ν Y ν ] = 0.24 with m t = 173.2 GeV. The parameters never exceed the upper limits coming from the unitarity bound. We have also confirmed that any other benchmark points wherever mentioned in our analysis satisfy the perturbativity unitarity limit. where in each cases neutrino Yukawa coupling Y ν has sizeable contributions. For this purpose, we consider m t = 173.2 GeV and M R = 10 8 GeV. From Fig.12, for SM + RH neutrinos, we see that stability can not be achieved. The metastability scenario is still valid in this case upto Tr[Y † ν Y ν ] < 0.26. Next we add a singlet scalar DM candidate with nonzero Higgs portal coupling to SM with RH neutrinos. Fig.14 (left panel) shows, for m DM = 300 GeV, stability of EW vacuum still remains elusive. On the other hand the metastability bound on Tr[Y † ν Y ν ] increases slightly from previous limit to 0.28. So DM with mass 300 GeV has mild impact on study of vacuum stability. Finally we add the extra scalar singlet with non zero vev to the SM with RH neutrinos and scalar DM. We have fixed the heavier Higgs mass m H2 = 300 GeV and sin θ = 0.2. Now in the combined set up of SM, scalar DM, scalar with non zero vev and RH neutrinos, the situation changes drastically from previous case as seen in Fig.17

Connection with other observables
In this section, we first discuss in brief the constraints on the parameters of the model that may arise from lepton flavor violating (LFV) decays. The most stringent limit follows from µ → eγ decay process. The branching ratio of such decay process in our set-up is given by [95][96][97] where α e = e 2 4π is the fine sructure constant, i runs from 1 to 3, The current experimental limit on LFV branching ratio is [4] Br(µ → eγ) < 5.7 × 10 −13 .
Using this limit, we therefore obtain bounds on |(Y † ν Y ν ) eµ | corresponding to a fixed M R value which can be converted to constrain Tr[Y † ν Y ν ] in our set up. In obtaining limits on Tr[Y † ν Y ν ] (for fixed M R ), first note that Y † ν Y ν remains function of M R and parameter a only (see Eq. (19) with O = I), once the best fit values of neutrino mixing angles [93,94] are used to evaluate U PMNS . Hence LFV limit basically constrains the parameter a which in turn is used to obtain Tr[Y † ν Y ν ]. This limit is shown on Fig.21 by the brown solid line, the left side of which is the disallowed region by LFV.
In the same plane of Fig.21 we also include the region of the parameter space allowed by both stability and metastability criteria. The green shaded region denotes the absolute stability of Higgs vacuum while the white region satisfies the metastability condition. We also indicate the instability region by pink patch in the same figure under discussion. For this purpose we have used m t = 173.2 GeV and m DM = 300 GeV, sin θ = 0.2, λ φH = 0.06 and λ χφ = 0.135 (corresponding to the benchmark point indicated by X in Fig.8). The brown shaded region is disfavored by the LFV constraint. Hence from Fig.21 we infer that for low M R , LFV constraints turn out to be stronger one and for high M R values, Tr[Y † ν Y ν ] is mostly restricted by the stability issue. It turns out that the proposed scenario does not provide any significant contribution to neutrinoless double beta decay [98][99][100][101][102][103] even for relatively low RH neutrino mass (∼ 10 3 GeV). This is in line with the observation made in [40]. Before concluding the section, it is perhaps important to comment on the possibility of explaining the baryon asymmetry of the Universe (BAU). The involvement of RH neutrinos would make the leptogenesis natural candidate to explain BAU from the completion point of view. However with the exactly degenerate RH neutrinos (we consider this for simplicity though), it is not possible. Once a small mass-splitting ∆M R between two heavy RH neutrinos can be introduced (for example by radiative effect [104][105][106]), resonant leptogenesis mechanism [107][108][109] can be succesfully implemented [110]. Apart from this, provided one can extend our vacuum stability analysis in presence of non-degenerate RH neutrinos [44] with DM and χ field, usual thermal leptogenesis can also be employed to explain the BAU of the universe.

Conclusions
We have considered an extension of the SM by three RH neutrinos and two scalar singlets with an aim to study the EW vacuum stability in a framework that can incorporate a stable light DM within the reach of collider experiments and to explain the light neutrino mass. A Z 2 × Z 2 symmetry is imposed of which Z 2 is broken from the vev of one of the scalars. It is known that with a real scalar singlet DM model, present experimental limits by LUX 2016 and XENON 1T rule out DM mass below m DM = 500 GeV. Also its presence does not modify the fate of EW vacuum much and hence keep it metastable only. Although metastability is acceptable, it however leaves some unwanted questions if we include primordial inflation in the picture. So an absolute stability of the EW vacuum is more favourable. On the other hand, introduction of RH neutrinos would have large impact on the running of the Higgs quartic coupling due to the neutrino Yukawa interaction. Provided the neutrino Yukawa coupling is as large as O(1) or more, it can actually destabilize the EW vacuum. Hence we have tried here achieving the stability of the EW vacuum in presence of RH neutrinos and DM. We also plan to find the possibility of a light scalar DM below 500 GeV. For this purpose, we have introduced additional scalar field which gets a vev. The other scalar among the two introduced does not get a vev and thereby is a good candidate for being a dark matter. The presence of the singlet with non-zero vev helps achieving the vacuum stability through a threshold like correction to λ H . So in this particular scenario i.e. SM extended by DM, three RH neutrinos plus one extra scalar, we have studied the Higgs vacuum stability issue considering large Yukawa coupling and variation of m t within 2σ range of uncertainty. We have found the stability region in the Tr[Y † ν Y ν ] − m t plane has been significantly increased in presence of χ. Simultaneously mixing of this extra scalar with SM Higgs doublet ensures its involvement in the DM annihilations. This mixing is effectively controlled by the Higgs portal coupling of the scalar which also enters into the running of the Higgs quartic coupling.
Hence an interplay between the two conditions: one is to achieve the EW vacuum stability and the other is to find a viable DM below 500 GeV, can actually constrain the parameters involved to some extent. Since the set-up involves several new particles, finding their existence in future and ongoing experiments would be an interesting possibility to search for. Here we have assumed the physical Higgs other than the SM one is heavier. The other situation where the second physical Higgs is lighter than the Higgs discovered at 125 GeV. However this case is not of very interest in the present study as following from Eq.(8), it can be seen that the effective Higgs quartic coupling becomes less than the SM one in this case and this would not help making EW vacuum stable. Also the sin θ allowed region for m H2 < m H1 /2 is almost excluded from the decay of H 2 → H 1 H 1 . Hence we discard this possibility. One interesting extension of our work could be the study of a SM gauge extension where the involvement of gauges bosons can modify our result. We keep it for a future study.