Effects of heavy neutrinos on vacuum stability in two-Higgs-doublet model with GUT scale supersymmetry

We analyse the implications of right-handed neutrinos on the stability of the electroweak vacuum in two-Higgs-doublet models with supersymmetry at high scale. It is assumed that supersymmetry is broken at scale $M_S = 2 \times 10^{16}$ GeV and effective theory below $M_S$ is two-Higgs-doublet model of type II with three generations of singlet neutrinos which induce small masses for the standard model neutrinos through type I seesaw mechanism. We study the high and low scale versions of seesaw mechanism. In both these cases, we show that the presence of right-handed neutrinos significantly improves the stability of electroweak vacuum if their Yukawa couplings with the SM leptons are of ${\cal O}(1)$ or greater. However, this possibility is severely constrained by the measured mass and couplings of Higgs and limits on the mass of the charged Higgs from the flavour physics data. It is shown that the stable or metastable electroweak vacuum and experimentally viable low energy scalar spectrum require $\tan\beta \lsim 2.5$ and the magnitude of neutrino Yukawa couplings smaller than ${\cal O}(1)$. The results obtained in this case are qualitatively similar to those without right-handed neutrinos.


I. INTRODUCTION
An embedding of the Standard Model (SM) into supersymmetric grand unified theories (GUTs) leads to an elegant and predictive framework which resolves several technical problems obstructing the extendability of the SM gauge theory to very short length scales. Supersymmetry (SUSY), close to TeV scale, stabilizes the electroweak scale against the large radiative corrections that arise due to higher fundamental scales present in the theory. TeV scale SUSY provides precision gauge coupling unification at the scale M GUT ≈ 2 × 10 16 GeV, although this can also be achieved if SUSY is broken at very high scale but some of the super-partners remain as light as few hundred GeV [1,2]. SUSY as a local symmetry is an essential ingredient in the superstring theory which provides a potential framework for unification of all the fundamental forces [3]. However, it is generically expected that the SUSY breaking scale in such a theory would be very close to the string scale. A similar situation arises in supersymmetric GUTs constructed in five or six dimensional spacetime in which the mechanisms used for breaking the unified gauge symmetry often break also the supersymmetry [4,5] 1 . In these frameworks, the effective theory below the GUT scale is described by non-supersymmetric SM which in some cases augmented by other possible light states as remnants of an underlying ultraviolet complete theory.
The experimental data collected in the first and second runs of the Large Hadron Collider (LHC) has shown no evidence for TeV scale supersymmetry [6,7]. If this trend continues in the future runs of LHC and other experiments then it would imply that supersymmetry is not an underlying mechanism for stabilization of electroweak scale. In this case, the assumption of existence of a weak scale supersymmetry is no longer necessarily required. Following this, in this paper we assume that SUSY exists in an underlying theory but it is broken at the scale much above the electroweak scale. It is assumed that supersymmetry does not play any role in stabilizing the electroweak scale and the electroweak scale remains finely tuned or some new dynamics in the high energy theory take care of the gauge hierarchy problem. Although an absence of SUSY at the low energy seems to make it less interesting from the phenomenological point of view, however its existence at high scale still leads to nontrivial consequences on the low energy theory. For example, the scalar potential of an effective theory below the SUSY breaking scale arises from the D-term potential of an underlying ultraviolet supersymmetric theory. Therefore, the electroweak symmetry breaking, stability of electroweak vacuum and mass spectrum of scalars in the effective theory are constrained by the high scale SUSY. Such consequences are already studied for an effective theory being only the SM [8][9][10][11], SM with additional Higgs doublet [12,13], and SM with higgsinos and gauginos [14].
It is observed that the SM alone as an effective theory cannot be matched to its minimal supersymmetric version (MSSM) when SUSY breaking scale is higher than 10 11 GeV because of the vacuum stability constraints [8,9] 2 . The SM with an additional Higgs doublet close to 1 Note that in [4], breaking of supersymmetry at high scale was prevented by introducing brane localized D-terms. Supersymmetry is broken at the GUT scale in absence of such terms. 2 The quoted limit is obtained by considering the mean value of the measured top quark mass. The limit the electroweak scale, known as two-Higgs-doublet model (see [15] for a review), is another possible effective theory with SUSY at high scale. The pair of scalar doublets in the two-Higgs-doublet model (THDM) can be identified as the two Higgs doublets of the MSSM. Such a matching between MSSM and THDM at high SUSY breaking scale has already been considered in [12,13,16,17]. It is found in [16] that THDM can be consistently matched to MSSM with SUSY breaking scale as high as the reduced Planck scale. This improvement over the SM is due to the presence of an additional Higgs doublet which modifies the stability conditions allowing more freedom in the effective potential. The stability or metastability of the electroweak vacuum however puts stringent constraints on the allowed values of tan β parameter which is the ratio of the vacuum expectation value (VEV) of two Higgs doublets.
In this paper, we investigate the effects of the so-called right handed (RH) neutrinos on the stability of scalar potential in THDM with supersymmetry broken at GUT scale. The SM augmented with such singlet fermions provides natural explanation for non-vanishing and tiny neutrino masses through type I seesaw mechanism [18][19][20][21][22][23][24]. In the GUTs based on SO(10) gauge group the RH neutrinos reside, along with the SM fermions, in three copies of 16-dimensional irreducible representation of the gauge group. In these models, if SUSY and gauge symmetry are both broken at the GUT scale leaving a pair of MSSM Higgs doublets light then the effective theory below the GUT scale is described by THDM with three generations of RH neutrinos. For example, this possibility is naturally realized in the GUT models based on flux compactification [5,25,26]. Motivated by this, we assume that the MSSM is broken at the scale M S = 2 × 10 16 GeV leaving THDM augmented with three generations of RH neutrinos as an effective theory below M S . The RH neutrinos obtain their masses through lepton number violating interactions which in turn induce tiny masses for the SM neutrinos through type I seesaw mechanism. We consider two distinct possibilities in which the mass scale of RH neutrinos is either close to M S or electroweak scale.
We find that RH neutrinos have considerable effects on the stability of electroweak vacuum. In particular, if these neutrinos are strongly coupled with the SM leptons then they lead to significant improvements in the stability of the scalar potential in THDM. It is shown that the stable vacuum can be achieved for almost any value of tan β if the magnitude of Yukawa couplings of RH neutrinos are larger than that of the top quark Yukawa coupling. However, the observed Higgs mass, measured couplings of Higgs with the gauge bosons and limits on the charged Higgs mass from flavour physics data severely constrain on this scenario.
The paper is organized as follows. We discuss the THDM framework with type I seesaw mechanism in the next section. The procedure of renormalization group (RG) evolution and matching at different scales have been discussed in section III. Numerical analysis and their results are discussed in section IV. The conclusion is presented in section V. Technical details related to renormalization group equations, threshold corrections at the high scale, extraction of the gauge and Yukawa couplings at top quark mass scale and dependency of results on the choice of top quark mass are elaborated in the Appendices.
is very sensitive with respect to the choice made for the value of top quark mass [10,11].

