Supersymmetric proton decay revisited

Encouraged by the advent of a new generation of underground detectors – JUNO, DUNE and Hyper-Kamiokande – that are projected to improve significantly on the present sensitivities to various baryon decay modes, we revisit baryon decay in the minimal supersymmetric SU(5) GUT. We discuss the phenomenological uncertainties associated with hadronic matrix elements and the value of the strong coupling αs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s$$\end{document} – which are the most important – the weak mixing angle θW\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _W$$\end{document}, quark masses including one-loop renormalization effects, quark mixing and novel GUT phases that are not visible in electroweak interaction processes. We apply our analysis to a variety of CMSSM, super- and sub-GUT scenarios in which soft supersymmetry-breaking parameters are assumed to be universal at, above and below the GUT scale, respectively. In many cases, we find that the next generation of underground detectors should be able to probe models with sparticle masses that are O(10)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {O}}}(10)$$\end{document} TeV, beyond the reach of the LHC.


Introduction
The time is ripe to review and re-evaluate the prospects for observing proton decay in supersymmetric models.
On the one hand, we anticipate that a new generation of large underground neutrino detectors will be constructed in the coming years [1][2][3][4], with heightened sensitivity to a wide range of possible proton (and neutron) decay modes, including many suggested by minimal supersymmetric grand unified theory (GUT) models. As reviewed in [5,6] and discussed in more detail in Sect. 2, the JUNO experiment is expected to have 90% CL sensitivity to the p → K + ν decay a e-mail: jlevans@sjtu.edu.cn (corresponding author) mode at the level of a lifetime > 2 × 10 34 year after operating for 10 years -see Section 10 of [1], the DUNE experiment is expected to have 90% CL exclusion sensitivity to p → K + ν at the level of ∼ 6.5 × 10 34 yr after 20 years -see Chapter 4 of Volume 2 of the DUNE CDR [2,3], and the Hyper-Kamiokande (Hyper-K) experiment is expected to have 90% CL exclusion sensitivity to p → K + ν at the level of ∼ 5 × 10 34 yr after 20 years -see Section III.2 of the Hyper-Kamiokande TDR [4]. These sensitivities reach around an order of magnitude better than the current experimental limit [7,8].
On the other hand, searches at the LHC have pushed up the experimental lower limits on the masses of many sparticle species [9][10][11][12]. For example, searches for stronglyinteracting sparticles by the ATLAS and CMS Collaborations have excluded them in mass ranges to ∼ 2 TeV or more [13,14] when their data are interpreted in simplified supersymmetric models. In parallel, direct searches have improved upper limits on dark matter scattering cross sections [15][16][17], at least in models with simplifying assumptions about the input soft supersymmetry-breaking masses. These developments have two relevant implications: the absence of direct evidence of supersymmetry gives added motivation to indirect searches, e.g., via proton decay, and proton decay matrix elements are generically more suppressed if the sparticle spectrum is heavier.
When interpreting these broad changes in the experimental perspectives for detecting supersymmetric proton decay, there are several phenomenological considerations that need to be taken into account in calculating the proton decay rate in any specific supersymmetric scenario, as we discus in this paper. Although our analysis is in the context of one specific model, namely the minimal superysmmet-ric SU(5) GUT, we expect our considerations to have broad applications to other models. First of all, we provide a treatment of the uncertainties in the hadronic matrix elements of the relevant nucleon decay operators, for which we use updated lattice calculations [18]. Secondly, we update the GUT matching of couplings, combining the current best estimate of α s (M Z ) and the latest experimental estimate of sin 2 θ W . As we discuss in more detail below, the masses of the coloured GUT Higgses, H C , are particularly sensitive to the value of sin 2 θ W , and the current value leads to a smaller value of M H C than some previous estimates, accelerating proton decay. We also consider the impacts of uncertainties in quark masses, and stress the importance of including one-loop mass-renormalization effects, which have an effect that is not negligible. We also review the implications of the present uncertainties in quark-mixing effects. The conventional Cabibbo-Kobayashi-Maskawa (CKM) matrix elements have experimental uncertainties that are small enough to be relatively unimportant for our purposes. However, we recall that in even the minimal supersymmetric GUT there are two additional CP-violating phases that are unconstrained by electroweak-scale measurements [19]. These can have significant effects on the magnitudes of the nucleon decay matrix elements, possibly even altering the expected hierarchy of proton decay modes [20].
The layout of this paper is as follows. After our discussion of the present and prospective experimental sensitivities to various possible nucleon decay modes in Sect. 2, we present a review of the minimal supersymmetric SU(5) model in Sect. 3.1, followed by a discussion of the basic elements in our proton decay calculation in Sect. 3.2. In Sect. 4 we discuss various phenomenological aspects of proton decay calculations and the relative uncertainties they introduce into the calculation of the proton lifetime.
We begin Sect. 4 with discussions of the dominant uncertainties originating in the hadronic matrix elements that feature in different GUT models for nucleon decay [18], and in the experimental input value of α s , in Sects. 4.1 and 4.2 respectively. The latter strongly affects the mass of the colour-triplet Higgs, whose strong sensitivity to sin 2 θ W through the GUT matching conditions is discussed in Sect. 4.3. The sensitivities to the second-generation quark masses, m s and m c , are considered in Sect. 4.4, and we discuss the effects of their one-loop mass renormalizations in Sect. 4.5. Then, in Sect. 4.6 we review the status of quark mixing, and emphasize in Sect. 4.7 the important uncertainties associated with the additional CP-violating phases that appear even in the minimal SU(5) GUT. Possible implications of the non-unification of Yukawa couplings in this model are addressed in Sect. 4.8. We apply these considerations to various GUT models in Sect. 5, including the constrained minimal supersymmetric Standard Model (CMSSM) in which the soft supersymm-etry-breaking sfermion and gaugino masses are assumed to be universal at the GUT scale in Sect. 5.1 [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35], 1 super-GUT models in which universality is imposed at a scale higher than the GUT scale [20,[37][38][39][40][41][42] in Sect. 5.2, and sub-GUT models in which universality is imposed a lower scale [32][33][34]43,46] in Sect. 5.3. Finally, Sect. 6 summarizes our analysis and discusses our conclusions. The main part of our analysis concentrates on proton decay mediated by dimension-five operators, but we add for completeness an Appendix concerning the dimension-six operators that mediate p → π 0 e + .

