Employing nucleon decay as a fingerprint of SUSY GUT models using \texttt{SusyTCProton}

While the observation of nucleon decay would be a smoking gun of Grand Unified Theories (GUTs) in general, the ratios between the decay rates of the various channels carry rich information about the specific GUT model realization. To investigate this fingerprint of GUT models in the context of supersymmetric (SUSY) GUTs, we present the software tool \texttt{SusyTCProton}, which is an extension of the module \texttt{SusyTC} to be used with the \texttt{REAP} package. It allows to calculate nucleon decay rates from the relevant dimension five GUT operators specified at the GUT scale, including the full loop-dressing at the SUSY scale. As an application, we investigate the fingerprints of two example GUT toy models with different flavor structures, performing an MCMC analysis to include the experimental uncertainties for the charged fermion masses and CKM mixing parameters. While both toy models provide equally good fits to the low energy data, we show how they could be distinguished via their predictions of ratios for nucleon decay rates. Together with \texttt{SusyTCProton} we also make the additional module \texttt{ProtonDecay} public. It can be used independently from \texttt{REAP} and allows to calculate nucleon decay rates from given $D=5$ and $D=6$ operator coefficients (accepting the required SUSY input for the $D=5$ case in SLHA format). The $D=6$ functionality can also be used to calculate nucleon decay in non-SUSY GUTs.


Introduction
Grand Unified Theories (GUTs) [1][2][3] continue to present an interesting framework for physics Beyond the Standard Model (BSM). The attractive features of this class of theories are the unification of all gauge interactions at some high energy scale, (partial) unification of matter into representations of a bigger symmetry, and an explanation for the quantization of hypercharge in the Standard Model (SM). The presence and type of low energy phenomena (e.g. light exotic matter) typically depends on the GUT model and as such these features are not universal predictions of gauge unification. An almost universal prediction of GUT models, however, is proton decay -and more generally, nucleon decay. The gauge bosons transforming as (3, 2, −5/6) under the Standard Model gauge group G 321 ≡ SU(3) × SU(2) × U(1) present in the adjoint representation of SU(5) violate baryon and lepton number (B and L) symmetries, and mediate the process of nucleon decay, generating B and L violating four-fermion operators, i.e. D = 6 proton decay. Other particles, such as scalar leptoquarks, can of course also mediate proton decay, but the main point is that as long as the unified group contains SU(5) as a subgroup, proton (nucleon) decay will be present.
One advantage of nucleon decay as a probe of new physics is that our experiments are able to measure decay rates corresponding to very large masses of the mediators, and current experimental bounds demand that the GUT scale M GUT must be higher than around 10 15 GeV. Since unification within field theory should happen below the Planck scale, this implies a window between roughly 10 15 GeV- 10 19 GeV for unification. The renormalization group (RG) running of gauge couplings in the SM can be made consistent with this window in the presence of appropriate intermediate states.
In supersymmetric (SUSY) theories, assuming the effective theory above the SUSY scale M SUSY to be the Minimal Supersymmetric Standard Model (MSSM), it is possible to construct B and L violating superpotential operators of four fields of the type qqql, where q stands for a quark and l for a lepton field. Such operators are generated, for example, when integrating out chiral supermultiplets of heavy color triplets (3, 1, −1/3) under G 321 . The resulting Lagrangian operators are dimension 5 operators, consisting of 2 fermions and 2 scalar superpartners of SM fermions. The D = 5 diagrams are then dressed at the SUSY scale to have SM external states, thus generating four-fermion nucleon decay operators. Contributions to nucleon decay generated in SUSY in such a way are referred to as D = 5 proton decay [4,5], and they typically dominate over D = 6 contributions, also usually requiring a higher GUT scale M GUT 10 16 GeV in SUSY GUTs compared to non-SUSY GUTs (due to the requirement of sufficiently large triplet masses).
A rough estimate of the total proton decay rate can be obtained using the formula where m p ≈ 0.94 GeV is the proton (nucleon) mass and the D = 5 amplitude A D=5 for the process can be estimated, see e.g. [6], as with α = g 2 /(4π), g the unified gauge coupling, y u and y d the Yukawa couplings of the up-and down-quarks, mw the gaugino mass, mf the sfermion mass, and m T the mass of the heavy color triplet mediating the decay. We assumed that the D = 6 contribution, e.g. from gauge bosons, is negligible, i.e. |A D=5 | |A D=6 |, since its amplitude is proportional to the square of the inverse heavy mass: with m X the mass of the gauge bosons mediating the B and L violating interaction. A thorough analysis of nucleon decay, however, should determine the decay rates into each channel separately, and both D = 5 and D = 6 contributions should be taken into account (we assume R-parity conservation, so no D = 4 proton decay). It is well known, see e.g. [7], that in SUSY theories the dominant decay channel is usually the kaon-neutrino channel p → K +ν . The proton lifetime bound for this channel currently comes from Super-Kamiokande [8]: 6.6 · 10 33 years at 90 % C.L. Future experiments are set to improve on this bound further: the expected 90 % C.L. bounds after 10 years of operation are 3.3 · 10 34 years for Dune [9], 3.0 · 10 34 years for Hyper-Kamiokande [10] and 2.0 · 10 34 years for JUNO [11].
To investigate the viability of a specific SUSY GUT model and determine its predictions, at least as far as nucleon decay is concerned, a more precise estimate of proton/nucleon decay rates is necessary than that of Eq. (1.2). While the procedure and expressions for the decay rates of the various nucleon decay channels are known, cf. [7,12], they are cumbersome to implement from scratch. A software package capable of the computation of nucleon decay rates in the various decay channels would thus be most welcome. Ideally, this tool would be able to handle all the complications in the computation that can arise, in particular both the sparticle spectrum and the flavor violating effects in the SUSY sector can play a prominent role. Furthermore, the effective D = 5 operators arising from integrating out the heavy triplets may also carry a flavor imprint from the Yukawa sector, since both the couplings to the heavy triplets and MSSM Higgs doublets arise from the same GUT operators. To be able to handle all of this automatically, the tool would thus ideally be built on top of an existing SUSY spectrum generator of sufficient sophistication.
In this paper we roll out exactly such a nucleon decay calculation tool for SUSY GUTs consisting of two in principle independent Mathematica modules SusyTCProton.m and ProtonDecay.m that we make publicly available: 4 1. The SusyTCProton package is a nucleon decay extension of the existing SusyTC package [13] built within the framework of REAP [14]. The standard SusyTC performs the RG evolution of the softly broken MSSM from the GUT scale to the Z-boson scale, automatically matching the MSSM and SM at an appropriate scale and generating all the SUSY sector information; the nucleon decay extension of this package performs, in addition, the 1-loop running of D = 5 operators from the GUT to the SUSY scale.
2. The ProtonDecay package takes as input the D = 5 decay operators and the sparticle masses and mixings, all at the SUSY scale, and dresses the relevant diagrams with SUSY sparticles to determine the effective four-fermion operators of nucleon decay, and after RG running to the nucleon scale computes the decay rates in 8 proton decay channels and 5 neutron decay channels (not counting different neutrino flavors, see Table 2 for the list). The input for this package is compatible both with a SusyTCProton output of internal values, as well as with an SLHA file [15] amended with the D = 5 operator information. This allows ProtonDecay to be used also with any spectrum generator that can output an SLHA file, such as SPheno [16,17], SOFTSUSY [18], SuSeFLAV [19] or SuSpect [20], provided that the user inputs D = 5 operator values from their own source. In addition, the ProtonDecay package can handle D = 6 contributions to nucleon decay, either as standalone terms (the non-SUSY GUT case) or in tandem with D = 5 contributions; the user needs to input them at a scale below which only the SM RG running is applied.
As a demonstration, we employ the above software packages in this paper to analyze the nucleon decay in two example toy models, cf. [21,22], both of which are of the flavor SUSY GUT variety. These models fit the SM fermion masses and mixings through particular choices of Yukawa textures and GUT-operators, which generate the Yukawa sector. The type of GUT operators used implies a choice of the Clebsch-Gordan (CG) coefficients between different sectors associated by the GUT symmetry. These same features impact also the D = 5 proton decay predictions, and it is precisely these flavor effects on nucleon decay that we focus on here. The magnitude of the decay rates in the various channels is subject to various potential ambiguities, such as the mass of the heavy triplet. These are overcome, at least in the models we consider here, by investigating ratios of decay rates, 4 A User's Guide for the installation and use of these packages is provided in Appendix C.
allowing to identify a nucleon decay "fingerprint" for each model. If nucleon decay events are seen experimentally in at least one decay channel, such a fingerprint technique, as we demonstrate in this paper, can be a powerful tool to either discriminate within a set of models with equally good Yukawa-sector fits and near-identical SUSY spectra, or to possibly rule out particular models of this type. The paper is organized as follows. We define the example flavor GUT models in Section 2, discuss how we implemented our numerical analysis in Section 3, present our results regarding nuclear decay in these models in Section 4 and give our conclusions in Section 5. The paper contains also 3 appendices. Appendix A presents a discussion on our doublet-triplet splitting setup and a few of its concrete realizations. Appendix B provides technical details on the nucleon decay computation and a note on conventions. Finally, Appendix C contains the documentation for our software packages SusyTCProton and ProtonDecay.
2 Choice of two example models for investigating nucleon decay

