Pseudo-Dirac Higgsino dark matter in GUT scale supersymmetry

We investigate a scenario in which supersymmetry is broken at scale $M_S \geq 10^{14}$ GeV leaving only a pair of Higgs doublets, their superpartners (Higgsinos) and a gauge singlet fermion (singlino) besides the standard model fermions and gauge bosons at low energy. The Higgsino-singlino mixing induces small splitting between the masses of electrically neutral components of Higgsinos which otherwise remain almost degenerate in GUT scale supersymmetry. The lightest combination of them provides a viable thermal dark matter if the Higgsino mass scale is close to $1$ TeV. The small mass splitting induced by singlino turns the neutral components of Higgsinos into pseudo-Dirac fermions which successfully evade the constraints from the direct detection experiments if the singlino mass is $\lesssim 10^8$ GeV. We analyse the constraints on the effective framework, arising from the stability of electroweak vacuum, observed mass and couplings of Higgs, and the limits on the masses of the other scalars, by matching it with the next-to-minimal supersymmetric standard model at $M_S$. It is found that the presence of singlino at intermediate scale significantly improves the stability of electroweak vacuum and allows a stable or metastable vacuum for almost all the values of $\tan\beta$ while the observed Higgs mass together with the limit on the charged Higgs mass favours $\tan\beta \lesssim 3$.


I. INTRODUCTION
If supersymmetry (SUSY) does not stabilize the electroweak scale then the scale of its breaking is not restricted to stay in the vicinity of the weak scale. Most of the other attractive features of low energy supersymmetry, such as precision gauge coupling unification and particle candidate(s) of weakly interacting dark matter (DM), could also be achieved in scenarios like split-supersymmetry [1,2] in which only some of the superpartners are required to be close to the weak scale. The essential role played by SUSY in superstring theories [3], which are potential candidates for the unification of all the fundamental forces, remains the same as long as it is broken at any scale below the string scale. In fact, it is typically expected that the SUSY breaking scale in such theories is close to the string scale. Similarly, in a class of models based on supersymmetric grand unified theories (GUT) in higher spacetime dimensions (see for example, [4][5][6][7][8][9][10][11][12][13]), the breaking of SUSY and unified gauge symmetry are often administered by a common mechanism which gives rise to the breaking of SUSY at the GUT scale.
In the post Higgs discovery era, scenarios of high scale SUSY are constrained by the measured Higgs mass and stability of the electroweak vacuum. For example, it is known that the standard model (SM) cannot be matched with its simplest supersymmetric counterpartthe minimal supersymmetric standard model (MSSM) -if the SUSY breaking scale is above 10 10 GeV because of the constraints arising from the stability of electroweak vacuum [14][15][16][17]. The scale can be raised up to the GUT scale if SM is replaced by the two-Higgs-doublet model (THDM) as a low energy effective theory [18]. The presence of additional Higgs doublet helps in achieving a stable electroweak vacuum in this case. The stability can further be improved if there exists strongly coupled right handed neutrinos below the GUT scale [19]. An effective theory involving the THDM and a pair of Higgsinos at the weak scale as remnants of GUT scale supersymmetry has also been studied in [18,20] and it is shown compatible with the vacuum stability constraints [18]. While the weak scale Higgsinos improve the convergence of gauge couplings, they can also provide a potentially viable candidate of DM given an unbroken R-parity in the underlying supersymmetric theory. However, the second possibility is disfavoured by the null results from the DM direct detection experiments.
The pure Dirac Higgsino DM is already ruled out because of its large elastic scattering cross section with nucleons [21]. The mixing of Higgsinons with bino, wino and/or the radiative corrections from supersymmetric particles can induce splitting between the masses of neutral components of Higgsinos making them Majorana fermions and hence the above constraints can be evaded. However, the required amount of mass splitting is obtained only if the mass scale of the other supersymmetric particles is < ∼ 10 8 GeV [22]. In this paper, we discuss a possibility in which pseudo-Dirac Higgsino DM can be reconciled with GUT scale supersymmetry. We consider an effective theory consists of THDM and a pair of Higgsinos augmented by a singlet fermion, namely the singlinoS. The effective theory is matched with the next-to-minimal supersymmetric standard model (NMSSM) at the SUSY breaking scale M S which is assumed close to the GUT scale. The assumed hierarchy among the different scales and corresponding mass scales of particles are depicted in Fig. 1. We find that the Higgsino-singlino-Higgs Yukawa interaction can induce required mass splitting between the neutral components of Higgsinos if the singlino mass scale is < ∼ 10 8 GeV. The same interaction also improves the stability of electroweak vacuum significantly. It is shown that viable Higgsino DM with stable or metastable electroweak vacuum and scalar spectrum consistent with the current experimental observations can be obtained within the framework of GUT scale supersymmetry.
The paper is organized as follows. In the next section, we outline the effective framework and discuss its matching with NMSSM at the SUSY breaking scale. We provide details of various phenomenological constraints applicable on the underlying framework in section III. The details of numerical analysis are described in section IV followed by the results and discussion in section V and summary in section VI. Some relevant technical details are provided in three Appendices.