II. THE FRAMEWORK
The scalar potential of the most general THDM can be parametrized as where H 1 and H 2 are two complex Higgs fields, each of them is a doublet under SU (2) L and carries hypercharge Y = 1. With an addition of three copies of fermion singlet to the SM fermion spectrum, the most general Yukawa Lagrangian of the model can be written as where i, j = 1, 2, 3 stand for three generations of fermions andH 1,2 = iσ 2 H * 1,2 . If THDM is assumed to be an effective theory, obtained from the MSSM after the SUSY is broken at the scale M S , then the scalar potential in Eq. (1) is matched to the MSSM Higgs potential at M S . The potential of the MSSM Higgs doublets H u and H d (with hypercharge Y = 1 and Y = −1, respectively) contains the D-terms of superpotential and the soft supersymmetry breaking terms. SettingH 1 = H d and H 2 = H u , tree level matching between the potentials leads to the following conditions at M S [13,27]: where g Y = 3/5 g 1 . Further, the terms involving H † d and H † u in the Yukawa Lagrangian are absent at the scale M S . This implies in Eq.
(2) at the scale M S . The conditions in Eqs. (4,5) imply that the theory at M S is essentially type II THDM [15] with additional boundary conditions on the scalar quartic couplings as given in Eq. (3). We assume that the supersymmetry breaking scale is very close to the GUT scale. If the type I seesaw mechanism is considered as an underlying mechanism to generate tiny masses for the SM neutrinos, it introduces new scales in the theory, namely the mass thresholds of RH neutrinos. The Majorana masses for the RH neutrinos can be written as We denote the physical masses of RH neutrinos by M R i with i = 1, 2, 3 and adopt a convention in which It is typically expected that M R i lie in between the electroweak and GUT scale. For M R 1 H u , the effective light neutrino mass matrix becomes The Y ν and M R cannot be fixed uniquely from the available experimental information of neutrino masses and mixing parameters. This lack of information is best parametrized by the Casas-Ibarra parametrization [28] in which the Dirac neutrino Yukawa coupling matrix in the diagonal basis of the charged lepton and RH neutrino mass matrices is expressed as Here U PMNS is the leptonic mixing matrix, and m ν i are the light neutrino masses. R is an unknown complex orthogonal matrix which parametrize the freedom in choice of Y ν allowed by the seesaw formula, Eq. (7). We however do not consider this general case but discuss two phenomenologically interesting limits as described in the following.