General setup in flavor SUSY GUTs
The models under consideration are flavor SU(5) SUSY GUTs of the type discussed in [21,22]. In this paper, we consider 2 toy models of this type and analyze their nucleon decay predictions.
The 3 generations of SM fermions are embedded in the standard way into SU(5) representations, where the index i runs from 1 to 3: The superpotential of the models can be separated into the Yukawa and breaking part where the breaking part consists of those terms, by definition, which do not contain the "fermionic" chiral representations 10 F i and5 F i . For each model we assume a particular texture of Yukawa matrices, to be specified later, and single operator dominance in each Yukawa entry, i.e. only a single GUT operator is present or dominates every Yukawa entry (see [22][23][24][25][26] for concrete realizations in the literature). This imposes strict relations between different fermion sectors and increases model predictivity. The Yukawa sector of each model thus consists of a set of (non-renormalizable) operators. We omit the high-energy model building in the neutrino sector. Also, we remain agnostic about the details of the breaking sector; the only requirement is that it breaks SU(5) symmetry to the Standard Model group G 321 := SU(3) × SU(2) × U(1) and includes the Higgs chiral supermultiplets assumed to be present due to their use in the Yukawa sector.
We are interested in models with as predictive D = 5 proton decay as possible. Since this mode of proton decay is mediated by color triplets (3, 1, −1/3) of G 321 , the details of doublet-triplet (DT) splitting come into the decay prediction. We thus consider as simple a scenario as possible for DT splitting, where only one triplet pair T 1 ⊕ T 1 couples to the fermionic sector, and for simplicity we assume that it is coming from 5 H ⊕5 H of SU (5).
This puts a limitation on which flavor GUT models can be considered in the context of this paper: the Higgs cannot be contained in the 45 or 45 of SU (5), and the non-renormalizable operators making use of such representations in the lists provided in [27,28] are not available, but only those containing 5 H or5 H .
At the GUT scale M GUT , the GUT models are matched to the softly broken MSSM as the effective theory, with the addition of the effective D = 5 proton decay operators in the superpotential. For the soft terms, we assume constrained MSSM (cMSSM, a.k.a. mSUGRA) boundary conditions at M GUT (see e.g. [29]), so they are parametrized by where m 0 is the universal soft scalar mass, M 1/2 the universal gaugino mass and A 0 the universal A-term factor. The choice of cMSSM is in principle not necessary, but it ensures that our analysis isolates the effects on nucleon decay from the flavor structure of the theory without added flavor ambiguities from the SUSY sector. At the level of the MSSM, the Yukawa sector (including the relevant operators for proton decay, but ignoring any neutrino terms) consists of where the color indices have been suppressed and the SU(2) contraction is written as Ψ·Φ = kl Ψ k Φ l , with kl the 2-index completely anti-symmetric tensor in the convention 12 = +1. Notice that we have written the Yukawa terms in the left-right (LR) convention. The triplet part of the Yukawa superpotential consists of whereâ,b,ĉ are SU(3) indices, while âbĉ and âbĉ are Levi-Civita tensors with 123 = 123 = 1. The SU(2) contractions are again labeled with a · dot. We refer to the 3 × 3 matricesỸ qq ,Ỹ eu ,Ỹ ql andỸ ud as quasi-Yukawa matrices in the remainder of this text.
This ansatz is specific to the models under consideration, since only the (anti-)triplets of 5 H ⊕5 H couple to the SM fermions; we thus omit any triplet indices I, J in the quasi-Yukawas. The general case with an arbitrary number of triplets is considered in Appendix B.1.1; the quasi-Yukawas written above would have I, J = 1, while all other quasi-Yukawas for I, J = 1 would be zero.
Integrating out the triplets and adapting the result from Appendix B.1.1 to our specific case, we obtain the following effective D = 5 proton decay operators: where the C 5 coefficients for our models equal denoting the "triplet effective mass". We provide a more detailed discussion on the form of the 2 × 2 matrix M T and the effective triplet mass in Appendix A. It is clear, however, that in this scheme the effective triplet mass does not necessarily correspond to any physical scale of particles, and as such can have a value above the Planck scale without losing consistency. In the SU(5) context, the Yukawa and quasi-Yukawa matrices come from the following GUT operators: where X is a placeholder for one or multiple Higgs fields in a Yukawa operator, possibly different for different matrix entries. The matrices arising from the same operator are related via SU(5) CG coefficients, which depend on the combination of Higgs fields X (entry by entry).

Specific models: their Yukawa and quasi-Yukawa sectors
We consider two models in this paper, referred to as model 1 and model 2. They in fact make use of the Yukawa textures of the models constructed in section 4 of [21]. We are interested only in the Yukawa sector here, so we remain agnostic about the details of doublet-triplet splitting (see Appendix A for some possibilities). The Yukawa textures of the two models are: model 1 : where * denotes a non-vanishing complex entry and Y 10 are always symmetric: Y 10 = Y T 10 . Note that in order to consider the two models on the same footing, we do not assume any spontaneous CP violation mechanism 5 of model 1, so the parameters are taken as complex.
The non-zero entries in these matrices are populated by the operators specified in Table 1  (chosen from Table 1 of [27]). The SU(5) representation in the index position indicates the type of contraction used between the pieces (analogous to integrating out the mediator in the index from a renormalizable interaction of that mediator with the fields in the parentheses). The vacuum expectation value (VEV) of the 24 H is at the GUT scale. This set of operators gives the following CG factor relations for the Yukawas and quasi-Yukawas, which we computed explicitly (the reader can also cross-check with [27]): • model 1: (2.16) • model 2: (2.18) The dot · indicates entry-wise multiplication of matrices.
3 Considerations and procedure for a numerical analysis