II. FRAMEWORK
We consider an effective theory below the SUSY breaking scale described by the following renormalizable Lagrangian: The first term denotes the Lagrangian of the most general two-Higgs-doublet model which contains a pair of weak doublet scalars H 1,2 with hypercharge Y = 1/2. The scalar potential and Yukawa interactions terms in L THDM can be written as and respectively. Here i, j = 1, 2, 3 denote three generations of SM fermions and H c 1,2 = iσ 2 H * 1,2 . The second term in eq. (1) represents free Lagrangian of a pair of fermions, namelyH u andH d , which are SU (2) L doublets with Y = 1/2 and −1/2, respectively. where The fieldsH u,d have the same gauge quantum numbers as those of Higgsinos in the MSSM [23] and their explicit forms are: Similarly, LS contains a Majorana mass term of the singlinoS: The last term in eq. (1) includes gauge invariant renormalizable Yukawa interactions involving the fermionsH u ,H d ,S and scalars H 1,2 . It is given as: The free Higgsino Lagrangian LH possesses a global U (1) symmetry under whichH u andH d have equal and opposite charges. However, this symmetry is broken either by the Yukawa interactions in eq. (7) or by the Majorana mass term in eq. (6) given H 1,2 remain uncharged under this symmetry. As a result, the Higgsinos become pseudo-Dirac fermions in this framework. It can also be noticed that L possesses a Z 2 symmetry under which the fieldsH u,d andS are odd while all the other fields are even. Origin of this symmetry in the effective theory can be attributed to the presence of unbroken R-parity in the underlying UV supersymmetric theory.
A. Matching with NMSSM at M S A minimal supersymmetric setup which can provide UV completion of an effective theory, described by L in eq. (1), is the well-known NMSSM. We therefore match the effective theory with NMSSM at the SUSY breaking scale and obtain constraints on various parameters. The most general superpotential of NMSSM is given as [24] where W MSSM is the standard MSSM superpotential (see for example [23]) andŜ is chiral superfield with singlinoS and a complex scalar S as its submultiplets. The corresponding soft SUSY breaking sector can be parametrized in terms of the following potential.
The trilinear scalar couplings in V soft are assumed to be proportional to the corresponding Yukawa couplings in the superpotential. The term linear in S is known to generate potentially dangerous quadratic divergences in supergravity [25] and therefore the coupling c S needs to be suppressed. Typically, this is achieved by introducing a symmetry which forbids such term. In the most popular versions of the NMSSM, the same symmetry is often utilized to solve the so-called µ problem, see [24,26] for examples. In our setup, we assume for brevity and consider the remaining parameters µ, µ S , λ, A λ and b S as real and positive. A phenomenologically consistent pseudo-Dirac Higgsino DM requires both µ and µ S well below the scale M S . We also assume b S ∼ O(µ 2 S ) and m S > ∼ M S which lead to b S M 2 S . We then compute the complete scalar potential of the theory from V soft and the D and F terms of superpotential W. After integrating out the singlet scalar S from the theory and keeping only the leading order terms in b S /M 2 S , the resulting effective potential is matched with the THDM potential given in eq. (2) using the identities H 2 = H u and H 1 = −iσ 2 H * d . This implies the following tree level matching conditions at the SUSY breaking scale M S : wherem S ≈ m 2 S + µ 2 S is physical mass of S. The parameters m 1 and m 2 in eq. (2) are generated by the soft masses of H d and H u respectively, while m 12 arises from the b µ term in V MSSM soft in eq. (9). Further, the matching between the Yukawa interactions in W and those in eq. (3,7) leads to y 1 = y 2 = λ , y 3 = y 4 = 0 , at M S .
It can be noticed from eq. (11) that the precise values of λ i at M S still depend on the details of SUSY breaking sector despite of the simplification achieved through ansatz in eq. (10). Withm S > ∼ M S , the couplings λ 5 and λ 6,7 are found to be suppressed by a factor of at least µ S /M S and µ/M S , respectively. Further suppression could be obtained if |A λ | M S . One finds that λ 6,7 and m 12 have vanishing values in the limit µ, b µ → 0. Together with conditions in eqs. (12,13), this implies Z 2 symmetry of L in eq. (1) under which the fields H 1 ,H d ,S, d i R and e i R are odd. The symmetry keeps the parameters which vanish at M S under control and prevents them from taking large values through renormalization group evolution (RGE) effects. In the effective theory, the Z 2 symmetry is softly broken by the terms proportional to µ and m 12 . Our assumption of real µ, µ S , λ, A λ , b µ and b S implies real values for couplings in eq. (2) 1 . The effective theory below M S , obtained from the NMSSM with the aforementioned conditions, resembles the well-known CP conserving type II version of THDM [27]. The effective theory, therefore, does not give rise to low energy flavour or CP violating effects additional to those already anticipated in type-II THDM.