A. High scale seesaw (HSS)
In SO(10) based GUTs, all the quarks and leptons of a given generation are embedded in a single 16-dimensional irreducible spinorial representation of the gauge group. In most of the situations, such a unification leads to an approximate equality between Y ν and Y u . For example, the renormalizable versions of supersymmetric SO(10) models with one or more 10-plet Higgs in the Yukawa sector always imply Y ν = Y u at the GUT scale [29,30]. The exact equality between Y u and Y ν at the GUT scale is broken if the underlying model contains higher dimensional Higgs representations, such as 126 and/or 120, or if the corrections from higher order non-renormalizable operators are taken into consideration. In many of these cases an approximate relation, Y ν ≈ Y u , still holds (see [31] for example). There also exists a possibility in which the hierarchy among the couplings in Y ν is widely different from those in Y u . This situation is known to arise from the orbifolded SO(10) GUTs in five or six spacetime dimensions [26,32,33]. If fermions are kept in the bulk and Higgs is localized on the brane, then the effective Yukawa coupling matrix in four dimensional theory is given by, is typically a matrix with elements of order unity, and F f L , F f R are diagonal matrices with elements representing the values of profile factors at the given four dimensional fixed point. The later decides the inter-generational mass hierarchies in a given fermion sector. In general, if the GUT symmetry is broken by orbifolding then F Q L = F L L and F u R = F ν R . Therefore, the resulting Y ν and Y u can have very different hierarchical structure and/or relative strength of magnitude [32,33].
To accommodate these possibilities, we generically parametrize Y ν at M S as: where , ξ and N are real numbers which determine the relative strength of Dirac Yukawa couplings of neutrinos with respect those of up type quarks. Since the largest coupling in Y u is already of O(y t ), very large value of leads to non-perturbative Y ν . We therefore consider ∈ [0.1 − 10], i.e. at most an order of magnitude difference between Y u and Y ν . The parameters ξ and N determine the hierarchical structure of couplings in Y ν . For ξ = 1, one obtains the hierarchy in Y ν same as that in Y u . Different hierarchical structure for Y ν can be obtained using suitably chosen values of ξ and N . For the above values of , Eq. (7) leads to the masses of RH neutrinos in the range 10 7 -10 16 GeV. This case is therefore named as the high scale seesaw (HSS) case. Since some of the couplings in Y ν are of O(y t ) or large, one expects considerable running effects from the RH neutrinos even though their masses are close to M S .

B. Low scale seesaw (LSS)
The running effects due to RH neutrinos are enhanced if the seesaw scale is close to the electroweak scale and they are strongly coupled with the THDM. A usual way to accommodate low seesaw scale is to consider 1 in Eq. (9) which in turn decreases the masses of RH neutrinos by a factor of 2 , as it can be seen from Eq. (7). Small however makes RH neutrinos very weakly coupled with the THDM and their effects on the running of couplings become negligible. There exists an alternate approach in which the low seesaw scale can be realized with O(1) couplings in Y ν [34][35][36]. In this case, the smallness of the SM neutrino masses is attributed to the flavour structure of Y ν and M R instead of the scale of RH neutrino masses or strength of couplings in Y ν . Due to the matrix structure of the seesaw formula, it is possible to choose the form of Y ν and M R such that the Eq. (7) leads to vanishing M ν . In the diagonal basis of RH neutrinos, they can be written as The above structures can be obtained from a global U (1) symmetry [34] or from a class of discrete symmetries [36]. One obtains M ν = 0 in this case irrespective of the values of y i , M and M 3 . It can be seen that the third RH neutrino does not couple with the SM leptons. The first two generations of RH neutrinos are strongly coupled if y i s are chosen to be of order unity. Viable neutrino masses can be generated by introducing perturbations to the above structure. The strength of these perturbations are found to be very small [36] and therefore their contribution to RG effects are negligible. The form of Y ν and M R given in Eq. (10) therefore provides a good description of low scale seesaw (LSS) with strongly coupled RH neutrinos.

III. RG EVOLUTION OF THE COUPLINGS AND CONSTRAINTS
The framework under consideration involves many hierarchically separated scales. We perform renormalization group evolution of the couplings of effective field theory between different scales and match their values at the boundaries. The couplings are evolved using 2-loop RG equations and matching at the thresholds are performed including 1-loop threshold corrections. The 2-loop RG equations are computed using a publicly available package SARAH [37] and they are listed in Appendix A. In the following subsections, we describe matching conditions, theoretical and phenomenological constraints on the couplings at various scales.