Implementation of the Yukawa sector
Following the considerations in Section 2, we implement models 1 and 2 at the GUT scale in the framework of the MSSM. The Yukawa sector is determined by specifying the two GUT Yukawa matrices Y5 and Y 10 , which is equivalent to specifying the up-and down-quark matrices Y u and Given the textures of the two models, Y u is a complex symmetric matrix, with the added constraint (Y u ) 13 = (Y u ) 31 = 0 on its entries in the case of model 1. Any symmetric complex 3 × 3 matrix Y can be decomposed, using Takagi decomposition, as where the singular values y i ≥ 0, and U is a 3 × 3 unitary matrix. For the parametrization of U, we shall make use of the definitions where U 12 , U 23 and U 13 represent unitary matrices of 2-D complex rotations, each in principle depending upon an angle and a phase. The most general 3 × 3 unitary matrix can be explicitly written by using the following family of real parameterizations: where the family parameter δ 0 is fixed, and the 9 free parameters describing the unitary rotation consist of angles θ ij and phases δ, φ i and η i : The first and last factor of the product in Eq. (3.5) are the two diagonal matrices containing the "outer phases" η i and the "inner phases" φ i (with inner/outer referring to their location in the product of Eq. 3.1). There are only 2 φ-phases, since a third such phase can be absorbed as a common phase factor into the η phases. Furthermore, the 1/2 factors for the φ-phases are there for later convenience; we shall view them as phases of the singular values in the Takagi decomposition rather than as part of the unitary rotation.
Due to greater convenience, we use different parameterizations in the two models. In model 1, we set δ 0 = π/2 and impose the additional (complex) condition (Y u ) 13 = 0 by solving for the real variables θ 13 and δ. For model 2 we set δ 0 = 0. The advantages of this choice of δ 0 will be discussed later for each model separately.
The use of Takagi decomposition for the parametrization of Y u has the advantage that it provides direct control over the singular values (Yukawa couplings) and mixing angles in the upquark Yukawa matrix. Also, note that we can redefine each fermionic representation 10 F i and 5 F i with an arbitrary complex phase, implying that 6 phases in the Yukawa sector are unphysical. In our convenient parameterization of Y u , a phase redefinition of 10 F i allows us to fix the outer phases without loss of generality: We shall use the phase redefinitions of5 F i to remove the phase of one complex entry in each column of Y d . In model 2, this implies that Y d is diagonal with (positive) real values, while in model 1 a complex phase remains in one 2nd column entry; we choose that entry to be (Y d ) 12 .
Putting all the above considerations together, the final form of the parameterizations of the Yukawa matrices Y u and Y d , which uniquely determine the entire Yukawa sector, are the following: • Model 1: According to the textures in Eq. (2.13), the parametrization of the down-and up-quark Yukawa matrix in the context of the MSSM at the GUT scale is given by where the values of θ uL 13 and δ depend on the other parameters of the up-sector; they are determined from demanding (Y u ) 13 = 0: 6 (3.12) Once the numeric values of the free parameters are specified, the C coefficients in Eq. (3.12) become complex numbers, allowing for the complex equation (3.11) to be solved numerically.
Given the hierarchical structure of masses and mixings in the quark sector, this equation can always be solved for δ and θ uL 13 in a numerically stable fashion. The values for the angle are very small: θ 13 10 −5 , which suggests a possible alternative which we mention here in passing: we could omit solving Eq. (3.11) and simply take θ uL 13 = 0 (δ dependence then disappears). This would be a relatively good approximation yielding |(Y u ) 13 | 10 −5 ≈ 0. 6 We omit an overall cos θ uL 13 factor in the (Yu)13 = 0 condition. Such a factor allows for an alternative solution θ uL 13 = π/2, which also sets (Yu)23 = 0, so we reject it.
The CKM matrix comes from diagonalizing both Y u and Y d . It turns out that for this particular texture, choosing the δ 0 = π/2 parametrization for the up-sector then already provides a good CKM phase for ϕ ≈ 0 in Y d .
Also, note that the complex phase φ here is in the (Y d ) 12 entry, in contrast to the class of models in [24] where a phase is in the (Y d ) 21 entry. Despite the same texture, the difference comes from removing the phases in Y d using5 F i in our case (removes one phase per column), while [24] prefers to remove phases via 10 F i redefinitions (removes one phase per row), so that they can make use of the5 F i freedom in their implementation of the neutrino sector, which we do not consider in the present model.
The charged lepton Yukawa matrix Y e and the quasi Yukawa matricesỸ qq ,Ỹ eu ,Ỹ ql and Y ud for model 1 are determined from Eq. (2.15) and (2.16).
• Model 2: Corresponding to the textures in Eq. (2.14), the parametrization of the down-and up-quark Yukawa matrix in the context of the MSSM at the GUT scale is given by where while Y e ,Ỹ qq ,Ỹ eu ,Ỹ ql andỸ ud are then determined from Eq. (2.17) and (2.18).
In model 2, the CKM matrix is controlled entirely by the up-sector: choosing the δ 0 = 0 parametrization in the up-sector thus provides the standard parametrization of the CKM matrix and the 3 angles and phases in U 2 are directly those of the CKM values at the GUT scale.

Model parameters
Given the implementation/parameterization of the Yukawa sectors of the two models in the previous subsection, we can now turn our attention to comprehensively analyzing the parameters and observables. The two models contain the following input parameters at the GUT scale: • Model 1 contains the input parameters where the Yukawa sector parameters y u,d , θ uL , ϕ and φ are used for the constructions in Eqs. (3.8)-(3.12).
• Model 2 contains the following input parameters at the GUT scale: where y u,d , θ uL , δ and φ are used in constructions of Eqs. (3.13)-(3.15).
Both models consist of 16 input parameters, 12 of which construct the Yukawa sector, and 4 parameters are intrinsic to the cMSSM: m 0 , M 1/2 , A 0 specify the mSUGRA boundary conditions at the GUT scale, and tan β is the ratio of the MSSM Higgs VEVs. Note that the two phases φ 1 and φ 2 in both models do not change the SM prediction of fermion masses and mixing angles, since they get removed in the construction of the CKM matrix, while they do have an impact on the quasi-Yukawa sector and hence proton decay predictions. These phases are typically referred to as "GUT phases" in the literature [33,34].
All initial parameters are real. Taking into account that singular values in a Takagi decomposition are non-negative, and all quark masses are non-vanishing, the parameters in the Yukawa sector are allowed to be in the following ranges (for both models): The cMSSM parameters can in principle also be varied in a fit. In this paper, however, our goal is an analysis which disentangles the flavor effects from the SUSY spectrum effects on proton decay. For this reason, we fix the cMSSM parameters to common example values in both models, where the SUSY threshold corrections allow for a good fit to (charged) fermion masses and the Higgs mass. We thus consider only a particular cMSSM benchmark point in our analyses.
To simplify the fitting procedure in the two models, we fix the gauge couplings at the GUT scale M GUT = 2 · 10 16 GeV: g 3 = 0.698, g 2 = 0.697 and g 1 = 0.704, which is consistent with the experimental values at low energies (assuming MSSM running above the SUSY scale of around 3 TeV). Furthermore, since the effective triplet mass scale M eff T only affects the nucleon decay widths by an overall scaling, i.e. Γ ∝ 1/(M eff T ) 2 , we can choose a fixed value M eff T = 10 19 GeV in the following analyses, and can easily relate the results to any other value. Finally, in order for the RGE of the ratio y µ /y s of the 2nd family to start at 6 at the GUT scale (based on our two models) and run to the experimental value at the scale M Z , we choose the sign of the µ coupling (connecting H u and H d ) in the MSSM to be sign µ = +1, such that the dominant gluino loop diagram in SUSY thresholds drives y s in the correct direction [27,35].