B. Higgsino dark matter
One of the main aims of this study is to show the existence of viable pseudo-Dirac Higgsino dark matter. We therefore discuss the Higgsino mass spectrum in detail. As already emphasized, the Higgsinos and their interactions in the effective theory below M S are well described by the terms in eqs. (4,6,7) with y 3,4 = 0. Further integrating out the singlino from the low energy spectrum, the effective Higgsino Lagrangian at scale below µ S is given by with the following matching conditions at µ S : The electroweak symmetry breaking then contributes into the masses for the neutral and charged components ofH u,d .
In the basis, N = H 0 d ,H 0 u T , the mass term for neutral components can be written as with and one obtains the following expressions for the tree level neutralino masses with a convention mχ0 1 < mχ0 2 : As it can be seen from eqs. (12,15), c 1,2 are real and positive for real values of λ. The unitary matrix representing neutralino mixing is Further, it is convenient to define the splitting between the masses of neutralinos as The charged components ofH u,d can be combined to form a Dirac fermionχ where mχ± = µ is mass of chargino at the tree level. As it can be seen from eq. (20), the contributions induced by dim-5 operators generates splitting between the masses of chargino and neutralino. In addition, it is known that loop corrections induced by the SM electroweak bosons make the chargino heavier than the neutralino [28]. The resulting mass splitting between the chargino and the lightest neutralino can be written as where ∆m rad ± is radiatively induced mass splitting. For µ ≈ 1 TeV, one obtains ∆m rad ± ≈ 341 MeV at one loop [28]. The contribution from the first two terms in eq. (24) remains positive for real λ.

III. PHENOMENOLOGICAL CONSTRAINTS
We now discuss the set of constraints which are imposed on the effective framework described in the previous section.

A. Dark matter
The nature of the lightest neutralinoχ 0 1 is almost pure Higgsino like for µ S µ in this framework. Assuming that it makes all of the DM produced thermally in the Early Universe, its relic abundance is estimated in [28][29][30] including the non-perturbative Sommerfeld corrections to DM annihilation cross sections. It is found that the observed relic abundance is obtained 2 if the mass of DM particle is close to 1 TeV. This implies also µ ≈ 1 TeV in the present framework.
Experiments based on the direct detection of DM are known to put stringent constraints on pure Dirac Higgsino DM. If ∆m 0 = 0 thenχ 0 1 andχ 0 2 can be paired to form a Dirac fermion, namelyχ. Asχ has vector coupling with Z boson, the elastic scattering with nucleon N (such asχ + N →χ + N ) can proceed through the exchange of the Z boson. The scattering cross section, which is unambiguously determined for a given mass of Higgsinos, turns out too large and therefore this case is disfavoured by the non-observation of any statistically significant signal in the direct detection experiments [21]. Nevertheless, this constraint can be evaded if ∆m 0 > 0. In this case,χ 0 1 andχ 0 2 are Majorana fermions and hence they do not have vector coupling with Z boson.
The inelastic scattering,χ 0 1 + N →χ 0 2 + N , is still subject to the constraints from the direct detection experiments for Majoranaχ 0 1,2 , if ∆m 0 is very small. Such processes arise through the t-channel exchange of Z boson. Considering this, a lower limit on ∆m 0 has been derived in [22] using the then available data from XENON 10 and XENON 100 experiments. We update this analysis for the latest available data, including those from XENON 1T. The details are described in Appendix A 1. We find that the present observations from the direct detection experiments lead to a lower bound on neutralino mass splitting, ∆m 0 ≥ 200 keV, at 90% confidence level. In the present framework, this bound translates into an upper limit on the singlino mass scale µ S . Using eq. (20) and assuming c i as O(1) numbers, we find A similar bound on gaugino mass scale was obtained earlier in the case in which the mixing of Higgsino with bino or wino was responsible for mass splitting [22,31]. Higgsino mass splitting can also be induced radiatively through stop-top loop if the stop mixing angle is nonzero [32]. This however also requires the stop masses < ∼ 10 8 GeV for ∆m 0 > 200 keV.
In our framework, all the super-partners can have GUT scale masses while the presence of singlino, with µ S < ∼ 10 8 GeV, can induce the required ∆m 0 . The spin-independent and spin-dependent elastic scattering processes, likeχ 0 1 +N →χ 0 1 + N , can also occur in the underlying framework through the exchange of THDM scalars or Z boson, respectively. We find that the scattering cross sections of the first type of interactions are suppressed due to µ S µ while those of the latter are negligible because of pseudo-Dirac nature of Higgsino DM. After determining the constraints on the scalar spectrum and couplings, we estimate these cross sections and show that the obtained results are in agreement with the current experimental limits. This is described in detail in Appendix A 2.
The Higgsino DM can give rise to indirect signatures through their pair annihilation into W + W − at tree level and ZZ, γγ, Zγ at loop level, which subsequently leads to the production of gamma rays and anti-protons. The latest constraints on almost pure Higgsino DM from the indirect searches are reviewed in [33,34]. The constraints from the observations of gamma-ray by Fermi-LAT [35] exclude Higgsino DM with mass ≤ 330 GeV while HESS [36] observations do not put any such limit. The strongest indirect search constraints on Higgsino DM arise from the latest AMS-02 results [37] on anti-protons which put a conservative limit mχ0 1 ≥ 500 GeV. Nevertheless, the almost pure Higgsino DM with mass ≈ 1 TeV considered in our framework remains unconstrained from the current indirect detection experiments.