A. Matching conditions at M S
We assume that supersymmetry is broken at the scale M S = 2 × 10 16 GeV and the theory below M S is an effective THDM with or without RH neutrinos at intermediate scales between M S and M t . As it is discussed earlier, a tree-level matching between MSSM and THDM leads to relations given in Eqs. (3,4). The one-loop threshold corrections to these matching conditions and to the Yukawa couplings, in the absence of RH neutrinos, are given in [13,27]. These corrections depend on the sparticle spectrum at M S and also on the values of trilinear couplings and µ parameter. For simplicity, we assume that where mq and ml represent degenerate squark and slepton masses respectively, µ is higgsino mass parameter, M i are gaugino mass parameters and A t,b,τ are the trilinear couplings of squarks and sleptons with relevant MSSM Higgs fields. The above assumption is realized in specific GUT based model [5]. With these assumptions, the threshold corrections induced by squarks and sleptons are suppressed by the degeneracy of their masses and also by vanishing trilinear couplings. It can be seen from the expressions given in [13], the one-loop threshold corrections to Yukawa couplings vanish entirely for the superpartner spectrum given in Eq. (11). The threshold corrections to the quartic couplings, induced through one-loop box and triangle diagrams, depend only on the Yukawa couplings and µ parameter in the limit of vanishing trilinear couplings. We also estimate one-loop threshold corrections to quartic couplings which arise from the Dirac Yukawa couplings of RH neutrinos with the SM leptons. The expressions of threshold corrections, used in our analysis, are listed in Appendix B. For the analysis presented in this paper, we have chosen µ = 0.1M S for definiteness.

B. Constraints at intermediate scales
While evolving gauge, Yukawa and quartic couplings from M S to M t , we adopt the following procedure. If the mass scale M R j of j th RH neutrino ν j R appears below M S , their running effects are taken into account by appropriate RG evolution. It is expected that for renormalization scale Q < M R j , the ν j R should be integrated out from the spectrum and it should not contribute in the running of couplings. We implement this decoupling of heavy neutrinos by switching off the Yukawa couplings of ν j R with the SM leptons at the scale M R j and below. In other words, the values of elements of jth column in the matrix Y ν is put to zero after the running scale Q crosses the scale M R j . This procedure is carried out sequentially for all RH neutrinos with masses between M S and M t .
We also consider various theoretical constraints on the couplings which should be satisfied at every scale. It is to be noted from Eq. (4) and RG equations that the couplings λ 5,6,7 vanish at all the scales. This happens because these couplings (as well as the Yukawa couplings inỸ f for f = u, d, e, ν) are protected by a softly broken Z 2 symmetry of an effective THDM theory and therefore if they are zero at one scale then they will not be generated by running 3 . A stability of scalar potential in Eq. (1) would require the following conditions to be satisfied by the remaining couplings [38] for M t ≤ Q ≤ M S . The above conditions are sufficient to provide absolute stability for the electroweak vacuum. One may also consider a phenomenologically allowed and a more conservative possibility in which the electroweak vacuum is not completely stable but it is metastable with lifetime greater than the age of universe ∼ 10 10 years. This replaces the last condition in Eq. (12) by a weaker condition [16] λ(Q) ≡ 4 λ 1 (Q)λ 2 (Q) λ 4 (Q) where The derivation of the above condition involves probability of tunnelling into the true vacuum which was estimated in case of single scalar field with φ 4 potential in [39] including the quantum effects. Following a similar approach, the metastability condition, Eq. (13), was derived in [16] after mapping the THDM scalar potential into single field potential using the first three conditions in Eq. (12) with convinient choice of gauge and field basis. More details about the derivation of Eq. (13) can be found in an Appendix of [16]. For most of the cases studied here, it is found that the first three conditions of Eq. (12) are always satisfied as a consequence of the boundary values set by supersymmetry at M S . Hence, the stability or metastability of the electroweak vacuum is solely decided by values of λ 4 and λ at the intermediate scales.