Observables
Both models are fit to the same 14 SM observables at low energies: , δ CKM , y e , y µ , y τ , m h , (3.22) corresponding to the masses of fermions in the up, down and charged lepton sector, the CKM parameters of the quark sector, and the SM Higgs mass m h . Although there are 16 input parameters in both models, 2 GUT phases φ 1 and φ 2 have no influence on the listed observables, but only on nucleon decay discussed later, so the fit has effectively only 14 relevant parameters. Since that number equals the number of listed observables, a reasonable expectation is that a good fit for both models can be obtained, with no flat directions in the χ 2 .
The values of the SM Yukawa observables at M Z and the error bars are taken from [36], where for the Yukawa couplings of the charged leptons we take a standard deviation of 1 % of the central value, which roughly corresponds to the precision of the running. The experimental value of m h is given in PDG (2018) [37], where we take a standard deviation of 3 GeV, namely 125.6 ± 3 GeV, due to the theoretical uncertainties present in the calculation (for which we use FeynHiggs, see Section 3.4).
Other observables of interest predicted by the model are the nucleon decay rates. We use our extended SusyTC code SusyTCProton to compute 13 different decay channels, 8 of which are proton decay channels and 5 are neutron decay channels (we sum over the neutrino flavor states, since the decay experiments do not distinguish them); they are listed in Table 2 along with their experimental bounds (we took the most strict bound we found across multiple experiments). We shall not consider these observables in a χ 2 fit, since the decay rates themselves are only predicted up to an overall factor, unless the effective triplet mass parameter M eff T is also known. Taking a high enough M eff T can always avoid all proton decay bounds for D = 5 contributions (we remind the reader that this parameter in our setup is not directly connected to any physical scale of triplets in the theory, so it can be taken above the Planck scale). Despite the lack of predictivity in nucleon rates themselves, our setup does allow for unambiguous predictions of branching ratios and, more generally, ratios of decay rates in different channels, e.g. a ratio between a proton and neutron decay channel. The 90 % confidence level exclusion bounds on proton and neutron decay from the Super-Kamiokande experiment for various (B − L)-preserving decay channels with a meson and a lepton in the final state. For each decay channel the corresponding lifetime bound τ /B and the translated partial decay width Γ partial are listed, where B is the branching ratio for the decay channel. The data is taken from Figure 5-3 in [38] and further updated with other sources [8,[39][40][41][42][43].

Calculational procedure
For given values of the parameters in model 1 and 2, listed in Eq. (3.16) and (3.17) respectively, the Yukawa matrices and the soft terms are implemented at the GUT scale M GUT = 2 · 10 16 GeV according to Section 3.1. Moreover, the dimension 5 operators which are relevant for nucleon decay are calculated at M GUT as well, using the formulas from Eq. (2.8) and (2.9). The RG running in the MSSM is performed by means of the Mathematica package SusyTCProton (cf. Appendix C.1), which is an extension of the Mathematica package SusyTC [13] used within REAP [14]. The running of the Yukawa matrices, gauge couplings and soft terms is computed at 2-loop, whereas the running of the dimension 5 operators is computed at 1-loop.
The SUSY scale M SUSY is determined dynamically in SusyTCProton by the geometric mean of the two up-type squark masses which have the largest mixing with the stop flavour eigenstatest 1 andt 2 , namely M SUSY = √ mt 1 mt 2 . The values of the Yukawa matrices, gauge couplings and soft terms at M SUSY are used to calculate the partial decay widths of the proton and neutron listed in Table 2 by means of the Mathematica package ProtonDecay (cf. Appendix C.2). In addition, the electroweak (EW) Higgs mass is calculated by using FeynHiggs 2.16.0 [44][45][46][47][48][49][50][51].
In SusyTCProton, the MSSM is matched to the SM by taking the SUSY threshold corrections into account. Furthermore, the running from the SUSY scale to the Z-boson mass scale M Z = 91.2 GeV is calculated at 2-loop. At M Z , the observables in the quark and charged lepton sector, namely the Yukawa couplings and the CKM parameters, are computed.
For the statistical model analyses we use the following procedures: • As a measure of goodness of fit we define the χ 2 function in the usual way: where the vector x consists of the input parameters of the models from either Eq. (3.16) or (3.17), and the index i goes over all observables. The y i denote the central values from the (experimental) data, σ i are their corresponding standard deviation errors, while f i ( x) are the predictions for the i-th observable given the parameter point x. Some observables may be equipped with asymmetric errors σ i+ and Each χ 2 i term in the χ 2 function is referred to as a pull associated to a particular observable; it will also be convenient to refer to a particular pull by replacing the index i with the label of the i-th observable.
• The best-fit point for a given model setup is calculated by using a differential evolution algorithm for global minimization using custom Mathematica code.
• To calculate the posterior density of the parameters and observables for a certain model setup, we perform an MCMC analysis by using an adaptive Metropolis-Hastings algorithm [52]. All prior probability distributions for the parameters are chosen to be flat. Furthermore, we calculate 12 independent chains with 5 · 10 4 valid points (beyond the "burn-in phase"), which gives 6 · 10 5 points in total for each dataset. All chains are started at the best-fit point of the corresponding model setup.

Results for nucleon decay from SusyTCProton
In this section, we present the results of our numerical analyses for the two models we considered.
Our main interest are the predicted nucleon decay rates for the points with a good fit to the low energy data.
Expressions for nucleon decay rates mediated by D = 5 operators contain many contributions, all of which depend on both the masses of the SUSY sparticles that are used to dress the D = 5 operators, as well as flavor effects both from the SM Yukawa sector and the sfermion sector. More details on the expressions used are given in Appendix B and references specified there.
The SusyTCProton.m and ProtonDecay.m packages developed for such a nucleon decay calculation, see Appendices C.1 and C.2, can handle all these effects. Our main interest in this paper are SUSY GUTs, and as such we are interested in the effect the Yukawa textures are having on proton decay. For this reason, we would like to isolate this effect from other contributions: taking cMSSM boundary conditions already insures no extra flavor violation in the SUSY sector, while the effects of variable sparticle masses can be removed by fixing the soft parameters.
We therefore choose in our analysis to fix the cMSSM parameters and consider only a benchmark point. We take the following values:  We use these values in both models, and they allow us to obtain good fits for both 7 , thus being able to compare them on an equal footing and isolate the flavor effects. We have checked that EW vacuum stablility for this parameter point by using Vevacious 1.2.03 [53], SPheno 4.0.3 [16,17], and the predefined model from SARAH 4.14.1 [54,55] with possible charge breaking via stau VEVs. The output shows that the produced EW vacuum is stable. Fixing the cMSSM parmeters leaves 12 free parameters to vary for a fit, with 10 having an effect on the fit, while varying the GUT phases φ 1 and φ 2 will only effect the nucleon decay rates.

Best-fit points and nucleon decay
Minimizing the χ 2 with respect to the flavor parameters, we obtain a best-fit point for each of the models. The parameter values for the best-fit points are shown in Table 3, while the χ 2 results and various dominant pulls are shown in Table 4. Most of the tension in the fit is coming from the pull in the s-quark mass.   We can now use these best-fit points to obtain predictions of proton and neutron decay rates in various decay channels. We postpone a more comprehensive comparison of all the channels to Section 4.2, and will only focus here on some selected and most interesting aspects of the results, namely the dependence of the rates and their ratios on the GUT phases φ 1 and φ 2 , since these phases are not determined by the fit and we thus have an embedded ambiguity of their values in the models. We plot the results in this section as contour plots in the φ 1 -φ 2 plane, and take a fixed value for the effective triplet mass M eff T = 10 19 GeV where needed to obtain specific numerical results. Each of these plots is equipped with it's own scale legend to the right. Note that despite the same color range for contour regions, the values represented by these colors change between different plots. Also, some of the plots provide information about the location and value of the minimum and maximum, which are added to the plots in red color. This information of the minimum and maximum is important, since it indicates that the predictions robustly lie in a compact region, and the decay rates cannot be made to vanish by a suitable choice of φ 1 and φ 2 .
To draw the plots, we computed the values on an equidistant 2D-grid with increments π/30 for both the φ 1 and φ 2 directions, while a separate minimization/maximization routine was run to obtain the extrema in the φ 1 -φ 2 plane.
We make the following observations of the results:  • The dominant decay channels of both the proton and neutron are those with neutrinos in the final state; furthermore, the decay channel with a kaon always dominates over that of the pion. We plot the GUT phase dependence of the dominant proton and neutron decay rate for each model in Figure 1.
We see that the rates for each process in a given model can differ by a factor of up to ∼ 3, depending on the unknown GUT phases, but can never vanish completely. Observing the decay rate values in each panel, we see that the dominant neutron decay channel dominates over the proton decay channel in the case of model 1, while the opposite is true in model 2. Also, for model 1 the φ 1 -φ 2 regions minimizing the proton decay rates roughly correspond to those minimizing the neutron decay rate; similarly, the proton and neutron decay rates are simultaneously maximal in a similar region. For model 2, on the other hand, minimizing the dominant proton decay roughly maximizes the neutron decay rate, and vice versa.
• We selected to show another decay channel, p → K 0 e + , in Figure 2.  lepton mass matrix is not diagonal, the relative changes of the rates are smaller, but the φ 2 dependence is not negligible.
• A result independent of the effective triplet mass can be obtained by considering ratios of decay rates. We choose to present two such ratios in the plots: the decay ratio Γ p→K +ν /Γ p→π +ν represents the ratio between the two most dominant proton decay channels, while the ratio Γ p→K +ν /Γ n→K 0ν is between the dominant proton and neutron channels. The GUT phase dependence of the ratios is shown in Figure 3.
We observe that the plots for ratios have a similar feature to those for the dominant rates: the ratios are simultaneously minimized or maximized for model 1, but minimizing one ratio roughly maximizes the other in model 2. This implies that the process p → π +ν , which is not shown in Figure 1, has a similar minimizing/maximizing region to the dominant proton decay process; the switch of the regions happens for model 2 only when migrating to neutron decays.
Another important observation is that the predicted range for a given ratio does not overlap when comparing the two models. This implies that experimentally determining such a ratio of nucleon decay channels would enable us to discriminate between the two flavor GUT models, despite the unknown GUT phases and despite the two models having the same soft SUSYbreaking parameters. This difference in the predicted ratios is thus purely a flavor effect, as could be expected; we provide additional discussion of using ratios for discrimination in Section 4.2.

MCMC results for nucleon decay
The nucleon decay results of the previous section involved varying only the two GUT phases φ 1 and φ 2 , while we keep the other flavor parameters confined to the best-fit point in each model. We now analyze the two models by varying all flavor parameters of Eq. (3.16) and (3.17), such that the fit is still close to the best-fit point; the cMSSM prameters are still fixed, however, to their benchmark values of Eq.   The plots in this section show HPD intervals of nucleon decay rates or ratios predicted by the two models. The HPD intervals are colored red (green) for model 1 (2). We show intervals corresponding to the 1-σ and 2-σ HPD regions as darker and lighter colored, respectively, and overlayed one on top of the other.
We make the following observations of the results: • Our SusyTCProton and ProtonDecay packages are able to compute 13 different decay channels (specified in Table 2). The MCMC results for the HPD intervals in each channel are shown in Figure 4. Note that we assumed the effective triplet mass to be M eff T = 10 19 GeV; the decay rate results for a different M eff T can be easily obtained by noting that Γ ∝ (M eff T ) −2 , implying that we linearly translate the intervals down by 2 on the log scale for each order of magnitude that M eff T increases over the reference value. The experimental bounds for the various channels are shown as blue bars; choosing a different M eff T can control whether we violate any of the experimental bounds or not with the predictions. For the reference value of M eff T , the experimental bound on p → K +ν would exclude model 2 but not model 1. We see that the dominant channels are those with a neutrino in the final state, while those with a charged lepton have decay rates at least 4 orders of magnitude below the dominant ones. The dominant proton decay channel is p → K +ν , while the dominant neutron decay channel is n → K 0ν . Which of the two is bigger differs between the two models. This result is consistent with the dominant decay channels when only the GUT phases were varied in Section 4.1.
Also, we notice that when a charged lepton is a decay product, the predictions for the rates are much sharper when the charged lepton is e + (the 2-σ ranges vary ∼ 20 % around the average), while the prediction for µ + in the final state typically spans 1 to 2 orders of magnitude. Given the decay rates of these processes, however, they are at least 4 orders of magnitude below the dominant channels and thus experimentally not within reach.
Comparing the HPD intervals between the two models process by process, we notice that for a given M eff T most channels are predicted with a similar decay rate, except for p → K +ν ; this process is thus most suitable for discrimination between the two models, assuming another channel has been measured as well. Luckily, the decay p → K +ν also has either the biggest or 2nd biggest decay rate, so it would likely be among the first experimentally measured channels if one of the considered models is realized in nature.
• The quantities independent of M eff T are the ratios of various decay rates. Note that a ratio for each parameter point in the MCMC dataset is computed separately, and only then is statistics performed on these values; this is an important point, since there could be additional correlations between various decay channels in a model not apparent when viewing only the HPD intervals of decay rates, i.e. when considering the numerator and denominator in the ratio separately.
In Figure 5 we show selected ratios which best discriminate between the two models. They effectively have no overlap of predicted 2-σ HPD intervals when comparing the two models, thus measuring even one such ratio has sufficient discriminating power. All ratios between the 13 channels not shown in the Figure all have some overlap in the 2-σ intervals. More broadly, given any specific model, predictions for a set of ratios effectively provides a fingerprint of that model, in principle allowing for it to be excluded if multiple channels of decays were measured experimentally.
The most interesting ratios experimentally are those involving the dominant decay channels, which would be the first ones measured. Also, a very big or small prediction of the ratio (compared to a 1 : 1 ratio) implies a very big discrepancy in the rates of the numerator and denominator, again making it unlikely that both processes could be measured in the near future. Given these considerations, the most experimentally relevant ratios are the first three, which compare the high rate processes with a neutrino in the final state.
All but one of the ratios in Figure 5 include the p → K +ν decay as one of the two processes, confirming our comment for Figure 4 that this process is a good discriminator between the two models. In terms of the sharpness of the prediction of ratios, the value for the fourth ratio in the figure has the smallest allowed range: the 2-σ HPD intervals for Γ p→K 0 e + /Γ p→η 0 e + are [5.7, 6.9] and [7.5, 9.3] for models 1 and 2, respectively. The decay channels in that ratio have very small rates and are experimentally out of reach in the considered scenarios; assuming an experimental measurement of such a ratio (or even just one channel and having bounds on the other) provides, however, an excellent way to exclude the model under consideration.
A final note on ratios: some ratios are predicted exactly the same with a value (1/ √ 2) ±1 for every point due to the form of the low-energy effective Lagrangian in chiral perturbation theory. In particular, these are the 3 ratios where only a replacement of π 0 with π ± takes place in the final state: These ratios are fixed and have no dependence on any high-energy or SM parameters. HPD intervals for nucleon decay rates Γ in different decay channels, M T eff = 10 19 GeV  Predicted HPD intervals for selected Γ 1 /Γ 2 ratios of nucleon decay channel rates Figure 5: The 1-σ (dark) and 2-σ (light) HPD regions predicted via MCMC for selected ratios of decay rates of two channels; we assume the same cMSSM benchmark point in both models.

Conclusions
While the observation of nucleon decay would be a smoking gun of Grand Unified Theories (GUTs) in general, the ratios between the decay rates of the various channels carry rich information about the specific GUT model realization. We investigated in this paper these nucleon decay fingerprints of two flavor SUSY GUT models based on the group SU (5). The two models differ in the texture of their Yukawa matrices, and the types of GUT operators generating the SM fermion masses, cf. Section 2. Both admit essentially equally good fits to the low energy Yukawa data of SM masses and mixing angles (with a χ 2 ≈ 5, cf. Section 4.1), even given the same sparticle spectrum. Yet their nucleon decay computation shows that the partial decay rates of various channels are good observables to discriminate between them, despite the decay ambiguity introduced by 2 unknown GUT phases. We investigated the dependence of the decay rate on the GUT phases at the best-fit point for the flavor parameters of the model in Section 4.1. The nucleon decay fingerprinting has been performed by computing the decay rates in 8 different proton decay channels and 5 different nucleon decay channel, cf. Table 2. We obtained their predicted ranges by varying the flavor parameters in a way still compatible with a low χ 2 of the low-energy observables using the MCMC algorithm. The main results of the analysis are presented in Section 4.2, see especially Figures 4 and 5. The general features are those expected in a SUSY GUT: the p → K +ν and n → K 0ν channels are the dominant modes of proton and neutron decay, respectively, and the rates with kaons in the final states are bigger than those with pions when the other decay product is the same. Other features found in the fingerprints are that the neutrino decay channels dominate over those with charged leptons in the final state and that the channels with e + in the final state are predicted much more sharply than those with µ + .
The decay rates invariably depend on the scale of the heavy color triplets that mediate D = 5 nucleon decay. We considered a doublet-triplet setup where the overall scale of the rates depends only on the triplet effective mass parameter M eff T . Robust and M eff T -independent predictions for both models are obtained by considering ratios of decay rates between different channels. The key connection of nucleon decay to flavor is that the couplings of the SM fermions to the D = 5 nucleon decay mediating triplets are related to the Yukawa coupling of the theory, since they arise from the same GUT operators, which is a generic feature of any GUT realization. The results for the considered models show that the impact of their different flavor structure is already sufficient to discriminate between them. In particular, the dominant proton decay mode p → K +ν turns out to differ relative to other dominant modes with neutrinos in the final states such as p → π +ν and n → π 0ν .
The computations of nucleon decay rates for the 13 channels in this paper were performed by use of the Mathematica software modules SusyTCProton and ProtonDecay, which we developed for the purposes of this paper and are rolling them out for public use. The SusyTCProton is an extension of the SusyTC package within the REAP framework, and allows to compute the RG running of D = 5 operators from the GUT scale to the Z scale alongside the running of all other parameters in a softly broken MSSM, performing the MSSM to SM matching and dressing the D = 5 operators at the SUSY scale. The ProtonDecay module admits as input either the SusyTCProton output or an SLHA file amended with D = 5 operator values at the SUSY scale, and computes the resulting nucleon decay rates in the various channels. The module is able to include D = 6 contributions to nucleon decay if the operator values are provided at some high scale where the SM is taken as a valid effective theory, so the package can be used for computing nucleon decay in non-SUSY GUTs as well. We list the documentation for the use of these software packages (a "User's guide") in Appendix C.
When used in tandem, the above software packages allow for the computation of the decay rates in any softly broken SUSY GUT if one is able to input the parameter values at a scale where the MSSM becomes valid (in our case this was immediately below the GUT scale). It computes the rates in all generality, including all the flavor effects from both the Yukawa and SUSY sectors, and considers all diagrams in these decays following [12]. This allows for the fingerprinting of nucleon decay within such a model, opening up a new avenue for the model to be tested in experiments. The fingerprint is most robust when considered at the level of ratios of decay rates for various channels. Experimentally, the measurement of decay events in just one channel either for the proton or neutron, alongside experimental bounds in the other channels, already allows for strategies for crosschecking with the fingerprint, since the model prediction might suggest that decay events in some other channels should have been observed already as well. A situation with decay events measured in multiple channels of course allows for even better tests of the model. These fingerprinting techniques can be readily used either to discriminate between otherwise indistinguishable models from the point of view of current experiments, or to possibly rule them out.
This paper represents a thorough analysis of nucleon decay in two flavor GUT toy models, demonstrating the efficacy of using their nucleon decay fingerprints to distinguish them. Alongside this analysis, we also roll out new software tools which can be employed by the community for further investigation of nucleon decay in models of various types, thus adding to the nucleon decay fingerprint collection.
In this appendix section we present two concrete realizations of such a scenario; despite multiple possibilities for the underlying dynamics, only the effective mass M eff T is the relevant parameter. Note that while this parameter has mass dimension 1, it does not (necessarily) represent a physical mass of a particle or a characteristic energy scale at which a process takes place, and as such can take trans-Planckian values. Indeed, for our numerical analysis in Section 4 we took for the effective triplet mass the benchmark value M eff T = 10 19 GeV, so the proposed mechanisms should be able to accommodate such a high parameter value.
We present two possibilities how to achieve a setup admitting high M eff T : a modification of the ordinary missing partner mechanism is discussed in Appendix A.1, while the double missing partner mechanism is described in Appendix A.2.
A.1 The modified missing partner mechanism A.1.1 The missing partner mechanism -the usual setup In the usual setup, the missing partner mechanism in SU(5) [57,58] requires that the color triplets T (or anti-triplets T ) come from the following SU(5) representations: The Higgs doublets (1, 2, ±1/2) on the other hand are present only in the 5 H ⊕5 H part; as mentioned earlier in Appendix A, the representations 50 (50) of SU(5) remarkably contains only a triplets (anti-triplet), but not an (anti-)doublet.
The doublet and anti-doublet directly correspond to the MSSM doublets H u and H d , respectively: The (anti-)triplets in Eq. (A.1) on the other hand are labeled by We assume then an SU(5)-level superpotential of the form where the G 321 singlet component in 75 H obtains a GUT scale VEV v = 75 . 8 Crucially, the explicit mass term 5 H · 5 H is missing (this can be achieved for example by imposing discrete symmetries). The superpotential of Eq. (A.5) yields the mass terms where the doublet and triplet mass matrices M D and M T are given by In general we have λ,λ, v, m ∈ C, but the phases of λv,λv and m can be absorbed into the overall phase definition of the irreps 50 H , 75 H and 50 H in Eq. (A.5), respectively. This implies that all entries in M T can be taken to be real numbers.
This situation in Eq. (A.7) clearly achieves DT-splitting, since the doublet remains effectively massless, while all triplets are roughly at the GUT scale due to the presence of the VEV v. This turns out to present a problem, since introducing a lot of matter below M GUT may increase the gauge coupling value to such a degree at M GUT that the Landau pole of the coupling is unacceptably close to the unification scale. Let's investigate this situation: 1. The solution to the 1-loop RG running of the gauge coupling g between two scales µ r and µ r is described by where α = g 2 /(4π) and β depends on the model. We remind the reader that in SUSY theories β = 3D 2 (V ) − D 2 (C), where D 2 is the Dynkin index, with V and C referring to all vector and chiral supermultiplets in the theory with masses below a given renormalization scale, respectively. Adding a lot of matter (chiral supermultiplets) makes β smaller (more negative), causing g to blow up sooner when running to higher scales.
2. Suppose we change the theory by moving some representations from the scale µ 0 to µ 0 /c (c > 1), thus changing the β-function between µ 0 /c and µ 0 from β to β . The new coupling α −1 and the old coupling α −1 should be the same at µ 0 /c, and using Eq. (A.9) we can relate them at the scale µ 0 via which depends on the difference of the two beta functions ∆β = β − β. For simplicity, we consider the unified coupling and a setup where an entire 50 ⊕ 50 is moved from µ 0 = M GUT to µ 0 /c: since it is a complete SU(5) representation, it does not change the unification scale µ 0 , but only the value of the unified coupling from α −1 (µ 0 ) to α −1 (µ 0 ).
3. The scale of the Landau pole µ L is defined implicitly by α −1 (µ L ) = 0, i.e. when g has a pole. We define C as the factor above the GUT scale where the Landau pole occurs, i.e. µ L = Cµ 0 . The Landau pole condition in the new theory α −1 (Cµ 0 ) = 0 can be related to α −1 (µ) by using the solution in Eq. (A.9) (assuming β f for the beta function above µ 0 and no new thresholds); when combined with Eq. (A.10) this relates c, C and the original gauge coupling α −1 at the GUT scale: . We now turn to our specific case of SU(5) and assume the full theory has at least the following chiral supermultiplets in the fermionic and Higgs sectors (no messengers, no flavons, we assume GUT breaking only with the 24 and 75): This yields β f = −57 (see [59][60][61] for tables of Dynkin indices), and assuming g(µ 0 ) ≈ 0. Even for a nearby Landau pole a factor only C = 10 above M GUT , we can move the 50s only a factor c ≈ 2.35 below the GUT scale.
The above calculation implies that to avoid a Landau pole too close to M GUT , the 50 cannot really be much lighter than the GUT scale. A small m 50 alone thus cannot be used to enhance M eff T by a large factor in the usual missing partner mechanism.
To avoid problems with perturbativity, we need to modify the missing partner mechanism so that the 2-2 element of M T has a richer structure. We can for example introduce representations 24 H and 75 H (possibly those that we already have of this dimension, depends on the full model), and two new terms to the superpotential W DT in Eq. (A.5): 24  The expression form 50 can now have all terms at the GUT scale, but the couplings λ 24 and λ 75 can be tuned such thatm 50 itself for the triplets is small, while keeping all SM representations in the 50 H ⊕ 50 H other than the triplet heavy. This can be seen from the (relative) CG factors in front of the 3 operators generating their masses that we computed explicitly in Table 5. We can see from this table that the CG factor pairs for the operators (λ 24 , λ 75 ) in each row are unique; since the mass of any SM irreducible representation in the 50 H ⊕ 50 H has the form of Eq. (A.16), but with CG factors added to the λ 24 and λ 75 terms, the uniqueness of the CG pairs implies that any one SM irreducible representations can be tuned to be light while keeping the others heavy; we choose to tune it to a low value for the triplets.
In this modified missing partner mechanism, we achieved what we wanted: it is possible to perform DT-splitting, as well as tune M eff T to arbitrarily high values. All physical masses of particles in the Higgs sector are around the GUT scale (except the MSSM doublet pair), while precision of the tuning in the parameterm 50 controls how high the effective triplet mass M eff T is: For λλv 2 ≈ (10 16 GeV) 2 , and effective triplet mass |M eff T | = 10 19 GeV used as a benchmark value in this paper, we requirem 50 ≈ 10 −3 M GUT , i.e. a relatively small tuning of one part in 10 3 , which retains the advantage of the missing partner mechanism over a brute force fine-tuning of a doublet from M GUT to the EW scale. Table 5: Clebsch factors (up to an overall VEV normalization factor) for all different SM representations in the 50 (or 50) of SU(5) in the three different operators generating their masses.

A.2 Double missing partner mechanism
A straightforward alternative to the modified missing partner mechanism discussed in Appendix A.1 is the double missing partner mechanism, cf. [21]. This also solves the DT splitting problem and provides only one effective triplet mass that can be trans-Planckian. In this scenario, the representations of the doublets/triplets are doubled: the DT sector consists of where the G 321 doublets (1, 2, ±1/2) are contained in .20) and the triplets T I ∼ (3, 1, −1/3) and anti-triplets T I ∼ (3, 1, +1/3) are in We then arrange the presence of such operators that the mass terms have the following specific form for the mass matrices: where α i are dimensionless couplings, while V could for example be a VEV V = 75 at the GUT scale (see [21] for more details on some realizations). In this set up, we require that only 5 H ⊕ 5 H (and not those with primes) can couple to the SM fermions in 10 F i ⊕ 5 F i , which yields the effective triplet setup with the effective mass parameter equal to To postpone the appearance of the Landau pole from the large (negative) contributions to the beta function, we can provide the representations 50 ⊕ 50 ⊕ 50 ⊕ 50 with large mass terms close to the Planck scale, i.e. we take the hierarchy We obtain from this setup the physical doublet masses to be 0 (the MSSM doublet) and µ , while the physical triplet masses are at scales η ±1 V (two triplets at each scale), where η ≡ V /M 50 . This constitutes a successful DT splitting yielding allowing for the effective mass to be large provided that µ is small. For V = 2 · 10 16 GeV, η 2 ≈ 10 −4 (so that M 50 and M 50 are close to the Planck scale) and |M eff T | 10 19 GeV, we need |µ | 4 · 10 9 GeV. Note that µ is associated to the scale of the second doublet (a rather small representation not irreparably spoiling the unification of gauge couplings), but not to any other particle, since the 5s contain beside the doublets only the already considered triplets. The effective superpotential of the dimension 5 operators, which are relevant for proton and neutron decay, reads whereâ,b,ĉ are SU(3) C indices with 123 = 123 = 1, and the contraction of SU(2) L indices is given by Ψ · Φ = ab Ψ a Φ b with 12 = 1. The quasi-Yukawa operators and the mass term for the (anti-)triplets are given by After integrating out the (anti-)triplets, the effective baryon and lepton number violating operators are obtained: where the dots indicate terms which are not relevant for nucleon decay. According to Eq. (B.1), the dimension 5 operators at the GUT scale are then calculated as follows:

B.1.2 Running of the dimension 5 operators
The running of the dimension 5 operators from the GUT scale to the SUSY scale is calculated in the framework of the MSSM. The β-functions of the dimension 5 operators in the MSSM at 1-loop in the DR scheme, by using the left-right convention for the Yukawa matrices, are given by The β-functions are computed by means of the general formulas in [62] for the 1-loop anomalous dimension matrix γ m i , and the non-renormalization theorem of the superpotential [63][64][65]. We use the GUT normalization for the U(1) Y gauge coupling g 1 ; the SM normalization can be recovered by g SM 1 = 3/5 g 1 . The RGEs are then given by where x ∈ {C ijkl 5L , C ijkl 5R } and µ is the renormalization scale.

B.1.3 Dressing of the dimension 5 operators
The dimension 6 four-fermion operators are built at 1-loop level by dressing the dimension 5 operators with gluino, chargino and neutralino exchange diagrams, which takes place at the SUSY scale. An example diagram is shown in Figure 6. The complete set of formulas to calculate the dressing of the dimension 5 operators is stated in [12], in particular Eqs.  Beside the dimension 5 operators, the calculation of the resulting dimension 6 operators involves the rotation matrices of the sfermions, charginos and neutralinos, as well as the mass eigenvalues of the sparticles. It is convenient to use the standardized SLHA [15] and SLHA2 [67] conventions to specify these quantities. In Section B.2 the translation of the conventions used in [12] into the SLHA conventions is stated.
Apart from integrating out the sparticles at the SUSY scale, we also switch to the EW symmetry broken phase of the SM. In the following u i , d i and e i (i ∈ {1, 2, 3}) represent the mass eigenstates of the up-type quarks, the down-type quarks and the charged leptons. In contrast, the ν i represent the neutrinos in the interaction (flavor) basis. The additional label L or R indicates whether a state is a left-or a right-handed Weyl spinor. The dimension 6 operators C (at the Lagrangian level), which are relevant for the calculation of proton decay, include only mass eigenstates of the fermions which are lighter than the nucleons and are given by Since only the lightest up-type quark u 1 is lighter than the proton and the neutron, only this state is present in Eq. (B.9) and the label i is neglected. 9