B. Vacuum stability and perturbativity
The matching conditions obtained in eq. (11) determine the values of quartic couplings in terms of the gauge couplings and several of NMSSM parameters. With the assumed hierarchies in the scales, i.e. µ < µ S m S ∼ M S , the couplings λ 5 and λ 6,7 are suppressed by factors of O(µ S /M S ) and O(µ/M S ), respectively. They are even more suppressed if A λ M S and/or b S ≈ M 2 S . Further, the approximate Z 2 symmetry of the effective theory forbids them from taking large values through running effects. Therefore, the contributions of the terms involving the couplings λ 5,6,7 in the scalar potential remain negligible at all the scales below M S . The remaining quartic couplings must satisfy the following conditions, at the intermediate scales between M S and M t , for the absolute stability of electroweak vacuum [38] In a more conservative approach, it is assumed that the scalar potential in eq. (2) may have multiple minima and the electroweak vacuum can decay into a more stable minimum through quantum tunnelling. However, the lifetime of electroweak vacuum is required to be greater than the current age of the Universe in order to make it phenomenologically viable. Such a minimum of potential is termed as metastable vacuum and its lifetime, in case of single scalar field, is estimated in [39] including the quantum effects. This was generalized in [18] for THDM scalar potential after mapping the underlying potential into single filed potential using the first three of the conditions given in eq. (26). The requirement of metastable vacuum replaces the last condition in eq. (26) with [18] where Q is the renormalization scale. The first three of the conditions listed in eq. (26) are found to be satisfied as a consequence of the high scale boundary conditions given in eq. (11). Note that with µ < µ S m S ∼ M S , the couplings λ 1,2 are positive while |λ 3 | √ λ 1 λ 2 [18,19]. Hence, it is the last condition in eq. (26) or eq. (27) which determine the stability or metastability of scalar potential, respectively. The same has been the case for THDM matched with MSSM at the GUT scale [18,19] in which the running of λ 4 dominantly decides the stability of vacuum. However, two important differences arise in the present framework. Firstly, the matching with NMSSM modifies boundary condition for λ 4 as can be seen from eq. (11). Unlike in the case of MSSM, λ 4 can be made positive at M S with appropriately chosen values of λ and A λ . Secondly, even if the tree level enhancement in λ 4 is absent (for example, if A λ ≈ M S ), the contribution from singlino-Higgsino loop modifies the running of λ 4 and drives its value toward the positive side while running from M S down to the M t as it can be seen from appropriate RG equations given in Appendix B.