C. Matching conditions and constraints at M t
The RG equations determine the values of the couplings of effective THDM at the scale M t . At M t , the gauge and Yukawa couplings are matched with their experimentally measured values while the quartic couplings determine the Higgs potential which is subject to the constraints imposed by consistent electroweak symmetry breaking and measurement of Higgs properties.
The electroweak symmetry breaking is governed by the VEVs of two Higgs fields, at the minimum of scalar potential. In our notation, v 1 = v d and v 2 = v u and they define electroweak VEV v and a parameter tan β as the following The breaking of electroweak symmetry gives rise to five physical Higgs bosons in the spectrum. These are two charged and CP-even (H ± ), two neutral and CP-even (h and H), and a neutral and CP odd (A) scalars. At the minimum of potential, the parameters m 2 1 , m 2 2 and m 2 12 can be replaced by the following tree-level expressions [38]: where M A is the mass of pseudo-scalar Higgs in MS renormalization scheme. With these replacements, the scalar potential given in Eq. (1) is completely specified by M A , v, tan β and the quartic couplings. The mass of charged Higgs is given by The CP-even scalar states mix with each other and it is convenient to work in so-called Higgs basis in which only one of the combinations of H 1 and H 2 , namely h 1 ≡ cos βH 1 + sin βH 2 , acquires a non-trivial VEV. The combination orthogonal to h 1 is identified as h 2 such that h 2 = 0. In the basis {h 1 , h 2 }, the squared mass matrix of CP-even neutral scalars is given as [38] where g 11 = λ 1 cos 2 β + λ 2 sin 2 β + 2(λ 3 + λ 4 ) sin 2 β cos 2 β , Performing another change of basis such that where m h and M H are masses of physical CP-even neutral Higgs bosons. These masses and the mixing angle β − α are computed from the above diagonalization. We assume that m h ≤ M H and identify the lighter state with the observed SM like Higgs. In order to make consistent matching between theory and data, we convert the running mass m h evaluated at the scale M t to the pole mass M h using the following formula where δm 2 h is the SM one-loop self-energy correction and its expression in terms of MS parameters is given in [10,13]. We do not include the contributions from the other scalars in δm 2 h and assume that they are sub-dominant compared to the SM contributions. Numerically we find that the correction to the Higgs mass induced by the second term in the above equation remains less than 0.6 GeV.
The spectrum of physical scalars and the angle β − α are subject to several direct and indirect constraints. We consider the experimentally measured value from [40] and allow a deviation of ±3 GeV from the central value to account for theoretical uncertainty in estimating the value of M h . In THDM of type II, the charged Higgs with mass up to 580 GeV is disfavoured by b → s + γ measurements at 95% confidence level [41] for almost any value of tan β. Further, it can be seen from Eq. (21) that the couplings of h with the weak bosons is proportional to sin 2 (β − α). These couplings are constrained by the signal strength of Higgs decaying into pair of vector bosons. The results of a recent global fit of THDM parameters indicate that the deviation from β − α = π/2 cannot be larger than 0.055 in the case of type II THDM [42]. The above constraints on the spectrum are summarized as We find that the limit on M H ± also puts a lower bound on M A , from Eq. (18), and hence a lower bound on M H as well, since all the quartic couplings are determined from the supersymmetry in this model. The bounds on the masses of THDM scalars obtained in this way are more stringent than the direct search bounds, see for example [42] and references therein. We also investigate the effects of these scalars on electroweak precision observables. For this, we estimate corrections to W boson mass M W and effective weak mixing angle θ lept eff in THDM and compare them with their measured values following the procedure adopted in [43]. These constraints are found to be always satisfied for the values of M A allowed by the constraints listed in Eq. (24). We use the experimental values of gauge couplings and fermion mass parameters measured at different scales and evolve them to the scale M t . This is described in Appendix C. These values are listed in Table I tree-level relations in a convenient basis: where V CKM is quark mixing matrix. For its elements, we use the latest values from the PDG [44]. In the standard parametrization, this matrix is given in terms of three mixing angles and a phase. We use their values: sin θ q 12 = 0.2251, sin θ q 23 = 0.041, sin θ q 13 = 0.0036 and δ CKM = 68.04 • . In the next section, we discuss our procedure for solving RG equations and present the obtained results in details.

IV. NUMERICAL RESULTS
We numerically solve the 2-loop RG equations for different cases and implement matching conditions and constraints as discussed in the previous sections. First, the gauge and Yukawa couplings are evolved from M t to M S using their values at M t and 1-loop RG equations. The running of these couplings do not depend on quartic couplings at 1-loop. The values of quartic couplings at M S are then obtained using conditions given in Eqs. (3,4) and oneloop threshold corrections listed in Appendix B. We then perform full 2-loop evolution of the gauge, Yukawa and quartic couplings form M S to M t and compare the obtained values of gauge and Yukawa couplings with their input values at M t . The couplings are then evolved again from M t to M S using 2-loop RG equations. This procedure is carried out iteratively until the values of couplings converge to their input values at M t . During the running, we check that the stability or metastability conditions, given in Eqs. (12,13), are satisfied at all intermediate scales. Once the convergence is obtained, we calculate the masses of physical scalars using Eqs. (18,19) as function of input parameter M A and tan β.
Before we discuss the results of our numerical analysis, we outline our method of obtaining the RH neutrino mass spectrum in the case of high scale seesaw. The RH neutrino masses are evaluated using the boundary condition given in Eq. (9) and seesaw formula Eq. (7) with the following replacement for M ν : where U PMNS is leptonic mixing matrix and it is parametrized by three mixing angles and three CP phases in the standard parametrization [44]. Since we consider a basis in which charged lepton Yukawa matrix is diagonal 4 , the above form of M ν leads to realistic lepton mixing. We assume normal ordering in the masses of neutrinos, and obtain the values of m ν 2 and m ν 3 in terms of m ν 1 from the solar and atmospheric squared mass differences,  [45]. We also fix the lightest neutrino mass m ν 1 = 0.001 eV for definiteness. Note that we use low energy values of neutrino masses and mixing parameters in seesaw formula to determine the RH neutrino mass spectrum. It is assumed that these parameters do not change significantly under RG evolution. Even if the running effects are taken into consideration for these parameters, the change in the masses of RH neutrinos obtained by Eq. (7) is small and therefore it has negligible effects on the results of our analysis.