B.1.4 Running of the dimension 6 operators
The running of the dimension 6 operators from the SUSY scale to the mass scale of the nucleons is calculated in the framework of the SM. Following [68], the main contribution to the running of the dimension 6 operators comes from the strong gauge coupling g 3 . At 1-loop the β-functions in the MS scheme of the dimension 6 operators C from Eq. (B.9), written in the general form C ijkl , and of the gauge coupling g 3 are given by where N F is the number of quarks whose masses are below the renormalization scale µ. The RGEs have the form as specified in Eq. (B.8). The running of the dimension 6 operators just corresponds to an overall scaling factor, which is also referred to as the long range effect on the effective operators C.

B.1.5 Partial decay widths
The formula for the calculation of the partial decay widths of a nucleon B i with a meson M j and a lepton l k in the final state is taken from [12] and is given by where m i and m j are the masses of the nucleon and meson, respectively. Moreover, f π = 0.130 GeV (cf. PDG [37]) is the pion decay constant. The amplitudes A ijk L and A ijk R for different decay channels are calculated by using the dimension 6 operators at the nucleon mass scale, and are specified in Table 1 in [12]. In that table, the nucleon mass m N is taken as the average of the proton mass m p and the neutron mass m n , and the baryon mass m B is given by the average of m Λ and m Σ . The additional constants with the corresponding (approximate) values are the constants from hyperon decay F ≈ 0.463, D ≈ 0.804 (cf. [7,69]), and the parameters for the proton decay matrix element α ≈ 0.0090 GeV 3 , β ≈ 0.0096 GeV 3 (cf. [7,70]).