C. Scalar spectrum
The matching of the effective theory with the NMSSM determines the values of all the quartic couplings which in turn provides useful correlations among the masses of THDM scalars. In order to check the viability of such correlations, we evaluate the mass spectrum of these scalars and impose various experimental constraints on the remaining free parameters of the potential. We closely follow the procedure and notations of [19] which is briefly outlined in the following.
The VEVs of the neutral components of H 1 and H 2 break the electroweak symmetry giving rise to five physical Higgs bosons: two CP even and electrically neutral (h, H), two CP even and charged (H ± ) and a CP odd and neutral (A). At the minimum, the parameters m 1 , m 2 and m 12 in eq. (2) can be replaced by appropriate functions of M A , tan β, v and quartic couplings [38]. Here, M A represents the physical mass of pseudo-scalar Higgs in MS renormalization scheme while tan β ≡ v 2 /v 1 . This replacement allows us to express the masses of CP even neutral and charged Higgs bosons in terms of the unknown parameters M A , tan β and known parameters v and λ i . The physical CP even neutral Higgs bosons can be obtained as linear combinations of neutral components in H 1 and H 2 such that H = cos α H 1 + sin α H 2 and h = − sin α H 1 + cos α H 2 . The state h is identified with the observed SM like Higgs and H is assumed to be heavier than h. For a consistent matching between the theoretically predicted mass of h and that of the observed Higgs, we covert the running mass into pole mass, namely M h , using the method described in [19]. The mixing angle α and the masses of heavy CP even Higgs (M H ) and charged Higgs (M H ± ) are also determined in terms of M A , tan β, v and λ i . Their expressions are given in [19] with all the necessary details.
We consider the following set of constraints on the masses of Higgs bosons and the mixing angle α.
For M h , the experimentally measured value from [40] is considered. In order to account for theoretical uncertainties in estimating the Higgs mass, we allow a deviation of ±3 GeV from the measured mean value. Further, the couplings of h with W ± and Z bosons are proportional to sin 2 (β − α) and therefore they are constrained from the observed signal strengths of h → W + W − and h → ZZ. This, through the results of the latest global fit performed in [41], implies that the deviation from the so-called alignment limit, i.e. β − α = π/2, cannot be greater than 0.055 in the case of THDM of type II which gives rise to the second constraint in eq. (28). The quoted lower bound on the mass of charged Higgs boson arises form the observed branching ratio of b → s + γ at 95% confidence level [42]. Note that this bound is applicable for almost all values of tan β in THDM of type II.
Since the masses of different scalars are correlated in the framework under consideration, we find that the lower limit on M H ± translates into lower limits on M A and M H , thus making all these scalars heavier than 580 GeV. Such heavy scalars already satisfy the direct search bounds and limits arising from the flavour observables such as B s → µ + µ − , see [41] for example and references therein. We also observe that for M A ≥ 580 GeV, the masses of the scalars H, A and H ± remain approximately degenerate in this framework. It is found that such a heavy and degenerate spectrum of scalars always satisfies the constraints imposed by the electroweak precision observables [19,43].

IV. NUMERICAL ANALYSIS
The viability of the proposed framework with respect to the various constraints, discussed in the previous section, is investigated by solving two-loop RGE equations numerically and implementing 1-loop corrected matching conditions. The one and two loop beta functions for the underlying framework are generated using publicly available package SARAH [44] and are listed in Appendix B. We obtain the values of gauge and Yukawa couplings at the top quark pole mass scale M t from the experimental values of gauge couplings and fermion masses measured at the different scales. The procedure of this with relevant details is described in our previous work [19]. The obtained values of gauge couplings and fermion masses at the scale M t are listed in Table I. The values of Yukawa couplings at M t are extracted from the fermion masses as described in [19]. For M t , we use the latest PDG average M t = 173.1 ± 0.9 GeV [45]. Following a similar procedure as described in [19], we first evolve the gauge and Yukawa couplings from M t to M S using one loop RGE equations as the quartic couplings do not contribute in their running at this level. We then obtain the quartic couplings at M S using the matching conditions given in eq. (11). The tree level matching of λ i should be corrected by one loop threshold corrections which are explicitly given in [20,46]. However, such corrections depend on the exact details of the SUSY spectrum and hence require explicit details of SUSY breaking. For simplicity, we assume all the squarks, sleptons and gauginos to be degenerate with mass ∼ M S and vanishing trilinear parameters. Although such a choice of SUSY spectrum is very specific, it is naturally realized in the models of SUSY breaking based on flux compactification, see for example [11]. With this choice of SUSY spectrum and from the fact that µ M S , the one loop threshold corrections in Yukawa couplings are found to be suppressed by the degeneracy in the masses of the superpartners or by the smallness of µ/M S as it can be seen from the relevant expressions given in [20]. The threshold corrections in quartic couplings are also suppressed by either vanishing trilinear couplings or (µ/M S ) 2 and hence they are negligible in the present framework.
After obtaining the values of quartic couplings at M S as described in the above, all the gauge, Yukawa and quartic couplings are evolved down to µ S using full two loop RGE equations. At µ S , we obtain the coefficients c i and d using the conditions given in eq. (15). After integrating out the singlino at this scale, the gauge, Yukawa and quartic couplings are evolved from µ S to M t using appropriate two loop RGE equations. All these couplings are again evolved from M t to M S and back to M t iteratively until the couplings converge to their input values supplied at M t . The stability and metastability of the potential is checked at the intermediate scales using the conditions given in eq. (26) and eq. (27), respectively. After the convergence is achieved, we calculate the masses of scalars and the Higgs mixing angle α as function of input parameters M A and tan β using the obtained values of quartic couplings at M t . We also evolve the effective operators given in eq. (14) from µ S to µ using the one loop beta functions: where β X = dX/d(ln Q). We have derived the above equations following a procedure similar to the one given in [47,48] in the context of neutrinos. The neutralino mass spectrum is then obtained by substituting the values of c i and d at the scale µ in eq. (17). We set µ = 1 TeV, as required by the observed DM abundance, and evaluate the vacuum stability constraints on the values of tan β and λ, for M S = 2 × 10 14 GeV and two sample values of µ S , simultaneously checking their viability to produce large enough ∆m 0 . Since A λ has direct implication to the boundary value of λ 4 , we perform this analysis for A λ = 0 and A λ = M S . We repeat all these cases for M S = 2 × 10 16 GeV as well. Once the consistent values of λ are found, we select some benchmark points and evaluate the low energy spectrum and constraints on M A and tan β. The results of the numerical analysis are discussed in the next section.