A. Without seesaw
We first analyse a case without RH neutrinos. In this case, it is observed in [16] that the quartic couplings λ 1 , λ 2 and λ 3 remain positive at all the scales because of their boundary conditions at M S and hence the first three of the conditions of Eq. (12) are always satisfied. λ 4 is negative at M S and it remains negative at the intermediate scales which in turn requires sufficiently large λ 1,2,3 in order to satisfy the last condition of Eq. (12) or Eq. (13) for stability or metastability respectively. For very small values of tan β, it is observed that the couplings λ 1 and λ 3 evolve slowly. The magnitude of λ 2 however increases rapidly while running from M S to M t because of large and negative contribution to 1-loop beta function from the top quark loop. In this case, it is large and positive λ 2 which makes λ 4 (Q) > 0 and leads to a stable scalar potential. With increasing tan β the top quark Yukawa coupling y t decreases which in turn slows down the running of λ 2 . In this case, λ 2 does not attain large enough value to make the potential stable. For very large values of tan β, the bottom quark Yukawa coupling y b becomes as strong as y t and hence the running effect in λ 1 becomes as strong as that in λ 2 and both are driven to positive and large values at scales below M S . Their combined contributions increase the value of λ 4 (Q) which drive electroweak vacuum towards metastability or stability. We obtain the values of quartic couplings at M t and compute the scalar mass spectrum and impose the constraints given in Eq. (24). The results are displayed in Fig. 1.
We find that the top Yukawa coupling becomes non-perturbative for tan β ≤ 1.36 and hence the perturbative approach of RG evolution breaks down. A stable scalar potential is achieved for 1.36 < tan β ≤ 1.77. As it is explained earlier, the large positive value of λ 2 makes λ 4 (Q) > 0 for all Q between M S and M t in this case. For 1.77 < tan β ≤ 3.18, the electroweak vacuum becomes metastable while the region 3.18 < tan β ≤ 38.5 is disfavoured by an unstable vacuum. The potential again becomes metastable for tan β > 38.5 because of large contribution of bottom quark in the running of λ 1 . It can be seen from Fig. 1 GeV. This together with constraint on Higgs mass restrict the values of tan β in a very narrow range, i.e. 1.36 < tan β ≤ 2.5. We find that our results are in qualitative agreement with the results obtained in [16] but for M S = 2 × 10 17 GeV.