B.2 Matching with SLHA conventions
In the following the relation between the conventions in [12] (later referred to as GN), which are used in the dressing of the dimension 5 operators, and the standardized set of conventions in SLHA [15] and SLHA2 [67] are specified.
• The Yukawa matrices are related as follows: • Both GN and SLHA use the SM normalization for the gauge coupling g 1 of U(1) Y , and they also use the same convention for the µ-term.
• In the EW symmetry broken phase the mass matrices of the sfermions are written in a basis where the sfermions are aligned with their SM superpartners.
- -SLHA uses the basis where all fermion mass matrices are diagonal, which corresponds to the super-CKM/PMNS basis.
The relation between the sfermion, the chargino and the neutralino mass matrices in GN and SLHA are then the following: where each entry represents a 3 × 3-block. In contrast, the mass matrix M 2 ν ≡ M 2 LL of the sneutrinos is only 3 × 3-dimensional, since only left-handed neutrinos are present.
Moreover, the 2 × 2-dimensional mass matrix M C of the charginos is written in the basis ( W ± , H ± u ), and for the 4 × 4-dimensional mass matrix M N of the neutralinos the basis ( B, W 0 , H 0 d , H 0 u ) is used. • The definitions and relations of the sparticle rotation matrices in GN and SLHA are given by: where "diag" indicates a diagonal matrix with real, positive entries. The definitions for the sfermion rotation matrices are taken from GN: Eq. (A.1.10) and SLHA2: Eqs. (28)-(31), whereas the definitions for the chargino and neutralino rotation matrices are stated in GN: Eq. (A.1.12) and SLHA: Eq. (12) and (15).