V. RESULTS AND DISCUSSIONS
The constraints on tan β and λ arising from the vacuum stability, perturbativity and neutralino mass difference are displayed in Figs. 2 and 3 for two example values of M S . In all these cases, the perturbativity of the effective theory disfavours λ > ∼ 1 and tan β < ∼ 1.2. The latter leads to non-perturbative values for the top quark Yukawa coupling. We find that the bound from direct detection experiment, ∆m 0 > 200 keV, implies λ ≥ 0.03 (0.3) for µ S = 10 5 (10 7 ) GeV for all the values of tan β. From the obtained results, we find an approximate relation: which is a direct consequence of eq. (25). In this case, the perturbativity constraint alone implies a robust upper bound, µ S ≤ 10 8 GeV, for a viable pseudo-Dirac Higgsino DM within this framework. The presence of Higgsino-singlino-Higgs Yukawa interaction also has interesting implications on the stability of the electroweak vacuum. It is known that in the case of pure THDM matched with MSSM at the GUT scale, only a narrow range in tan β ∈ [1.1, 1.8] is allowed by the absolute stability of the electroweak vacuum [18,19]. This conclusion is changed in the present framework for large enough values of λ as it can be seen from the stability, metastability and instability contours displayed in Figs the stability of potential in two ways. For A λ = 0, λ directly modifies the boundary value of λ 4 , making it positive for λ 2 > g 2 2 /2, and hence rescues the electroweak vacuum from instability. This direct effect becomes feeble if A λ ≈ M S . However, the contribution of λ in the running of λ 4 is still able to improve the stability of the vacuum. This indirect effect requires relatively larger values of λ in order to achieve metastability or absolute stability in comparison to the case with A λ = 0. We find that change in µ S , within the range 10 5 -10 8 GeV, has negligible effects on the results of vacuum stability. From Figs. 2 and 3, it can also be observed that the instability regions grow when M S is increased because of the relatively longer running.
We now discuss the consequences of the low energy constraints, listed in eq. (28), on the parameters of the underlying framework. For this, we select two benchmark values of λ, allowed by ∆m 0 ≥ 200 keV, from each of the panels in Fig. 2 and obtain the scalar spectrum at M t as function of parameters M A and tan β. The results are displayed in Figs. 4 and 5. As it can be seen, λ ≥ 0.5 for A λ = 0 is disfavoured by the low energy constraints. In this case, large value of λ results into relatively heavier SM like Higgs through λ 4 [19]. The observed range in M h requires either tan β < 1 or M A < 400 GeV. The first is ruled The discovery prospects of light Higgsinos at colliders, with the other SUSY particles much heavier, have been widely studied, see for example [22,34,[49][50][51][52][53][54][55]. In the absence of coloured superpartners, the Higgsinos are produced through only weak interactions and therefore their production rate is relatively smaller. Further, the collider searches are known to be highly insensitive if ∆m 0 < 5 GeV [50,52]. The decay of chargino into neutralino and charged pion is also difficult to detect if ∆m ± > ∼ 300 MeV [22,56]. Therefore, the present framework seems to remain unconstrained from the current collider searches of Higgsinos. It is more feasible to probe almost-degenerate Higgsinos at future lepton colliders with sufficient centre-of-mass energy through tagging of the photon produced in association with a pair of neutralinos or charginos [49]. The future electron-proton colliders may also provide a potential platform for the searches of almost-degenerate Higgsinos [57]. A dedicated analysis in these directions, for the ranges of masses and couplings of Higgsinos obtained in the underlying framework, would be useful.
Although the underlying framework is shown to accommodate TeV scale pseudo-Dirac Higgsino DM consistent with GUT scale SUSY breaking, it is not complete in all aspects. One of the pertinent issues with the GUT scale SUSY breaking is that it does not lead to precision unification of the SM gauge couplings. The convergence of the couplings, at one-loop, requires [58,59] where b i for i = 1, 2, 3 are effective beta function coefficients of U (1), SU (2) L and SU (3) C gauge couplings, respectively (see [59] for more details). One obtains R ≈ 0.528 and 0.556 for the SM and THDM, respectively. The THDM together with a pair of weak scale Higgsinos lead to R = 0.673, bringing it closer to the value as required in eq. (31). Therefore it improves the convergence of the gauge couplings significantly. However, the scale of unification turns out to be of the O(10 14 ) GeV which is one or two orders of magnitude lower than what is naively expected from the latest experimental limit on the proton lifetime [60]. In this case, achieving the precise unification and/or satisfying constraints from the proton lifetime require new effects at the scale M S or below. This may change the derived results depending on the exact nature of these new effects. Further, if the type-I seesaw mechanism is invoked for non-vanishing neutrino masses then the presence of singlet neutrinos may also modify the vacuum stability constraints if they are strongly coupled [19]. All these cases require dedicated analysis similar to the one presented in this paper in order to obtain a precise limit on the parameters of the underlying framework.

VI. SUMMARY
There exists a class of models in which supersymmetry is broken at the GUT scale giving rise to an effective theory which contains the SM with an additional Higgs doublet and a pair of Higgsinos. If R-parity is unbroken, the lightest of the neutral components of Higgsinos can be dark matter candidate. The observed thermal relic abundance can be achieved if its mass is ∼ 1 TeV. However, in the absence of any other new physics below the GUT scale, the neutral components of Higgsinos are almost degenerate with very tiny mass splitting between them, ∆m 0 ≈ O(0.1) eV. Such an almost pure Dirac-like Higgsino DM is already disfavoured by the direct detection experiments.
The problem can be circumvented minimally by introducing a gauge singlet Majorana fermion, namely the singlino, which is odd under R-parity. The Higgsino-singlino-Higgs Yukawa interaction gives rise to the Higgsino-singlino mixing through electroweak symmetry breaking, which induces large enough ∆m 0 to evade the direct detection constraints. We show that the same Yukawa interaction also improves the stability of the electroweak vacuum. In order to quantify these effects, we match the effective theory with NMSSM which minimally accommodates the singlino with MSSM. We derive matching conditions between effective theory parameters and those of NMSSM and perform two loop RGE analysis with one loop threshold corrections in order to check the viability of effective theory with respect to the various phenomenological constraints.
We find that viable pseudo-Dirac Higgsino DM puts an upper bound on the singlino mass scale, µ S < ∼ 10 8 GeV. The O(1) Yukawa coupling between Higgsino, singlino and Higgs makes the electroweak vacuum metastable or stable at all the scale up to M S for all the values of tan β allowed by perturbativity constraints. However, the observed Higgs mass together with constraints on the charged Higgs mass and Higgs to gauge boson couplings strongly disfavour tan β > ∼ 3 in almost all the cases discussed in this paper. It is shown that the viable pseudo-Dirac Higgsino DM, consistent with other phenomenological constraints, can be obtained within the underlying framework of GUT scale superymmetry for µ S ∈ [10 5 , 10 7 ] GeV, µ ≈ 1 TeV, tan β ∈ [1.2, 3] and with O(1) Higgsino-singlino-Higgs Yukawa coupling.
A dark matter particle χ of mass m χ , streaming with velocity distribution f (v) in our galaxy can interact with target nucleus X of mass m T . This interaction can be recorded by observing the recoil energy, E R , of the nucleus in the direct detection experiments, such as [61][62][63][64]. The differential rate of such events is given by [65] where µ T = m χ m T /(m χ + m T ) is the reduced mass of the nucleus-DM system, N T is number of target nuclei in the detector, ρ χ = 0.3 GeV/cm 3 is the local DM density in our galaxy.
v min is the minimum velocity of DM that can trigger a nuclear recoil of energy E R . v max is the maximum velocity of DM in the galaxy which is equal to v esc , i.e. the minimum velocity required for DM particle in order to escape from the galaxy. σ is the scattering cross section of DM particle with nucleus which is determined by the underlying particle physics framework. The integral in eq. (A1) is estimated in [65] using Maxwellian distribution for velocities. The obtained result depends on v obs in addition to v min and v esc where v obs takes into account for the relative motion between the earth and the rest frame of DM. We use the result of [65] for estimation of v obs . In the case of elastic scattering with a nucleus, i.e. χ+X → χ+X, the minimum velocity of DM particle required to produce nuclear recoil with energy E R is given by If χ is a pure Higgsino DM with m χ ∼ 1 TeV, the number of events estimated using the above formulae turns out to be much larger than the total observed events in the experiments such as Xenon 10 and Xenon 100. Therefore, the case of pure Dirac Higgsino DM is disfavoured.
In the case of pseudo-Dirac Higgsinos discussed in the paper, we have Majorana neutralinosχ 0 1 andχ 0 2 with mass difference ∆m 0 as given by eq. (22). In this caseχ 0 1 can scatter inelastically off the nucleus giving rise to a process:χ 0 1 + X →χ 0 2 + X. If ∆m 0 is small enough, the nuclear recoil can also be observed in direct detection experiments. The value of v min that can give rise to nuclear recoil with energy E R is then given by and the scattering cross section σ in eq. (A1) is replaced by inelastic cross section, σ Inelastic ≡ σ(χ 0 1 + X →χ 0 2 + X). In the framework under consideration, the dominant contribution to σ Inelastic arises from the exchange of Z bosons as the couplings betweenχ 0 1,2 and THDM scalars are suppressed by O(v/µ S ), see eq. (14). Using eq. (18) and kinetic term of Higgsinos in eq. (14), one obtain the following effective interaction between theχ 0 1,2 and quarks after integrating out the Z boson: where T q 3 (Q q ) is 1/2 (2/3) for q = u, c, t and −1/2 (−1/3) for q = d, s, b. G F is Fermi constant and θ W is weak mixing angle. We have sin 2θ ≈ 1 for pseudo-Dirac Higgsinos in our framework. The effective spin-independent inelastic cross section between the neutralinos and proton or neutron are estimated by summing over the valance quark contributions using the above results. This gives the following result for the cross section between the neutralinos and nucleus where N n (N p ) is number of neutrons (protons) in a given nucleus and F (E R ) is the nuclear form factor which parametrizes the distributions of protons and neutrons in the nucleus [66]. We estimate the total number of events R using eqs. (A1,A3,A5) for given range of recoil energy E R . The results are then compared with the data collected by the experiments: Xenon 10 with 58.6 live days of data and a target mass of 5.4 Kg [62], Xenon 100 with 224.6 live days of data and a target mass of 34 kg [63] , and Xenon 1T with 278.8 live days of data with target mass of 1300 Kg [64]. The values of parameters, mχ0 1 and ∆m 0 , are excluded with 90% confidence level if theoretically estimated number of events are greater than the observed number of events in accordance with poisson statistics. The results are displayed in Fig. 6. We find a conservative limit ∆m 0 ≥ 200 keV for mχ0 1 ≈ 1 TeV in order to evade the constraints from direct detection experiments on pseudo-Dirac Higgsino DM.

Elastic scattering
In this section we discuss the elastic scattering of pseudo-Dirac Higgsino DM in the underlying framework. The relavant terms can be read from the effective Lagrangian given in eq. (14) and are given by where g Z = g 2 Y + g 2 2 , g h = −c 1 sin α cos β + c 2 cos α sin β + d cos(α + β) , g H = c 1 cos α cos β + c 2 sin α sin β + d sin(α + β) .
The first term in eq. (A6) gives rise to spin-dependent (SD) contribution while the remaining terms induce spin-independent (SI) contributions to the elastic scattering. Upon integrating out the Z boson and THDM scalars, one obtains the following relevant effective operators: where b q h = −b q H = − sin α/ sin β for q = u, c, t and b q h = b q H = cos α/ cos β for q = d, s, b. The SD cross section is proportional to cos 2θ hence it is negligible in our framework.
Using the first in eq. (A7), we compute SI cross section of DM with a single proton which is given by and m p is mass of proton. The factors f T q are obtained from the lattice computation [67]. We estimate σ SI p for the values of parameters favoured by our results, i.e. tan β = 1.5, β − α = π/2, M h = 125 GeV, M H = 1 TeV and µ S = 10 5 GeV.We find σ SI p ≈ 10 −57 cm 2 which is much smaller than the sensitivity ( > ∼ 10 −44 cm 2 ) of current generation experiments [61]. Typically, the neutrino background, with cross section ∼ 10 −48 cm 2 , dominates over the DM signal [68]. Hence, it remains challenging to constraint the pseudo-Higgsino DM discussed in this framework from their elastic scattering signatures.

Appendix B: Renormalization group equations
We list 2-loop renormalization group equations for various couplings appearing in our framework which are obtained using publicly available package SARAH [44]. The couplings evolve according to the following equation: where C represents gauge, Yukawa or quartic couplings and Q is the renormalization scale. The step function Θ(Q − Q 0 ) = 1 for Q > Q 0 and it vanishes otherwise. The one and two-loop beta functions for the different couplings are as the following.