B. Case: HSS
The effects of high scale seesaw mechanism on vacuum stability is studied by incorporating RH neutrino thresholds in the RG evolution and by considering the SO(10) GUT inspired boundary conditions, Eq. (9), at the scale M S . The parametrization is chosen in such a way that ξ and N control inter-generational hierarchy in Y ν . For example, the couplings in Y ν are as hierarchical as those in Y u for ξ ≈ 1. Further, from the extrapolated values of couplings in Y u we find that all the three RH neutrino species couple with the SM leptons with equal strength for N ≈ 2.07 and ξ ≈ 340. The parameter sets overall scale of Y ν and hence the scale of RH neutrinos. For example, the RH neutrino masses obtained using Eqs. (7,9), as explained in the beginning of this section, are displayed Fig. 2 for ξ = 340 and N = 2.07. It can be seen that for small values of tan β, all the neutrinos decouple from the effective theory for ≥ 9. For ∈ [0. 1,9], at least one neutrino remains in the spectrum below M S . We fix N = 2.07 and evaluate the effect of neutrino(s) below M S on the vacuum stability in two different cases. First we fix = 1 and vary ξ in the range from 1 to 400 to investigate the effects of different hierarchical structure of Y ν . In the second, we fix ξ = 340 which corresponds to all the three neutrinos coupled with approximately equal strength and vary in the range 0.1 to 9. The constraints on tan β obtained by the consideration of vacuum stability for these cases are displayed in Fig. 3.
The presence of RH neutrinos below M S modifies the running of quartic couplings in the following way. As it can be seen from the expressions of beta functions given in Appendix A, RH neutrinos contribute to the running of λ 2 with a term proportional to −Tr . This contribution increases the value of λ 2 from its already positive value at M S while evolving from M S to M R i . As a result of this, λ 2 takes relatively larger value at the scales between M R i and M t in comparison to the case without RH neutrinos. This helps in obtaining increased λ 4 (Q) and hence this effect can drive the potential towards stability depending on the magnitude of couplings in Y ν . As it can be seen in the right panel in Fig. 3, this happens only when > 1.5 for which RH neutrinos are very strongly coupled with SM leptons. One finds that all the values of tan β between 1.2 and 50 are allowed by stability constraints for ∈ [1.5, 8]. For > 8 and large tan β, the quartic coupling λ 2 suddenly becomes negative resulting into an unstable region displayed in the upper right corner of tan β-plane in Fig.  3. At two-loop, the contribution from RH neutrino Yukawa couplings in the running of λ 2 comes with opposite sign in comparison to that of one-loop as it can be seen from the expressions of beta functions given in Appendix A. For large and just before the coupling Y ν becomes non-perturbative, the two-loop beta function dominates over the one-loop. This effect is further enhanced by large tan β through large Y d and Y e . This results into large positive value of the total beta function which rapidly drives λ 2 from its initial positive value at M S to negative value at the scale below M S . For ≤ 1.5, the heavy neutrinos do not significantly change the stability constraints obtained in the case without RH neutrinos as it can be seen from both the panels in Fig.  3. Although all the neutrinos are not decoupled in this case, they do not have couplings strong enough to make any significant change in the running of quartic couplings. Even for ξ = 340, N = 2.07 which leads to all the three neutrinos with couplings of O(y t ), one obtains similar constraints on tan β as in the case without RH neutrinos.
In order to check the viability of above results with Higgs mass and other low energy constraints, we evaluate the scalar spectrum for three different benchmark points. These are: (i) = ξ = 1, (ii) = 1, ξ = 340 and (iii) = 4, ξ = 340. The choices (i) and (ii) are motivated from a consideration that at least one or all the three neutrinos have couplings as large as O(y t ) while (iii) corresponds to very strongly coupled neutrinos which lead to significant improvements in stability constraints. The results are displayed in Fig. 4. For = ξ = 1 (or equivalently Y ν = Y u at the GUT scale), the results are almost identical to the ones obtained in the case without RH neutrinos as it can be seen from Figs. 1 and 3. The results do not change significantly also for = 1 and for any value of ξ between 1 and 340. For these benchmark points, the stability or metastability of electroweak vacuum and constraints on the masses of scalars imply 1.2 ≤ tan β ≤ 2.2 and M A > 580 GeV as viable ranges in which the effective THDM theory can be extrapolated to M S . > 1.5 leads to stability of the scalar potential for almost all the values of tan β, however they face stringent constraints from M h and M H ± as it can be seen from the bottom panel in Fig. 4. Strongly coupled neutrinos in this case increase the magnitude of λ 2 at M t which results into relatively higher Higgs mass for a given value of M A . Hence to obtain M h < 128 GeV, one needs lower M A which is already disfavoured by the constraints on the charged Higgs mass. This disfavours the values of > 1. The above results imply that an existence of strongly coupled heavy neutrinos below M S improves the stability of THDM scalar potential, however this scenario is very strongly constrained by the observed Higgs mass and branching ratio of b → s + γ. in Eq. (10). The orange region corresponds to unstable scalar potential while the grey region in the bottom is where the perturbativity of the couplings is lost.

C. Case: LSS
We perform a similar analysis for a low scale seesaw scenario discussed in the section II. For simplicity, we assume that y 1 = y 2 = y 3 ≡ y ν in Eq. (10). The third RH neutrino does not couple to THDM as it can be seen from the flavour structure of Y ν in Eq. (10). We take two sample values for mass of the remaining degenerate heavy neutrinos, M = 10 3 and 10 9 GeV. The results are displayed in Fig. 5. As it can be seen, strongly coupled neutrinos lead to significant changes in the stability of the scalar potential. The quartic couplings receive significant contributions from such neutrinos in this case which in turn helps in making the scalar potential more stable. The stability constraints on y ν and tan β do not depend considerably on the mass scale of RH neutrinos.
We analyse the low energy scalar spectrum for three benchmark values of coupling: y ν = 0.1, 0.3 and 0.5. The results are displayed in Fig. 6. For y ν ≤ 0.1, the obtained constraints on tan β and M A are very similar to those obtained in the case without RH neutrinos. In the regime of strong y ν , these results change significantly. For y ν = 0.3, one always obtains M h > 125 GeV if constraint on the charged Higgs mass is considered. Very strongly coupled neutrinos further push the value of M h on the higher side and such cases are disfavoured by limits on M H ± and β − α.
In the analysis discussed above, we consider the top pole mas M t = 173.5 GeV. The change in top mass is known to have considerable effects on the RG evolution. To evaluate