Baryon decay: experimental status and prospects
The current lower limits on the lifetimes for many possible proton and neutron decay modes provided by the Super-Kamiokande experiment are summarized in Refs. [4,8,47,48]. As we review below, proton decay in the minimal supersymmetric SU(5) GUT is predominantly induced by the dimension-five effective operators [49,50] generated by color-triplet Higgs exchange. In this case, the most important decay mode is generally p → K + ν, for which the current lower decay lifetime limit is 6.6 × 10 33 yr [7,8]. 2 In the models we study, the decay n → K 0 ν is expected to occur at a similar rate, as these two decay channels are related through an SU(2) isospin rotation. However, the current lifetime limit on this decay mode is much weaker, τ (n → K 0 ν) > 2.6 × 10 32 year [48,52], and we expect the future sensitivity also to be much less than that for p → K + ν. As we discuss later, there are (small) regions of parameter space where the decay p → π + ν may occur at a rate comparable to p → K + ν, but the current sensitivity to this decay mode, τ ( p → π + ν) > 3.9 × 10 32 year [53], is also considerably less than that for p → K + ν. Although there have been no dedicated studies of the sensitivities of future detectors to this decay mode, we expect them also to be significantly weaker than those to p → K + ν, and hence uncompetitive.
In the second column in Table 1 we summarize the current 90% CL limits on nucleon decay modes in units of 10 33 year. For other decay modes, see Refs. [4,8,47,48].
There are several approved experiments with interesting prospects for detecting nucleon decay. Among them, the DUNE experiment [2,3] is expected to offer the best sensitivity to p → K +ν . The DUNE far detector is a liquid-argon time-projection chamber (LArTPC) with a 40 kt fiducial mass, located ∼ 1.5 km underground at the Sanford Underground Research Facility in South Dakota, USA. This type of detector is advantageous for identifying K + , compared to water Cherenkov detectors such as Super-Kamiokande, since a slow K + produced in the decay of a proton has a high ionization density. This information, as well as the reconstruction of the decay products associated with the K + , enables a LArTPC to identify the K + track with high efficiency and thus to have a good sensitivity to nucleon decay channels with a K + in the final state. The expected 90% CL limit on the p → K +ν channel in 10 (20) years is ∼ 3.3 (6.5) × 10 34 yrs [2,3], where the initial run with 10kt and adding another 3 There is a short discussion of p → e + π 0 in the "Appendix". 4 A recent update of this limit is reported in Ref. [51]: τ ( p → e + π 0 ) > 2.0 × 10 34 year. 5 A preliminary update, τ ( p → μ + π 0 ) > 1.2 × 10 34 year, is given in Ref. [51]. 10 kt each year for four years is assumed. The 3-σ discovery reach of DUNE is estimated in Ref. [4] based on the expected signal efficiency and background rates given in Refs. [2,3]: τ ( p → K +ν ) ∼ 3 (5) × 10 34 yrs in 10 (20) years.
Hyper-Kamiokande (Hyper-K) is a water Cherenkov detector with a fiducial mass of 187 kt [4]. The candidate site of the Hyper-Kamiokande detector is in the Tochibora mine in Gifu prefecture, Japan, which is located 8 km south of Super-Kamiokande. It is challenging to detect the K + from nucleon decays in a water Cherenkov detector since its momentum is below the Cherenkov light threshold, and thus only the decay products of K + can be used for the identification of K + . The expected 90% CL limit on the p → K +ν channel in 10 (20) years is ∼ 3.2 (5) × 10 34 years [4]; if a second tank is installed after 6 years, the reach becomes ∼ 4 (7) × 10 34 years in 10 (20) years. The 3-σ discovery reach is ∼ 2 (3) × 10 34 years in 10 (20) years. We note that Hyper-K has not yet reported prospective sensitivities for p, n → π +,0ν .
In the third and fourth columns in Table 1, we summarize the future 3-σ discovery and 90% CL limit sensitivities, respectively, for 10 (and, where available, 20) year operations of the future experiments. As can be seen, considerable improvements in sensitivities are expected for many decay channels, including p → K +ν in particular.
Since the coloured triplets in the H and H representations are required to be heavy (as they mediate proton decay), while the corresponding doublets need to be light (they correspond to the MSSM Higgs doublets), we must split the masses of the doublet and triplet, which requires a fine-tuning condition μ H − 3λV V . In this case, the color-triplet Higgs states have masses M H C = 5λV . In addition, the masses of the color and weak adjoint components of are equal to M = 5λ V /2, while the singlet component of acquires a mass The Yukawa couplings (h 10 ) i j and h 5 i j in Eq. (1) have redundant degrees of freedom, some of which can be eliminated by field re-definitions of i and i . The coupling (h 10 ) i j is a symmetric matrix, and thus has six independent complex components, while h 5 i j has nine. The field redefinitions form an U(3) ⊗ U(3) transformation group, and thus the number of physical degrees of freedom is 12+18−(9×2) = 12. Among them, six correspond to quark mass eigenvalues and four to the CKM matrix elements, so there are two additional GUT phases [19]. In this paper, we follow Ref. [89] by adopting the basis in which where h 10,i and h 5, j are the eigenvalues of (h 10 ) i j and h 5 i j , respectively. The phase factors ϕ i satisfy the condition and thus two of them are independent. We discuss the effect of these GUT phases on the nucleon decay rates in Sect. 4.7. The eigenvalues h 10,i and h 5,i are obtained from the quark Yukawa couplings at the GUT scale. We use the same matching conditions as used in Ref. [20], namely for the third-generation Yukawa couplings, 6 while for the first-and second-generations we use where and f e i (M GUT ) are the uptype, down-type, and charged lepton Yukawa couplings at the GUT scale M GUT , respectively, and V i j are the CKM matrix elements. We note that there is an ambiguity in the determination of h 5,i , and discuss the implications of this ambiguity in Sect. 4.8.
The MSSM matter superfields are embedded into the SU(5) matter multiplets as so Eq. (1) leads to: 7 where a, b, c denote the colour indices, and the dots indicate contractions of SU(2) L indices with the anti-symmetric tensor. We note that the new GUT phase factors (3) appear only in the couplings of the colour-triplet Higgs multiplets, as expected since they are unobservable at the electroweak scale.
In the super-GUT models we discuss below, the supersymmetry-breaking soft mass terms for the MSSM fields are determined at the GUT scale (defined as the energy scale when the electroweak gauge couplings are unified) by a set of matching conditions. The soft supersymmetrybreaking terms in the minimal supersymmetric SU(5) GUT are where ψ i and φ i are the scalar components of i and i , respectively, the λ A are the SU(5) gauginos, and we use the same symbols for the scalar components of the Higgs fields as for the corresponding superfields.
In the CMSSM and its generalizations considered here we impose universality conditions for the soft-mass parameters at a soft supersymmetry-breaking mass input scale M in : 7 To simplify the following expression, we replace ( f d3 + f e3 )/2 by f d3 in this equation, but we use the full expression in our calculations.
The bilinear soft SUSY-breaking therms B and B H are discussed below. When M in = M GUT , the above conditions are equivalent to those in the CMSSM and the renormalization group equations (RGEs) are run between the weak and the GUT scale. When M in < M GUT , as in sub-GUT models [43][44][45], the RGEs are run only up to M in . In contrast, in the super-GUT scenario where M in > M GUT [20,[37][38][39][40][41][42], RGE running occurs for MSSM parameters between the weak and GUT scale, where they are matched into SU(5) GUT parameters which are then run up to M in . Boundary conditions are set at both the electroweak scale (e.g., for the gauge and Yukawa couplings) and at M in (e.g., for the soft supersymmetry-breaking parameters). For a more complete discussion of the RGE evolution in super-GUT models, see Ref. [20].
In the super-GUT scenario with M in > M GUT , after the parameters listed in Eq. (11) are run down to M GUT they are matched onto the corresponding MSSM parameters. Of particular importance are the matching conditions for the gauge couplings and gaugino masses. We use the DR scheme for the gauge couplings, and include the one-loop threshold corrections due to the GUT-scale fields. We can then take linear combinations of the matching conditions to obtain the following convenient expressions [20,[89][90][91]: where g 1 , g 2 , g 3 , and g 5 are the U(1), SU(2), SU(3), and SU(5) gauge couplings, respectively, and Q is a renormalization scale taken in our analysis to be the unification scale: Q = M GUT . 8 The last terms in equations (12) and (14) represent the contribution from the dimension-five operator where W ≡ W A T A denotes the superfields corresponding to the field strengths of the SU(5) gauge vector bosons V ≡ V A T A . Since V /M P 10 −2 , these terms can be comparable to the one-loop threshold corrections, and thus should be taken into account when discussing gauge-coupling unification [94]. In what follows, we use the notation The matching conditions for the gaugino masses are given by [59,94,95]: Again we see that the contributions of the dimension-five operator can be comparable to those of the one-loop threshold corrections.
In the absence of the dimension-five operator (i.e., when = 0), Eqs. (12)(13)(14) provide three conditions on the masses of M H C , M , and M X , as well as g 5 . In addition, we can relate the three masses to the GUT Higgs vev V through the couplings λ, λ , and g 5 respectively. This gives us six constraints on seven quantities: the three masses, the three couplings, and V . Thus only one of the two GUT couplings, λ or λ can be chosen as a free parameter. On the other hand, if = 0, λ and λ can be chosen independently with the following condition on the dimension-five coupling: which can be obtained from Eq. (12) with g 1 (M GUT ) = g 2 (M GUT ). We note that this relation for can also be used in the CMSSM (in which M in = M GUT ) if λ and λ are specified, even if no running above the GUT scale is considered. As we will see, 'turning on' enables the coloured Higgs mass to be increased, and thus increases the proton lifetime. The remaining MSSM soft supersymmetry-breaking mass terms and trilinear couplings are related at M GUT by The MSSM μ and B terms are [96] with In the absence of a more elegant solution for the separation of the GUT and weak scales, we must tune |μ H − 3λV | to be O(M SUSY ). In practice, μ and B are determined at the electroweak scale by the minimization of the Higgs potential as in the CMSSM. These are then run up to the scale where Eqs. (22) and (23) are applied. Then, using = 0 (which is stable against radiative corrections as shown in Ref. [97]), we can solve for B H and B at the GUT scale, which are needed in the matching conditions for the gaugino masses.
As was pointed out in [20], Eqs. (23) and (24) have no solution unless we take A λ 8m 2 . For super-GUT theories with A 0 = 0, like in the focus-point case we consider below, this conditions is generally not satisfied. However, these types of models can satisfy this condition if the following Giudice-Masiero terms [98] are added to the Kähler potential: Although the above terms only give a small correction to the B terms, they give an important contribution to (23), since where m 3/2 is the gravitino mass which sets the scale for m 0 and m 1/2 . The additional contributions in the above equations make it trivial to satisfy (23) even with A 0 = 0, see [99] for more discussion.
To summarize, the well-studied CMSSM is characterized by four parameters and one sign: In sub-GUT models we must in addition specify the input universality scale: and in super-GUT models we also need to specify the Higgs couplings: which are fixed at Q = M GUT . As noted above, CMSSM models with non-zero values of also require fixing these GUT Higgs couplings.

Nucleon decay in minimal supersymmetric SU(5)
We now review the calculation of the proton decay lifetime in the minimal supersymmetric SU(5) model [33,59,73,89,[100][101][102][103][104]. As we have seen in Sect. 2, the forthcoming experiments are expected to offer great sensitivities to nucleon decay. To make the best of these experiments, therefore, it is desirable to formulate a precise and systematic method for the computation of nucleon decay lifetimes. To that end, we adopt the method of effective field theories. In this method, the fundamental theory is matched onto a low-energy effective theory at a high-energy scale (the GUT scale in our case), where the effect of heavy states (the GUT-scale fields) is included into the Wilson coefficients of higher-dimensional effective operators. We then run the Wilson coefficients down to the hadronic scale by using RGEs, which allows us to resum large logarithmic radiative corrections due to the large hierarchy in energy scales. The long-distance QCD effect is taken into account through the calculation of hadronic matrix elements. With this procedure, we can separate the short-and long-range contributions to the decay amplitude in a consistent manner. Note that this prescription is the same as those used for the calculation of precision physics observables, such as flavour observables [105] and the dark matternucleon scattering cross section [106,107]. As we discussed above, the most important decay mode is p → K +ν , which is induced by the exchange of the colourtriplet Higgs multiplets [49,50]. The effective Lagrangian for this contribution is where the effective operators O 5L i jkl and O 5R i jkl are defined by We note that antisymmetry with respect to the colour indices requires that the operators include at least two generations of quarks. For this reason, the dominant decay modes generally contain a strange quark in the final state, such as p → K +ν [108,109].
The Wilson coefficients C i jkl 5L and C i jkl 5R are run down to the supersymmetric scale M SUSY using the RGEs where Q denotes the renormalization scale. At the scale M SUSY , sfermions are integrated out via the wino-or Higgsino-exchange one-loop diagrams. The low-energy effective Lagrangian below the supersymmetric scale is then given by where the effective operators have the form with i = 1, 2, j = 2, 3, and k = 1, 2, 3, and their Wilson coefficients are evaluated as with From the supersymmetric breaking scale to the electroweak scale, we use the RGEs [110] μ d dμ where y u j denotes the SM up-type Yukawa couplings. Below the electroweak scale, the effective interactions that give rise to the p → K +ν k decay mode are described by where the coefficients of these interactions are obtained at the electroweak scale as We then use the two-loop RGE given in Ref. [111] to evolve these coefficients down to the hadronic scale μ had = 2 GeV, and finally obtain the partial decay width of the p → K +ν i mode as where m p and m K are the masses of proton and kaon, respectively. Since the experiments cannot determine the flavor of the neutrino, below, we will use the notation that p → K +ν represents the sum of the decays to all neutrino flavors. The decay amplitude A( p → K +ν i ) is the sum of the products of Wilson coefficients with hadronic matrix elements: The hadronic matrix elements are evaluated at the scale μ had = 2 GeV using QCD lattice simulations, which we discuss in detail in Sect. 4.1.
The dimension-five effective interactions in Eq. (32) also induce other nucleon decay modes, such as p → π + ν and n → π 0 ν. The calculation of the decay rates of these mode is the same as for p → K + ν above the electroweak scale. The effective interactions for these decay modes below the electroweak scale are where the matching conditions for the Wilson coefficients are Using these interactions, we compute the partial decay widths of p → π + ν and n → π 0 ν as where m n and m π are the masses of neutron and pion, respectively, and Again, since the experiments are unable to determine the flavor of the neutrino, we will use p → π +ν and n → π +ν to represent the sum of decays to all neutrino flavors. The dimension-six nucleon decay can also be computed in a similar way, which we show in the "Appendix", for completeness.

Uncertainties in nucleon decay calculations
In this section we discuss various uncertainties in the calculation of the proton decay rate for fixed values of the supersymmetric GUT model parameters. We start with the uncertainties in the hadronic matrix elements and the experimental input value of α s , which are generally the most important. 9 We also consider the effects of uncertainties in the weak mixing angle, quark masses, loop corrections to quark Yukawa couplings, and quark mixing parameters. The experimental inputs for α s and the quark masses that we use are listed in Table 2 [48]. 10 As many of the contributions to the proton lifetime uncertainty are small compared to those arising from the matrix elements and the value of α s , when we compute the "total" uncertainty in τ p ≡ 1/ i ( p → K +ν i ), we propagate only the effects of the matrix elements and the direct effect of α s on M H C , which is the most important source of sensitivity to the value of α s . 11

Hadronic matrix-element uncertainties
The hadronic matrix elements that determine directly the most relevant proton partial decay rates have been updated recently [18]. The improvement in the update provides a total accuracy of 10 to 15 % of the matrix elements. Table 3 lists the 9 The effect of the uncertainties in α s is greatly diminished when the operator (15) is considered, as we discuss below. 10 We note that there the PDG [48] lists an updated value of α s = 0.1179 ± 0.0010, a change that is a small fraction of the uncertainty and does not affect our results significantly. 11 For example, changes in the renormalization-group (RG) running below M Z shift the Yukawa couplings of the light quarks at the GUT Scale by about 2 percent, a little more for charm. The shifts in the Yukawa couplings and the Wilson coefficients from the RG running above M Z are estimated to be smaller, as are the effects on the soft masses.
values of the most relevant baryonic decay matrix elements calculated previously in [112] (second column) and recently in [18] (fourth column), including both the statistical and systematic uncertainties, which are indicated by (...)(...). Also shown in the third and fifth columns are the corresponding total errors after combining these uncertainties in quadrature. We see that in some cases the calculated matrix elements have changed significantly between [112] and [18], and that the uncertainties have been reduced substantially in every case. Both of these two simulations utilize the same gauge ensemble of N f = 2 + 1 domain-wall fermions with the same lattice spacing a = 0.11 fm, lattice volume (2.65) 3 fm 3 , and the pion mass 0.34-0.69 GeV, and directly compute the three-point function of the nucleon-to-pseudoscalar transition with an insertion of the baryon-number violating operators ("direct method"). 12 In Ref. [18], they use an algorithm called all-mode-averaging (AMA) [113][114][115], with which the statistical error has significantly been reduced. In addition, an error in the renormalization-scheme matching factors was corrected in Ref. [18]; this leads to a 6-7% change in the matrix elements (see Footnote 2 in Ref. [18]). Notice that these simulations are still performed with an unphysical pion mass. A simulation at the physical point is on-going [116], which is expected to reduce the systematic uncertainty associated with the chiral extrapolation. All in all, a precision with < 10% uncertainty is expected to be achieved in 5 years [117].
We focus in this paper primarily on the p → K +ν decay mode, which is determined by the four matrix elements that are featured in Fig. 1. This figure illustrates the sensitivities of the p → K +ν decay rate to variations of these four hadronic matrix elements, which are each given in units of the total uncertainties of the new elements given in column 5 of the Table. Thus, σ = 0 corresponds to the current central value of each of the four matrix elements and σ = ±1 corresponds to adding (subtracting) the 1-σ total uncertainty to (from) the corresponding hadronic matrix element. The previous values of these matrix elements would lie at σ = −0.83 which are significant changes, especially for the last two. These large changes in the matrix elements are responsible for the bulk of the reduction in the p lifetime relative to values found in previous work [20]. The sensitivities to the four matrix elements shown in  15) to be nonzero for the same stop-coannihilation point. In this case, we see that while the relative sensitivity to the matrix elements is similar, the lifetime is significantly increased. Finally, in the lower right panel we show a super-GUT example with M in = 10 17 GeV. Because the super-GUT running tends to reduce the stop mass relative to the other sfermion masses, the proton lifetime becomes much more sensitive to the Wilson coefficient arising from Higgsino exchange. This alters the sensitivity to the hadron matrix elements.
We illustrate in Fig. 2 the effect of the uncertainties in the hadronic matrix elements on the allowed ranges of the input parameters in a representative (m 1/2 , m 0 ) plane of the CMSSM. Here and in the remaining figures in this section, we consider the (m 1/2 , m 0 ) plane for the fixed values tan β = 5, A 0 /m 0 = 3 and μ > 0. The input universality scale is taken to be the GUT scale, defined as the renormal-ization scale for which g 1 = g 2 . There are in general two dark red shaded regions in each figure. The lower region where m 1/2 m 0 is excluded because there the lighter stau is the lightest supersymmetric particle (LSP), whereas the upper region with m 0 m 1/2 is excluded because there the lighter stop is the LSP. Along the boundary of the stop LSP region, there is a very thin blue strip where stop coannihilation [118][119][120][121][122][123][124][125][126][127][128][129] is effective in reducing the LSP relic density to match the cold dark matter density determined by Planck [130]. 13 We allow the relic density to vary between 0.01 < χ h 2 < 2.0 to enhance the visibility of the strip in the figures. The points chosen in the left panels of Fig. 1 lie on this strip at m 1/2 = 9.79 TeV. We note that regions of the (m 1/2 , m 0 ) plane that would correspond in a conventional cosmological scenario with adiabatic expansion to a cold dark matter density larger than that determined by Planck [130] would, however, be allowed in scenarios with late entropy generation. 14 The bulk regions of the displayed planes could therefore be allowed in such a case. The red dotdashed lines are contours of constant Higgs masses between m h = 122 and 130 GeV in increments of 1 GeV as calculated using FeynHiggs [132][133][134][135][136][137].
The solid black curve in Fig. 2 shows the contour of constant p → K +ν lifetime set at its current lower limit of 0.066 × 10 35 yrs for central values of the matrix elements 13 The region where the LSP accounts for cold dark matter appears as a thin strip here, largely because we display only a particular slice of the parameter space with fixed tan β and A 0 . In a more general exploration of the CMSSM parameter space, the region of the (m 1/2 , m 0 ) plane where the desired dark matter density can be attained is broader, see, e.g., [31]. 14 For a recent analysis in such a scenario, see [131].  [18] and other parameters. Its location is significantly higher than in previous work [20], due to the numerous updates incorporated here. These include updates to the hadronic matrix elements, the value of α s and the value of sin 2 θ W , as well as the one-loop correction to the charm quark Yukawa coupling. The dotted black curve shows the shift in the lower limit induced by a 1-σ tot increase in the p → K +ν lifetime where, as described above, σ tot includes only the contributions from the variations in the matrix elements, σ had , and the effect of the variation in α s on M H C , σ 2 τ p (discussed further in the following section) added in quadrature: Here σ had is the combined uncertainty due to all the hadronic matrix elements entering into τ p , which can be read from the amplitude of Eq. (44). The contour corresponding to a 1-σ decrease in the proton lifetime is off the scale of the plot. The pair of dot-dashed curves show the ±1σ uncertainty stemming from the hadronic matrix elements alone. That is, we are plotting the contours where τ p ± σ had = 0.066 × 10 35 years.  Table 3 (dot-dashed lines), and for a one-standard-deviation variations in both the hadronic matrix elements and α s (dotted). The red dot-dashed contours show the lightest Higgs mass as calculated using FeynHiggs [132][133][134][135][136][137]. In the upper red shaded region, the lighter stop is the LSP, and in the lower red shaded region, the lighter stau is the LSP. In the blue strip along the stop LSP region, the LSP has an enlarged relic density range, 0.01 < χ h 2 < 2.0, to enhance its visibility 4.2 Dependence on the strong coupling Despite the impressive reduction of the uncertainty in the experimental value of α s [48], the proton lifetime still varies drastically for a 1-σ change in α s . This can be understood from Eqs. (12)(13)(14). In particular, Eq. (12) indicates that the mass of the coloured Higgs field, M H C is exponentially sensitive to changes in α s = g 2 3 /4π [91]. 15 As the relevant Wilson coefficients are inversely proportional to the coloured Higgs mass, the proton decay rate is proportional to M −2 H C , and so is quite sensitive to variations in α s .
To estimate numerically the sensitivity of the proton decay width to variation in α s , we solve Eq. (12) for M H C and 15 However, as we discuss in more detail below, this sensitivity is substantially reduced when the dimension-five coupling c in Eq. (15) is allowed to be non-zero and vary while λ and λ are kept fixed. In the lower set of curves, the coupling associated with the dimension 5 operator is non-zero, and the meanings of the curves are similar. Also shown is the propagated uncertainty in α s alone, shown by the dotdashed curves, which are now clearly distinct from the dotted curves using σ tot where α s is the 1-σ uncertainty in α s and we have assumed that the variation of M H C is dominantly determined by g 3 . This then gives However, a variation in α s also leads to changes in the light quark masses of roughly 2% for a 1σ change in α s , which leads in turn to an 8% change in the proton decay rate. As indicated above, when we consider the 'total' proton lifetime uncertainty, the variation connected to the quark masses is not included, though it is included when we show the isolated effect using α s ± α s . For this reason, the 1-σ spread due to α s shown below appears larger than the 'total' uncertainty.
The upper set of curves in Fig. 3 have been produced assuming that the coupling of the dimension-five operator c vanishes. As in Fig. 2, the solid black curve in Fig. 3 shows the contour of constant p → K +ν lifetime set at its current lower limit of 0.066×10 35 yrs for central values of the matrix elements and model parameters. Similarly, the dotted black curve shows the shift in the lower limit induced by a 1-σ tot decrease in the proton decay rate.
We also show as dashed black lines in Fig. 3 the result of varying α s within its uncertainty ±0.0011 when determining all the supersymmetric spectrum and other observables, including τ ( p → K +ν ). The lower black dashed line in Fig. 3 corresponds to τ ( p → K +ν ) = 0.066 × 10 35 yrs when calculated with α s = 0.1192, whereas the upper black dashed line corresponds to calculations with α s = 0.117, i.e., ±1 − σ excursions in α s . As noted above, the strong dependence on α s can be understood largely from the fact that varying α s affects M H C , but also from the variations in the values of the light quark masses when evolved both between M Z and M GUT and between 2 GeV or m c and M Z .
In particular, there are two aspects that affect the evolution between M Z and M GUT . First, the running of the light quark masses is completely controlled by α s due to the relatively small values of the couplings f d , f s and f c : where g = g 1 √ 3/5. This contrasts with the running of the b and t Yukawa couplings, e.g., Secondly, the one-loop corrections, m q , to light-quark masses due to gluino loops are also controlled by α s , e.g., for m s [139]: where B 0 and B 1 are the Passarino-Veltman functions given in Ref. [139], and θ m s represents the mixing angle between the strange squarkss 1 ands 2 with mass eigenvalues of ms 1 and ms 2 , respectively, which is determined by . Therefore, for larger α s there is a bigger change in m q , and hence a bigger change in the value of the decay amplitude for the proton lifetime.
We note that the sensitivity to α s is more pronounced for μ < 0 than for μ > 0, simply because the angle θ m s above changes sign due to the change of sign in μ, giving a bigger contribution to m s for μ < 0. For example, in the region of m 0 ∼ 10 TeV and M 1/2 ∼ 10 TeV, the contribution to m s from gluinos is of the order of 26% for μ negative, while for μ positive it is of order 13%.
As for the values of the light quark masses at M Z , it is well known that α s increases considerably between M Z and m c . Hence a 1-σ change in α s affects by about 2% the estimation of the light quark masses at M Z that we use as inputs to our calculations at M Z , as shown in Table 4.
In contrast to the above analysis, we show in the lower set of curves in Fig. 3 the contours of τ ( p → K +ν ) = 0.066 × 10 35 yrs when the dimension-five operator coupling, c, is allowed to be non-zero with λ and λ fixed, e.g., here we set λ = 0.6 and λ = 0.0001. In this case, Eq. (13) fixes the combination M 2 X M that in turn fixes the the coloured Higgs mass: Then Eq. (53) becomes Thus we expect the uncertainty in the proton lifetime to be significantly less sensitive to the uncertainty in α s , by a factor ∼ 1/15. This is seen in the lower set of curves in Fig. 3. Since M H C is now essentially fixed, the proton lifetime is substantially larger, and the limit contour (shown again as the solid curve) appears at lower m 1/2 and m 0 . On either side of the central curve, we show 3 sets of curves displaying the uncertainty due to α s . Nearest the centre, the dot-dashed curves correspond to the propagated variation due to α s alone. In the upper set of curves, we did not show this, as it would have been indistinguishable from the dotted curve showing the total sensitivity, which was dominated by α s . Here we see clearly that, with c = 0, the uncertainty due to α s is greatly diminished. The dashed curves show the shift in the limit contour when the values α s = 0.1170 and 0.1192 are used for the supersymmetric spectrum and all other observables, as was done in computing the dashed curves in the upper part of the figure with c = 0. Finally, the dotted curves show the total propagated uncertainty, which is now dominated by the uncertainty in the matrix elements.

Dependence on the weak mixing angle
It was assumed in previous work [20] that sin 2 θ W = 0.2325 in the MS prescription, and this value was taken as an input condition at M Z . However, the precision of measurements of electroweak symmetry breaking (EWSB) observables warrants paying careful attention to the precise input value of sin 2 θ W , and in our calculations here we specify sin 2 θ W in the DR scheme, which can be extracted from where the quantities M DR (W,Z ),susy (M Z ) contain one-loop corrections to the W and Z boson masses calculated in the DR scheme. Alternatively, sin 2 θ W | DR can be extracted from where the current measurement of sin 2 θ W,eff (M Z ) is 0.23155(4) [48], and where Z γ ( p 2 ) is the mixed self-energy of Z and γ at the momentum scale p 2 , and V is a function of the DR quantities,ŝ ≡ sin θ DR W,susy (M Z ) andĉ 2 = 1 −ŝ 2 [138,139]. In the expression Eq. (63) above, as a first approximation c on the right-hand side of the equation can be taken as c = cos θ W,eff (M Z ), rather thanĉ = cos θ DR W,susy (M Z ). We choose to use sin 2 θ W,eff as our starting-point, since the MSSM corrections to sin 2 θ W,eff (M Z ) have been studied in some depth. In particular, for supersymmetric masses bigger than 1 TeV these corrections are known to be O(10 −4 ) and always negative [139].
The factorκ in Eq. (63) can be interpreted as the conversion factor to the DR scheme, where typically 1/Reκ represents a decrease by another amount of O(10 −4 ). Therefore we expect sin 2 θ DR W,susy (M Z ) to vary in the range (0.2312, 0.2315), depending on the supersymmetric contribution to sin 2 θ W,eff (M Z ). Since the proton lifetime is relatively insensitive to variations of O(10 −4 ) in sin 2 θ W | DR , as we see below, we consider a precise computation for each point of the parameter space to lie beyond the scope of this work. We illustrate in Fig. 4 the sensitivity of the p → K +ν lifetime calculation to varying sin 2 θ W | DR over the range 0.2312 to 0.2315, corresponding to the upper and lower black dashed lines. As in Fig. 3, the solid black contour shows the position of the limit τ ( p → K +ν ) = 0.066 × 10 35 yrs for central values of the inputs, and the dotted black contour shows the shift in this limit due to the σ tot uncertainty. We see that the induced uncertainty associated with sin 2 θ W is significantly smaller than that due to the hadronic matrix elements and α s . On the other hand, using the previous value of 0.2325 would have given quite different results, as illustrated by the solid brown curve in Fig. 4.

Sensitivities to quark masses
The amplitudes A( p → K +ν i ) for i = e, μ, τ are the following sums of products of the Wilson coefficients with hadronic matrix elements: Wino exchange contributes to the Wilson coefficients C L L (usdν i ) and C L L (udsν i ), which may be approximated by where m d i are the masses of the down-type quarks. On the other hand, as can be seen from Eq. (42), Higgsino exchange contributes only to A( p → K +ν τ ), via C RL (usdν τ ) and C RL (udsν τ ), which are given approximately by We find that the total decay width ( p → K + ν) = i=e,μ,τ ( p → K + ν i ) is dominated throughout the plane by the contributions as a result of the dependences on quark masses, CKM elements and phases that we describe in this and the following sections.
We first discuss the sensitivity to m s in the range m s = 93 +11 −5 MeV when tan β = 5 and A 0 /m 0 = 3, assuming GUT-scale universality. When the GUT phases are zero, all the contributions in Eq. (64) are of the same order and the contribution of the second term in C L L (usdν 2 ) is maximized, see Eq. (65), rendering this contribution of the same size or, in most of the parameter space, even larger than that proportional to C RL (udsν τ ).
Since C L L (usdν 2 ) is proportional to m s and the uncertainty in m s is between −5% and +12%, any change in m s affects the p → K + ν lifetime more than the other quark masses. We show in Fig. 5  The second most important quark-mass sensitivity is that to m c , which contributes to the second term in Eq. (65). We vary m c in the range m c = 1.27 ± 0.02 GeV, and find that, when tan β = 5, A 0 /m 0 = 3 and M in = M GUT the sensitivity to m c is less than that due to the uncertainty in m s as seen in Fig. 6.
The sensitivity to the masses m d , m b and m t is very mild. We can see from Eq. (67) that C RL (udsν τ ) ∝ m d , but the contribution to τ ( p → K + ν) from this Wilson coefficient is suppressed when the GUT phases are chosen so as to maximize the contribution to A( p → K +ν μ ) proportional to C L L (usdν μ ) K + |(us) L d L | p . Given the precision in the measurement of the top mass, m t = 172.9 ± 0.4 GeV [48], its effect on the p-lifetime is negligible, even though all the dominant Wilson coefficients, Eqs. (65) and (67), are proportional to f t . The sensitivity to m b = 4.18 +0.03 −0.02 GeV is also negligible, because it does not enter in either of the leading contributions, C RL (udsν τ ) K + |(us) R d L | p and C L L (udsν μ ) K + |(ud) L s L | p , to the total decay amplitude.

Sensitivities to one-loop mass renormalization effects
In Sect. 4.2 we detailed how α s enters, and controls, the 1loop corrections. In Table 5 we illustrate the effects of varying α s in the 1-loop corrections to the quark masses with exam- In the previous Section, we discussed how m d , m s , and m c enter the total decay width, as they contribute to Eqs. (65)(66)(67). In Figs. 5 and 6 we show the effects of the 1-loop corrections to m s and m c . In Fig. 5, the dot-dashed line shows the position of the lifetime limit when the 1-loop corrections to m s are ignored. As one can see, the curve lies above the nominal central contour (where the correction is included), indicating that the correction to m s increases the lifetime and weakens the limit (allowing lower sparticle masses). The effect of the correction to m d is qualitatively similar but less important and is not shown. In contrast, the dot-dashed line in Fig. 6 shows the p-lifetime limit calculating m c (M Z ) without loop corrections. In this case, we see that the one-loop correction significantly decreases the proton lifetime, making the constraint stronger so that the central limit lies at higher sparticle masses. This effect has a bigger impact than the 1.6% variation due to the uncertainty in m c .

Quark mixing uncertainties
The minimal SU(5) GUT does not contain a way to describe fermion mixing, but we know that any additional part of the theory which can describe it must reproduce at low energy the CKM elements within their experimental error. We there-  fore explore the sensitivity to this uncertainty through the fitted values of the Wolfenstein parameterisation of the CKM matrix [48].
The CKM phase δ plays no role, so the 3 relevant parameters are A, ρ and λ. Of these, by far the greatest sensitivity is to A, as we illustrate in Fig. 7 assuming A = 0.836 ± 0.015, tan β = 5, A 0 /m 0 = 3, M in = M GUT and μ > 0.
To understand the sensitivity to the uncertainty on A, we see from Eq. (65) that C L L (udsν m u) ∝ V * us V td V ts . 16 In terms of the Wolfenstein parametrization this can be written as C L L (udsν m u) ∝ A 2 λ 7 . We can then see from Eq. (67) that 16 The second term in C L L (udsν μ ) is in principle proportional to V * us V cd V cs . However, this term is proportional to the phase factor e i(φ2−φ3) , and we scan over phases so as to minimize the rate, bringing this term as close as possible to −1, thus concealing the dependence on V cd V cs .
Hence the dependence on the quark mixing matrix reduces to those on the A and λ parameters, which have uncertainties of 1.8% and 0.2%, respectively. It is not a surprise, then, that the sensitivity to A is comparable to that of that of m c , whose uncertainty is 1.6%, see Fig. 6.

GUT phases
We now discuss the uncertainties associated with the GUT phases (3). Since the two terms in Eq. (65) have comparable magnitudes, the Wilson coefficients C L L (usdν i ) and C L L (udsν i ) may be suppressed in certain ranges of the GUT phases. A general overview of the dependence of the lifetime for p → π +ν in the plane of the two GUT phases (φ 2 , φ 3 ) for the CMSSM parameter choices tan β = 5, A 0 /m 0 = 3, m 0 = 14.13 TeV and m 1/2 = 9.79 TeV is shown in Fig. 8. The maximum value of the lifetime is indicated by a green triangle.
We mentioned in Sect. 4.4 that Higgsino exchange contributes only to A( p → K +ν τ ), via C RL (usdν τ ) and C RL (udsν τ ). These coefficients are approximately given by Eqs. (66) and (67), respectively, where we see that, unlike the coefficients in Eq. (65), their absolute values do not change when the phases vary. However, the difference in the phase structure from that in Eq. (65) contributes to the GUT phase dependence of A( p → K +ν τ ), which is different from that of A( p → K +ν e,μ ). This feature is seen in the left panel of Fig. 9, where we choose tan β = 5, A 0 /m 0 = 3, m 0 = 15.75 TeV and m 1/2 = 11 TeV as in Fig. 8, and φ 3 is chosen to maximize approximately the p → K +ν lifetime. We see that the ratio between the rates for p → K +ν e,μ (green dotdashed line and blue dotted line, respectively) is independent of the GUT phase φ 2 , whereas the rate for p → K +ν τ (red dashed line) has a quite different dependence on φ 2 . The solid black line is the total p → K +ν decay rate.
Another potentially important proton decay mode is p → π +ν i , whose decay amplitude is where C RL (uddν i ) is non-vanishing only for i = 3. The neutron decay mode n → π 0ν is also potentially important, and is given by the same Wilson operators: As seen in the right panel of Fig. 9, the rate for p → π +ν (blue dotted line) is smaller than that for p → K +ν for all values of φ 2 , though it becomes comparable for φ 2 ∼ 200 o . 17 The rate for n → π 0ν (red dashed line) is always smaller than that for p → π +ν for the central values of the hadronic matrix elements that we use. As mentioned earlier, we do not consider the experimental searches for p → π +ν and n → π 0ν , as the current limits on these decays are significantly weaker than those for p → K +ν , and no detailed studies are yet available for the next-generation detectors.

Yukawa non-unification
As can be seen in Eq. (9), in minimal SU(5) the lepton and down-type quark Yukawa couplings should be equal at the GUT scale. However, when running the physical values up from the electroweak scale, one finds that they are quite different for the first two generations, whereas Yukawa unification is a good approximation for the b and τ . The differences for the lighter generations can, however, easily be compensated by effects from physics above the GUT scale.
In particular, operators of higher mass dimension induced at the Planck scale that contribute to the Yukawa couplings may account for this difference [140][141][142][143][144]. 18 Among such Fig. 9 Left panel: sensitivity of the p → K +ν partial lifetimes to the GUT phase φ 2 for the same point considered in Fig. 8. Right panel: comparison of partial lifetimes for p → K +ν (black solid line), p → π +ν (blue dotted line) and n → π 0ν (red dashed line) as functions of φ 2 operators, those of lowest dimension are which yield Yukawa interaction terms when acquires a vev.
In particular, the first operator in Eq. (71) splits the lepton and down-type quark Yukawa couplings by the product of the superpotential coupling with V /M P ∼ 10 −2 , which is sufficient to explain the differences in the lepton and downtype quark Yukawa couplings for all of the three generations. The operators in Eq. (71) also modify the couplings of the colour-triplet Higgs fields to the quark and lepton fields, and thus directly affect the proton decay amplitude. Our ignorance of the coefficients c h in Eq. (71) leads to ambiguity in these couplings, which then results in the uncertainty in the Wilson coefficients in Eq. (34).
The range of this uncertainty is indicated by the differences between the quark and lepton Yukawa couplings. In the previous Sections we have chosen the quark Yukawa couplings, i.e., f s and f d . Since f s < f μ , we would expect that in general using the strange-quark Yukawa coupling may yield a longer lifetime than using the muon coupling whereas, since f d > f e , using the down-quark Yukawa coupling may give a shorter lifetime than using the electron coupling. These expectations are borne out in tests we have made using a CMSSM GUT point with tan β = 5, A 0 /m 0 = 3, m 1/2 = 9.8 TeV, m 0 = 14.1 TeV and μ > 0. Our default choice of Yukawa couplings, f s,d , yields a proton lifetime 5.4 × 10 33 y, whereas using f s,e yields a lifetime 5.6×10 33 y, a 4% difference. On the other hand, replacing f s by f μ yields a lifetime that is 23 times smaller. Thus, our choice f s,d is quite conservative, and the most conservative choice f s,e would have resulted in an insignificant difference.
In principle, couplings of the type (71) could also modify the pattern of quark mixing in GUT Higgs triplet interactions. However, we would not expect this to modify the generic prediction that the dominant proton decay mode should be into K + ν, which results from the combination of colour and flavour antisymmetry in the effective dimension-five interaction [108,109]. Nevertheless, it is clear that more detailed studies of this ambiguity in specific GUT models are warranted, though they lie beyond the scope of this paper.

Results
In this section we display (m 1/2 , m 0 ) planes for various choices of the supersymmetric model parameters. For CMSSM models with GUT scale universality, we show two sets of proton decay limit contours. Those in black are for the minimal supersymmetric SU(5) GUT, and those in green are calculated assuming that the dimension-five operator in Eq. (15) is present with c = 0. In both cases, the solid lines correspond to the proton decay lifetime limit of 0.066 × 10 35 yrs using the Standard Model inputs described in the previous Section. We also show dot-dashed lines corresponding a lifetime τ ( p → K +ν ) = 5 × 10 34 year, corresponding to the estimated 3-σ discovery sensitivity of the DUNE experiment after 20 years of operation (see Table 1). The dashed contours surrounding the solid contour correspond to the 1σ tot uncertainty in the position of the limit. As in the previous section, σ tot takes into account the propagated uncertainties from the hadronic matrix elements and the strong coupling as it affects M H C . For super-GUT CMSSM models, 1-σ variations. The green lines are corresponding results including the dimension-5 contribution discussed in the text. When present, the green dot-dashed curves correspond to the DUNE discovery sensitivity. The red lines are the indicated contours of m h the dimension-five operator is needed to satisfy the boundary conditions, and only one set of contours are shown and coloured black. At each point in the supersymmetric space, we choose the unknown GUT phases so as to minimize the p → K +ν decay rate.
In each (m 1/2 , m 0 ) plane, we show contours of m h calculated using FeynHiggs 2.14.1 [137] that are consistent with the measured Higgs mass within the estimated calculational uncertainties. These are shown as red dot-dashed contours. Regions of the planes that are shaded brick red are excluded because there the LSP would be charged. Typically, in such regions at large m 0 , the LSP is the lighter stop, and when present, brick red regions at lower m 0 contain a stau LSP. Regions shaded pink are excluded because there is no consistent electroweak symmetry-breaking vacuum. In addition, there are very narrow strips shaded blue where the LSP density calculated in standard Big Bang cosmology falls within the range allowed by Planck and other measurements. Here, to make these good relic density regions visible on the scale plotted, we allow the relic density to vary between 0.01 < χ h 2 < 2.0. In other regions of the (m 1/2 , m 0 ) planes the LSP would generally be overdense in the absence of some scenario for modified cosmological evolution with entropy generation (see, e.g., Ref. [131]). Table 6 Lifetimes at points with χ h 2 ≈ 0.12, and m h ≈ 125 GeV. Masses are in TeV and lifetimes in units of 10 33 years

The CMSSM
We begin the discussion of our main results with the CMSSM. We recall from Eq. (29) that the CMSSM is defined by four parameters given at the GUT scale, defined to be where the two electroweak gauge couplings are equal. Because we are primarily interested in calculating proton decay rates, we need to determine the mass of the coloured Higgs triplet, M H C , which we obtain from the matching conditions in Eqs. (12)(13)(14). As discussed earlier, in the CMSSM with M in = M GUT , we do not run the RGEs above the GUT scale, and no additional matching to GUT scale parameters is needed. As a result, we can define CMSSM models with the dimension-five operator turned off, i.e., c = 0. When this operator is turned on, fixing M H C requires specifying the SU(5) Higgs couplings λ and λ . In all figures below with A 0 = 0, we have fixed λ = 0.6 and λ = 0.0001. For A 0 = 0, we take λ = 0.1 instead, since otherwise the focus-point region would be pushed to very large values of m 0 where the RGE running becomes unstable. For more on the dependence of τ p on these two GUT couplings, see [20]. We show in Fig. 10 four examples of CMSSM planes. In the two left panels, we assume tan β = 5 with A 0 /m 0 = 3, whereas in the right panels we take A 0 /m 0 = −4.2 for the same value of tan β. In the upper two panels we take μ > 0, whereas μ < 0 in the lower panels. These values are chosen so as to bring the relic density strip (shaded blue) in a position to intersect with experimentally viable values of the Higgs mass (allowing for uncertainties in the theoretical calculation of the Higgs mass).
The upper left panel of Fig. 10 corresponds to the same choice of parameters as used in the previous section. Indeed, this panel is essentially a simplified version of that shown in Fig. 3, keeping only the central contour limits (for both c = 0 (black) and c = 0 (green)) along with the variation of these contours by ±1σ tot . (We recall that the stronger limit lies off the scale shown in the plot when c = 0.) Here and in all the other panels, the limit on the parameter space is greatly weakened when c = 0, as the coloured Higgs mass is much larger, being determined by Eq. (59) with our choices of λ and λ . In this panel, as in subsequent panels, we also show the location of the DUNE sensitivity τ p = 5 × 10 34 yrs by the dot-dashed curve with c = 0. When c = 0, the contour often lies beyond the range shown. It is found in the upper right corner of the panels of Fig. 10, shown by the dot-dashed green curves, except for the case of A 0 /m 0 = 3 and μ < 0 (lower left panel), where it is outside the parameter ranges shown.
As an example, consider a point near (m 1/2 , m 0 ) = (9.8,14.1) TeV. It corresponds to a bino LSP with mass of roughly 4.8 TeV that is nearly degenerate with the lighter stop. The Higgs mass is close to 125 GeV. With c = 0, this point has a lifetime which is slightly less than the experimental limit, but is well within one σ of the limit. However, when c = 0 it lies safely above the lifetime limit. More specifically, with c = 0, we find τ p = 5.4 ± 4.6 × 10 33 years, whereas for the same choice of parameters when c = 0, we find = 8cV /M P = .0024, and τ p = 3.3 ± 0.6 × 10 34 years, a factor of over 6 times larger. This and other examples discussed in this section are summarized in Table 6.
When μ < 0 as in the lower left panel of Fig. 10, we see that the most important change is in the Higgs mass, which now requires significantly lower values of (m 1/2 , m 0 ) to obtain m h = 125 GeV with the correct relic density. In this case, unless c = 0, the proton decay lifetime limit is badly violated. The DUNE sensitivity lies beyond the range shown, implying that this experiment should be able to explore fully this parameter range. In the right panels of Fig. 10 where A 0 /m 0 = −4.2, the 125 GeV Higgs mass contours intersect the relic density strip at intermediate values of (m 1/2 , m 0 ). The bands surrounding these lines represent ±1σ tot uncertainties. The horizontal black lines are the current limit (solid) and the expected future (dot-dashed) 90% CL sensitivity for p → K + ν from DUNE [2,3]. The dot-dashed red lines are contours of m h evaluated using FeynHiggs 2.14.1 [137] (legends on the right axes), and the horizontal shaded band shows where the m h calculation agrees with the experimental measurement within the estimated uncertainty of ±3 GeV However, both examples require c = 0 to be compatible with proton lifetime limit. In Fig. 11, we show the proton lifetime calculated for c = 0 (black) and c = 0 (green) as functions of m 1/2 along selected stop coannihilation strips in the CMSSM and super-GUT models, with the corresponding scales on the left axes. The shaded bands surrounding the curves show the ±1σ tot uncertainties in our calculations. The red dot-dashed curves show the Higgs mass calculated with FeynHiggs 2.14.1 [137], with the corresponding scales on the right axes. The horizontal shaded region shows the estimated ±3 GeV theoretical uncertainty in the calculated Higgs mass. The two horizontal lines show the current limit on the proton lifetime (solid) and expected 20-year DUNE [2,3] sensitivity limit (dot-dashed).
The upper left panel of Fig. 11 is computed in the CMSSM using tan β = 5, A 0 /m 0 = 3, M in = M GUT , and μ > 0 as in the upper left panel of Fig. 10. In this case, we see that the current proton lifetime limit already excludes m 1/2 8 TeV when c = 0, but there is a portion of the parameter space extending to m 1/2 15 TeV that is also consistent with the experimental value of the Higgs mass. When c = 0, the current proton lifetime limit allows the range of m 1/2 5 TeV where m h > 122 GeV remains allowed. We also see that DUNE [2,3] should be able to explore the entire m 1/2 range shown. The upper right panel assumes the same CMSSM input parameters, but with μ < 0 as in the lower left panel of Fig. 10. In this case, when c = 0 the calculated Higgs mass exceeds 131 GeV when the proton lifetime is sufficiently long. In contrast, with c = 0, as in the case of μ > 0, a range of m 1/2 3 TeV and 7 TeV is compatible with the measurement of m h as well as the current proton lifetime limit. Here too, DUNE [2,3] should be able to explore the entire range of m 1/2 allowed by m h . A different region of the CMSSM parameter space where the relic density is acceptable is found when A 0 /m 0 = 0 and m 0 is large, namely the focus-point region where μ → 0 [146][147][148][149]. The value of m 0 at a focus point is sensitive to λ. For λ = 0.6, the value of m 0 needed to drive μ close to zero is so large that the RGE running becomes unstable. Therefore, in our analysis of the focus-point region we take λ = 0.1 and λ = 0.0001. Two examples of (m 1/2 , m 0 ) planes exhibiting the focus-point region are shown in Fig. 12, where we have chosen tan β = 3.25 (left panel) and tan β = 4.0 (right panel). In both of these examples, the LSP is Higgsino-like along the focus-point strip, and its mass is therefore close to 1.1 TeV everywhere along the strip [150,151]. The strip in these panels appears relatively thick because we have (as in previous plots) shaded the region where 0.01 < χ h 2 < 2. We note also the appearance of a red shaded strip below the focus point, where the chargino is the LSP.
In contrast to the previous examples, the proton lifetime constraints are not monotonic. Consider for example, the left panel of Fig. 12 with c = 0. We see two solid black contours corresponding to the proton lifetime limit of 0.066 × 10 35 years. One of the two spans the figure at m 0 between 10 and 20 TeV. The 2nd contour is found inside the blue shaded region, just below the boundary where there is no radiative EWSB. Between the two, the proton lifetime is found to be greater than the limit. Also within the blue shaded region, there is the 1-σ limit on the lifetime contour (dotted black). Above this dotted line, the proton lifetime is too short. The unusual suppression of the proton lifetime near the focus point is due to the reduction in M H C as μ decreases [104], which enhances the decay rate along the focus-point strip. In addition, the cancellation between the Higgsino-and Winomediated pieces (important at values of m 0 between the solid black lines) disappears due to the suppression of the Higgsino contribution, leading to a shorter proton lifetime. These two effects combined give a proton lifetime that is too short along the focus-point strip. (See, for example, the entry in Table 6 for this figure.) Another interesting feature in this Figure is the fact that the c = 0 proton lifetime constraint is sometimes stronger than that for c = 0. This is because we have fixed the values of λ and λ , which effectively fix the coloured Higgs mass. If this value is smaller than the value determined by Eq. (59) for c = 0, then in Eq. (16) suppresses the proton lifetime. In addition, for the regions with smaller m 1/2 there is an enhancement in M H C for c = 0 relative to the case with c = 0, so the corresponding constraint (shown in green) is weaker. Additionally, in the upper left corner of Fig. 12, both the Wino and Higgsino mass go to zero, further suppressing the proton decay width, whereas the opposite is the case for Outside the loop the lifetime is smaller than the expected DUNE reach. In the right panel of Fig. 12, which has larger tan β, nearly the entire region displayed is excluded by the proton lifetime constraint. The weaker constraint in this case is for c = 0, for which some region of parameter space is allowed if we consider a 1 − σ tot variation in the lifetime. Table 6 gives details for a point corresponding to each panel where m h = 125 GeV and χ h 2 = 0.12.

Super-GUT models
We next consider super-GUT models in which M in > M GUT , using for illustration M in = 10 17 GeV. Since the RGEs must now be run above the GUT scale, we must apply all of the boundary conditions discussed in Section 3.1. In particular, since we fix λ = 0.6 and λ = 0.0001 as in the previous section, we must take c = 0 in order to satisfy simultaneously Eqs. (12)(13)(14) and Eq. (59). The gaugino mass, M 5 , scalar masses, and A-terms are fixed by Eq. (11), which are matched to MSSM parameters at the GUT scale using Eqs. (17)(18)(19) and Eq. (21). As in the CMSSM, we do not specify either of the GUT B-terms at M in . Instead, again as in the CMSSM, B (and μ) are determined by the minimization of the Higgs potential at the weak scale. These are run up to the GUT scale and B H , B (and μ H ) are given by Eqs. (22) and (9) with = 0. We recall that μ = λ V /4, and note that B H and B are needed in the gaugino matching conditions at M GUT .
For better comparison with the CMSSM results in the previous Section, we choose the same input values of tan β and A 0 , and our illustrative super-GUT (m 1/2 , m 0 ) planes are shown in Fig. 13. In each of the planes, we see pink shaded regions where the minimization of the Higgs potential fails to provide a solution for μ. In the upper left panel, where tan β = 5, A 0 /m 0 = 3 and μ > 0, we see that the proton lifetime sensitivity is relatively weak and in most of the stop coannihilation strip (shaded blue) the proton lifetime is sufficiently long. DUNE will be able to explore much of this strip, extending to m 0 ∼ 17 TeV if proton decay is not seen. However, for this choice of parameters, the 125 GeV Higgs mass contour does not intersect the stop coannihilation region, which extends beyond the range shown. For this reason, no entry is given in Table 6 for this case. However, we note that the 123 GeV Higgs mass contour, which is acceptable given the uncertainties in the calculation of the Higgs mass, does intersect the stop coannihilation strip.
In contrast, the 125 GeV Higgs contours do intersect the relic density strips in all the three other panels of Fig. 13, and the locations of these intersections are summarized in Table 6. In the upper right panel, with A 0 /m 0 = −4.2 and μ > 0, the lifetime is well beyond the current limit, but well within the reach of DUNE. In the lower two panels, with μ < 0, the current limit on the proton lifetime intersects the stop coannihilation strip for values of the Higgs mass that are consistent with experiment, with the calculational uncertainties.
The proton lifetime profiles as a function of m 1/2 along the stop coannihilation strip for the two super-GUT models with tan β = 5 and A 0 /m 0 = 3 are shown in the two lower panels of Fig. 11. As c = 0 is necessary to satisfy the boundary conditions, only one lifetime profile is shown in each panel. When μ > 0, as noted earlier, the Higgs mass is low for the range of m 1/2 shown, whereas m h is consistent with experiment for a range of m 1/2 2 TeV and 7 TeV for μ < 0. Much of this range is compatible with the current proton lifetime limit. For both signs of μ, DUNE [2,3] should be able to explore most of the range of m 1/2 15 TeV.
We also show in Fig. 14  the Higgsino and Wino masses are small and the sfermion masses are large. However, the Higgs mass tends to be too light in these regions, even when considering the calculational uncertainties on the Higgs mass. Thus, DUNE should be able to probe all the viable parameter space in this illustrative example. Fig. 16 The range of p → K + ν lifetimes found in the CMSSM (see Fig. 11) for the cases c = 0 and c = 0 (blue bands) compared with the sensitivities of the JUNO, Hyper-K and DUNE experiments. We also show results for p → π 0 e + in the CMSSM (see Fig. 17) with c = 0 (green band). The gray shaded areas are excluded by the Super-Kamiokande experiment [8,72]

A sub-GUT model
Our final example is a subGUT model with M in = 10 11 GeV. As in the CMSSM with M in = M GUT , in this case it is possible to set the dimension-five coupling c = 0, since the universality scale is below the GUT scale. The proton lifetime limit and its uncertainty for c = 0 are shown by the solid black contours, and the 1 − σ tot line is dotted. For c = 0, the limit and the 1σ tot uncertainties are in green and the proton decay constraint is significantly weaker.

Summary and discussion
We have analyzed in this paper the uncertainties associated with various phenomenological inputs in the calculation of the nucleon lifetime. We have used the minimal SU(5) GUT for this analysis, motivated by its relative simplicity, but in full knowledge of its shortcomings and the existence of more attractive alternatives. The considerations we have developed here could also be applied to any other specific GUT model, e.g., flipped SU(5) [131]. We have found that the largest uncertainties are those associated with lattice calculations of hadronic matrix elements, which have recently found significant changes in central values and reduced errors, and the strong coupling α s . We have also stressed the importance of using the appropriate value of sin 2 θ W in GUT calculations, while noting that its present uncertainty is of lesser importance. The most important quark mass uncertainty is that associated with m s , followed by m c , and we stress the importance of including one-loop mass renormalization effects. The most important CKM mixing uncertainty is that associated with the Wolfenstein parameter A. However, much larger uncertainties are associated with the GUT phases that are not observable in electroweak interactions, which can modify not only the dominant p → K + ν decay rates, but also modify significantly the branching ratios for other decay modes such as p → π + ν. However, our overall conclusion is that p → K + ν is the most promising decay mode for the next generation of underground detectors, particularly DUNE. We have also commented on the ambiguities in the proton decay predictions associated with the discrepancies between the masses of the charged leptons and charged-1/3 quarks, which warrant detailed study in specific models.
In this paper we have applied our analysis to variants of the minimal supersymmetric GUT with universality of the soft supersymmetry-breaking parameters imposed at the GUT scale (the CMSSM), above it (super-GUTs) and below it (sub-GUTs) (Fig. 15). The uncertainties reviewed in the previous paragraph, combined with our lack of knowledge of the possible masses of supersymmetric particles, make it impossible to be specific about the nucleon lifetime, even in such well-defined models as those as we have studied. However, our analysis shows that in all these models large regions of model parameter space with sparticle masses O(10) TeV can be explored with the upcoming generation of underground detectors. This is illustrated in Fig. 16, where we display the ranges of p → K + ν lifetimes found in the CMSSM (see Fig. 11) for the cases c = 0 and c = 0 (blue bands) compared with the sensitivities of the JUNO, Hyper-K and DUNE experiments. The gray shaded area is excluded by the Super-Kamiokande experiment. We also show results for p → π 0 e + in the CMSSM with c = 0 (green band), where we see that Hyper-K has some sensitivity, as discussed in the "Appendix" (see Fig. 17).
It will be interesting to apply the considerations presented here to other GUT models such as flipped SU(5), in which dimension-five baryon decay operators are absent, and the leading baryon-number-violating operators have dimension 6. Some remarks about dimension-6 proton decay are presented in the Appendix. We have also highlighted the ambiguities associated with models of the first-and secondgeneration quark and lepton masses. Interest in these and other issues in baryon decay will surely increase in the coming years as the start-up dates of JUNO, DUNE and Hyper-Kamiokande get closer. We trust that this paper will serve as a useful contribution to this coming trend. Lodovico

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: There is not data associated with these figures effectively. The are files with data in them but they just contain numbers and are not meaningful to anyone except the authors of this publication. Furthermore, the numbers are just the output of code and are effectively plotted all meaningful data.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix: Dimension-six proton decay
We describe in this appendix the calculation of the rate of proton decay induced by the exchange of the SU(5) gauge bosons, just for completeness. In this case, the relevant effective interactions are expressed by the following dimensionsix effective operators: i jkl + h.c., where O 6(1) We note that these coefficients have the identical phase factor, e iϕ i . As a result, this phase factor affects only the overall phase of the decay amplitude, and thus the decay rate is independent of this phase. At the one-loop level, 19 the RGEs of these coefficients can easily be solved [153,154]. The coefficients are then matched at the electroweak scale onto the effective operators where We again use the two-loop results given in Ref. [111] for the QCD RGEs. The partial decay width for p → e + π 0 is then given by with where A 1 and A 2 are the renormalization factors: The two-loop RGEs are given in Ref. [152].