C User's guide to Mathematica packages SusyTCProton and ProtonDecay
This section provides a guide for the usage of the Mathematica packages SusyTCProton.m and ProtonDecay.m. Visit the web page http://particlesandcosmology.unibas.ch/downloads/ protondecay.html to obtain the packages. The zip file ProtonDecay.zip there contains the README.md file with installation instructions, the .m package file and a directory containing some example notebooks for evaluation. For SusyTCProton, follow the link on the web page specified above. Mathematica package requirements: REAP for SusyTCProton, none for ProtonDecay. The documentation below uses the typewriter font to specify functions and variables. In particular, the dimension 5 operators C ijkl 5L and C ijkl 5R from Eq. (B.1) are labeled as CL5qlqq and CR5udeu, respectively.

C.1 SusyTCProton.m
The package SusyTCProton.m has the same functionality as the package SusyTC.m from [13], with the following extensions/modifications:

• RGEGetSolution:
This function has the same functionality as described in [13], where in addition the values of the dimension 5 operators RGECL5qlqq and RGECR5udeu can be output at a given energy scale.
The output of the function consists of a list with elements of the form {label,decayWidths}, where label specifies the decay channel (e.g. p → π 0 e + k ). Furthermore, decayWidths is the list of the corresponding decay widths in GeV at the nucleon mass scale, where k ∈ {1, 2} if the lepton in the final state is a charged lepton, and k ∈ {1, 2, 3} if the lepton is a neutrino. Note that the charged leptons are considered in the mass eigenbasis, whereas the neutrinos are considered in the interaction (flavor) basis.