V. CONCLUSION
We have investigated the viability of two-Higgs-doublet model, with type I seesaw mechanism incorporated in it, as an effective theory below the GUT scale in the presence of its minimal supersymmetric completion at the same scale. Supersymmetry is assumed to be broken at the scale M S = 2 × 10 16 GeV leaving only the SM particles, three copies of RH neutrinos and an additional Higgs doublet below M S . It also leaves its imprints on the effective scalar potential of THDM by relating quartic couplings with gauge couplings at the SUSY breaking scale. Because of these constraints provided by SUSY, the effective potential of THDM can be written in terms of only two free parameters: the ratio of vacuum expectation value of two Higgs doublets tan β, and the mass of pseudo-scalar Higgs M A . The same potential governs the electroweak symmetry breaking and determines the masses of physical scalars in the theory. The requirement that the electroweak vacuum is stable or metastable (i.e. its lifetime is greater than the age of universe) and constraints on the scalar spectrum from various direct and indirect experimental searches restrict the allowed values of tan β and M A . In the absence of RH neutrinos, these constraints favour the values for tan β in the range from 1.4 to 2.5 and M A > ∼ 580 GeV. Existence of RH neutrinos below the scale M S causes significant changes in the above results if they are very strongly coupled with the SM leptons. We have studied this possibility in the context of high and low scale type I seesaw mechanism. In the case of high scale seesaw, we find that if one or more neutrinos have Dirac type Yukawa couplings greater than y t at the GUT scale, then they considerably modify the running of quartic couplings and allow stable or metastable electroweak vacuum for almost all the values of tan β for which the effective theory remains perturbative. The RH neutrinos improve the stability of scalar potential in this case unlike in the case of the SM without high scale SUSY in which the strongly coupled RH neutrinos are known [9,[46][47][48][49][50][51][52] to destabilize the scalar potential. The difference between the two cases arise mainly because of the presence of SUSY at the high scale which ensures the stability of scalar potential in ultraviolet completion. The difference is also attributed to the presence of additional Higgs doublet in THDM which modifies the stability conditions allowing more freedom in the effective potential. We find that the observed Higgs mass and limit on the charged Higgs mass from the flavour physics data put stringent constraints on tan β and M A when RH neutrinos are strongly coupled. Similar results are obtained in case of low scale seesaw when neutrino Yukawa couplings are O(1).
It is observed that if Yukawa couplings of RH neutrinos are of O(y t ) or smaller at the GUT scale then they do not have significant impact on vacuum stability constraints in THDM. In the case of high scale seesaw, particular cases have been explored in which it is assumed that the neutrino Yukawa coupling matrix Y ν is equal to up-type quark Yukawa coupling matrix Y u or all the couplings in Y ν are O(y t ) at the GUT scale. This class of boundary conditions are often realized in SO(10) based GUTs. It is found in these cases that RH neutrinos have negligible effect on stability of the scalar potential and scalar spectrum at low energy. A stable or metastable scalar potential consistent with low energy constraints is obtained for tan β ∈ [1.2 − 2.5] in these cases. These results are qualitatively similar to those obtained in case without right-handed neutrinos in [16].
The results obtained in this paper are useful in constraining a class of models in which supersymmetry is broken at the GUT scale and an effective theory below this scale is THDM of type II with type I seesaw mechanism. It is shown here that the stability or metastability of electroweak vacuum and a consistent low energy scalar spectrum are achieved for only small tan β and for neutrino Yukawa couplings of O(y t ) or smaller. In concrete ultraviolet models, see for example the ones given in [25,26], the neutrino Yukawa couplings and/or tan β are often determined from the enhanced symmetry structure of underlying theory. Therefore, the results obtained in this paper can be used to constrain such frameworks.
Our analysis provides a generic understanding of the effects of RH neutrinos. The numerical results obtained in this paper are subject to change for different choice of values of M S and/or GUT scale threshold corrections. We have assumed a particular value of SUSY breaking scale, M S = 2 × 10 16 GeV, throughout our analysis. It is known that stability of vacuum potential is considerably sensitive to the choice of such a scale [16]. In particular, the constraints on the values of tan β imposed by unstability of potential become feeble (stronger) for relatively smaller (larger) values of M S . Another important issue is that the precise unification of gauge couplings does not occur at the scale M S as the effective theory below this scale is pure THDM with singlet neutrinos. Exact gauge coupling unification may require either new fields at intermediate scales (for example, a pair of TeV scale Higgsinos [16]) or sizeable threshold corrections at the GUT scale [11]. For all these cases, a dedicated analysis would be required in order to derive quantitative constraints on tan β and M A . However, we anticipate that the qualitative effects of RH neutrinos will remain the same. The presence of one or more strongly coupled RH neutrinos below M S improves the stability of vacuum in an underlying framework. In this appendix, we provide 2-loop renormalization group equations for THDM of type II. They are obtained using publicly available package SARAH [37]. Note that we use different convention for the quartic couplings λ 1 and λ 2 in Eq. (1) in comparison to the one used in SARAH. The RG equations listed below are therefore modified accordingly. The same equations are also listed in [13] but with considering only the third generation of fermions.
The couplings evolve according to the following equation: where C represents gauge, Yukawa and quartic couplings and µ is the renormalization scale. The one and two-loop beta functions for the different couplings are as the following.