Flavour Structure of GUTs and Uncertainties in Proton Lifetime Estimates

We study the flavour aspects of proton lifetime estimates in simple Grand unified models paying particular attention to their inherent fragility due to the notorious lack of control of some of the key parameters governing the relevant hard-process amplitudes. Among these, the theoretical uncertainties in the flavour structure of the baryon and lepton number violating charged currents due to the potential higher-order effects afflicting the matching of the underlying Yukawa couplings to the low-energy data often play a prominent role. Focusing on the minimal variants of the most popular unified models we study the potential instabilities of the corresponding proton lifetime estimates based on the renormalizable-level Yukawa fits with respect to the Planck-scale induced flavour effects. In particular, we perform a detailed numerical analysis of all minimal $SO(10)$ Yukawa sector fits available in the literature and show that the proton lifetime estimates based on these inputs exhibit a high degree of robustness with respect to moderate-size perturbations, well within the expected"improvement window"of the upcoming proton decay searches.


I. INTRODUCTION
Mortality of protons is one of the most prominent smoking-gun signals of the idea that strong and electroweak interactions may be just different facets of a unified gauge dynamics at a super-large (YeV) scale. Since its conception in mid 1970's [1] there has been a number of attempts to estimate proton lifetime at vastly different levels of accuracy characterized, namely, by the steadily improving quality of the input data, progress in the field theory calculation techniques, better understanding of the hadronic matrix elements and so on.
With the upcoming generation of dedicated experimental searches planned with the Hyper-K and/or DUNE facilities [2,3], which should be able to push the current lower limits (e.g., τ p > 1.4 × 10 34 y in the "golden" p → π 0 e + channel [4]) by as much as one order of magnitude, the importance of a good quality prediction becomes particularly pronounced.
From this perspective, the current status of the theory affairs is far from satisfactory. Barring the nonperturbative nature of the hadronic layer, even at the level of the underlying "hard" processes, i.e., with amplitudes featuring quarks rather than hadrons as initial and final states, good quality calculations turn out to be an endeavour of enormous complexity. Indeed, the simple structure of the basic baryon and lepton number violating (BLNV) vector currents 1 even in the minimal * helena.kolesova@uis.no † malinsky@ipnp.troja.mff.cuni.cz 1 In this study we will focus predominantly on the vector-bosonmediated amplitudes as those mediated by the colored scalars are often sub-leading due to the usual suppression of their couplings to the first-generation quarks and leptons.
Georgi-Glashow SU (5) model [1], as simple as it reads, encompasses a great deal of arbitrariness (in the mass of the leptoquark X µ , to be identified with the GUT scale M G , and, in particular, in the flavour structure emerging when these currents are recast in the quark and lepton mass basis) that may be only partially reflected in the currently accessible low-energy observables.
Concerning the relative impact of uncertainties in these basic parameters on the proton lifetime estimates, the most critical of these is the value of M G which is determined from the requirement of a proper coalescence of the three Standard Model (SM) gauge couplings at (about) that scale. To this end, note that the logarithmic nature of the gauge running makes even a small error in the low-energy boundary (or high-scale matching) conditions propagate into M G exponentially. This, in turn, calls for 2 higher-loop account of the running effects including the appropriate-level threshold corrections both at M Z as well as at M G (and other intermediate scales, if present); needless to say, this is a highly technically demanding task in practice.
Second in the row is the high degree of uncertainty in the flavour structure of the baryon and lepton number violating (BLNV) charged currents which is namely due to the generic lack of low-energy access to the currentto-mass-basis rotations in the sector of right-handed fermions (as the CKM and PMNS matrices are combinations of the left-handed ones only). If no extra information (such as, e.g., symmetry features of some of the fermionic mass matrices) is available, the total freedom in these unitary transformations is usually enough to spread the outcome of the proton lifetime calculation over many orders of magnitude.
In this respect, it is remarkable that the classical showstopper of the past, namely, the uncertainties in the hadronic matrix elements, have recently got tamed to such a degree (with typical errors pulled down to few tens of percent) that, nowadays, they can be safely placed as only third in the row, see, e.g., [5] and references therein.
With this basic hierarchy at hand, one can perform a simple classification of the robustness of the most commonly followed strategies in predicting proton lifetime: 1) The first attempt usually consists in the renormalization group (RG) analysis of the gauge unification constraints which provides information about M G but often ignores the flavour structure of the BLNV vector currents, typically because the scalar sector of the model is not fully fixed or analyzed. Hence, the uncertainties of thus obtained proton lifetime estimates are generally huge, stretching over many orders of magnitude. This, however, to a large extent hinders the prospects of discrimination among different scenarios. 2) Sometimes, a great deal of information may be derived from the symmetry features of the effective fermion mass matrices even without performing their detailed fit (usually very demanding), see, e.g., [6][7][8]. In specific scenarios like, e.g., in the minimal realistic SU(5) models, this may be enough to draw rather accurate conclusions about at least some of the partial decay widths (though often not for the "golden" channel p → π 0 e + ); the potential to discriminate among such models is obviously much higher then.
3) The ultimate achievement would be clearly a full-fledged combined analysis of the running together with a detailed Yukawa sector fit. This, however, is very difficult in practice and only very few such attempts have been undertaken in the literature, see, e.g., [9].
Nevertheless, even in the most favourable situation of case 3) above there is often an extra source of large and essentially irreducible uncertainties plaguing any proton lifetime estimate obtained in the realm of the simplest renormalizable models, namely, the effects of the higher dimensional effective operators, especially those including the scalar field(s) (to be denoted S l ) responsible for the GUT-scale symmetry breaking, i.e., the ones with the vacuum expectation values (VEVs) of the order of M G .
At first glance, there is a number of such structures to be considered at the d = 5 level (with the ordering reflecting their expected "nuisance" power), e.g., In the formulae above 3 X µν stands for the gauge field tensor, f i denote matter fermions, Φ is a generic scalar field, H k are scalars over which the SM Higgs doublet is spanned, and κ n denote (generally unknown) O(1) couplings; for the sake of simplicity the spinorial structure has been suppressed. It is important to notice that not all of these are, however, independent structures from the low-energy effective theory point of view: O 3 and O 4 may be removed from the effective operator basis by use of equations of motion and/or by integrations by parts, see, e.g., [10]. Hence, in what follows, we shall focus entirely on the O 1 and O 2 types of the d = 5 structures. Despite that, in the broken phase, they usually affect their renormalizable-level couterparts (namely, the gauge-kinetic forms and the Yukawa couplings) at only a relatively small -order 1% -level (given by the typical ratio of S ∼ M G ∼ 10 16 GeV and the Planck scale M P l ∼ 10 18 GeV) they can have truly devastating consequences for the robustness of the renormalizable-level results 4 , see, e.g., [11]. To this end, the most dangerous is O 1 which has been studied thoroughly in many works, see e.g. [12,13]. Its main effect, i.e., an inhomogeneous shift in the high-scale gauge matching conditions, can inflict a significant shift in the exponent of the functional dependence of the M G /M Z ratio which, even for 1% shifts in the matching, may change M G by as much as an order of magnitude and, hence, alter τ p by several orders.
Concerning O 2 , there are many studies in the literature (like, e.g., [14][15][16]) in which non-renormalizable contributions to the Yukawa couplings have been added to an originally renormalizable Lagrangian on purpose, usually with the aim to save a renormalizable model suffering from a badly non-realistic Yukawa sector (like in the minimal SU (5) model [1]). Needless to say, this approach is orthogonal to the line of thoughts we want to pursue here, namely, focusing on the potential impact of a-priori unknown Planck-suppressed operators on the existing renormalizable-level predictions, hence, testing their overall robustness.
In this study, we concentrate on the stability of the tree-level gauge-boson-mediated contributions to the proton decay in the simplest renormalizable unified models based on the SO(10) [17], flipped-SU (5) [18][19][20] and SU (5) [1] gauge groups with respect to several types of uncertainties, either due to the lack of any (or part of the) information about their flavour structure, or due to the presence of only mildly suppressed (up to order 1%) Planck-scale induced operators of the O 2 type above.
To this end, we first (in Section II) recapitulate the generic analytic observations made by Dorsner and Fileviez-Perez [6][7][8] in the realm of the simplest SU (5), flipped-SU (5) and SO(10) scenarios and complement them with an explicit numerical simulation of the relevant formulae revealing, e.g., an extra room for large cancellation effects in SO(10) GUTs; we will show how these can in some cases boost the uncertainties beyond naive expectation. In Section III, we consider a few specific types of renormalizable-level proton lifetime estimates and assess their generic robustness with respect to the effects inflicted by the possible presence of the M P l -induced Yukawa-altering d = 5 non-renormalizable operators. To this end, we especially focus on a thorough numerical analysis of the stability of the flavour structure of the BLNV currents corresponding to a variety of existing renormalizable-level Yukawa-sector fits [21][22][23][24][25] in the minimal SO(10) and its variants.

A. Partial decay widths
Focusing on tree-level amplitudes mediated by heavy vector bosons arising from terms like (1) in the Lagrangian the relevant partial proton decay widths can be written as 5 [26] where incoherent summation over the neutrino flavours is performed, since the neutrinos in the final state are not detected; for similar reasons it is also summed over the chirality of the charged leptons in the final state. 5 We assume here that neutrinos are Majorana and some form of a seesaw mechanism is in operation, hence ν R is too heavy to be produced in proton decay. Should neutrinos be Dirac the sensitivity of the proton widths is less pronounced, see e.g. [6].
The definition of the flavour independent prefactors C π , C η , C K and B 1,2,3 is postponed to Appendix A, let us focus here on the flavour structure of the partial widths determined by the c-amplitudes where k i = g G /( √ 2M i ) with g G denoting the universal gauge coupling at the GUT scale and M 1,2 encoding the masses of the heavy vectors with the SU (3) c × SU (2) L × U (1) Y quantum numbers (3, 2, 5/6) and (3, 2, −1/6), respectively. In the flipped-SU (5) scenario only the latter is present and, hence, k 1 = 0; similarly, k 2 = 0 in case of the ordinary SU (5). In the SO(10) GUTs both k 1,2 are non-zero. The unitary rotations entering the coefficients (8)-(10) are defined as Up to a possible multiplication by phase factors the LH rotations in (8)-(10) are correlated to the physical CKM and PMNS matrices via where the tilded quantities correspond to their "raw" form, i.e., the form before the freedom in the phase redefinition of the fermionic fields 6 has been exploited. Note that, in general, these are also the only experimental constraints on the rotations in (11) one has; without extra information about the flavour structure of a particular model U C , D C and E C (assuming Majorana ν's) are completely free.

B. Model independent constraints
With this information at hand, number of semianalytic observations about the relative sizes of the uncertainties plaguing the partial widths (3)-(7) due to the lack of grip on most of the flavour structures therein has been made [6,7] even without any extra model-dependent assumptions on the shape of the mixing matrices in (11). In particular, the question whether the total proton decay width 7 can be zeroed out, i.e., whether the proton decay can be "rotated away" was addressed.
On the contrary, in the standard SU (5) GUTs where k 1 = 0, k 2 = 0 the non-zero value of |(V CKM ) 13 | element forbids to rotate the proton decay away [7]; however, the small size of this parameter admits up to about O(10 −3 ) suppression of the amplitudes (8)- (10) and, hence, up to some O(10 −6 ) suppression of the total decay width (see FIG. 1). Consequently, the proton decay can be "hidden" from the current experiments even if the unification scale would be as low as 10 14 GeV [6].
In case of SO(10) unifications with k 1 = 0, k 2 = 0 one would naively expect similar lower bounds on the amplitudes (8)-(10) which, in turn, might suggest the same O(10 −6 ) maximum suppression of the total proton decay width. However, with both k 1 and k 2 at play, destructive interference effects may sometimes occur in all the coefficients (8)-(10) which would make the total proton decay width even smaller than that, see FIG. 1 for a numerical simulation (with M 1 = M 2 assumed for simplicity 8 ).

C. Minimal renormalizable settings
In scenarios in which the scalar sector is specified, extra correlations among the flavour rotations (11) are often in operation. If, for instance, some of the Yukawa matrices happen to be symmetric the RH and LH rotations are strongly correlated and significant simplifications may occur in formulae (8)- (10).
Namely, in the original Georgi-Glashow SU (5) model [1], the symmetry of the fermionic 10 ⊗ 10 bilinear in 7 Here it is implicitly assumed that the rates into final states with higher spin mesons will be suppressed; moreover, their flavour structure is essentially the same as for the spin zero modes (3)- (7) and, hence, qualitative changes of the results are not expected. 8 This is justified by the fact that for M 1 M 2 or M 1 M 2 one of the two previously discussed cases is effectively recovered. the flavour space implies that the up-quark Yukawa matrix is symmetric and, hence, U = U C . In this case the partial widths with (unidentified) neutrinos in the final state become entirely driven by the CKM matrix elements [27] and the uncertainty in the total proton lifetime shrinks considerably as depicted in FIG. 1. On the other hand, the basic SU (5) flavour structure considered above should be extended in order to deal with the down-type-quark-charged-lepton degeneracy issues (and other notorious problems of the simplest SU (5) unifications concerning, e.g., the unification of gauge couplings or non-zero neutrino masses). New fields (such as 45dimensional scalar representation [28,29]) and/or higherdimensional operators [14,16,30] are usually employed for that sake. In either case the exact symmetry of Y u is lifted. For this reason, the red bar for SU (5) unifications in FIG. 1 is to be taken as purely illustrative.
In the flipped SU (5) scenario, the RH quark field d C is swapped with u C and, hence, it is the down-type Yukawa that gets symmetric in the minimal settings which implies D = D C . This leads [27] to very simple relations and, consequently, to a significant reduction of the uncertainty in the proton lifetime estimates, see FIG. 1. Let us emphasize that in case of flipped SU (5) the constraint D = D C is satisfied also by fully realistic models including all necessary ingredients like, e.g., non-zero neutrino masses [31,32]. In case of the renormalizable SO(10) unifications, the minimal potentially realistic choice of the scalar fields shaping the Yukawa sector corresponds to a 10dimensional vector and a 126-dimensional 5-index antisymmetric self-dual tensor. Both these yield symmetric Yukawa couplings and, thus, all RH rotations in (11) are strongly correlated with the LH ones. This leads to [8] Again, the uncertainty in the total proton lifetime shrinks enormously, see FIG. 1.

D. Two-body p-decay amplitudes with a charged lepton in the final state
Unlike for the (anti)neutrino channels above the amplitudes of the two-body partial proton decay widths with a charged lepton in the final state are generally driven by non-trivial combinations of the mass-diagonalization matrices (11) with only an indirect 9 connection to the low-energy flavour observables.
Hence, in order to get any theoretical grip on these channels, one must resort to a specific model and construct a detailed map of all possible phenomenologycompatible flavour patterns, i.e., a complete set {F, F C } with F = U, D, E, N defined in (11). Technically, this information can be obtained from a thorough analysis of the fits of the Yukawa structure underlying the quark and lepton mass matrices M f . For each such setting the relevant amplitudes can then be fully reconstructed and, if desired, extremized over the entire set of such configurations.
This, however, is a highly non-linear game and, thus, in most cases [21][22][23][24][25], the authors resort to the renormalizable-level approximation in which one typically deals with a limited number n of independent Yukawa matrices Y k of the d = 4 Yukawa lagrangian entering the renormalizable-level matching conditions for the effective low-energy couplingsỸ f in the form (no summation over f = u, d, , ν, generation indices suppressed): Here v denotes the electroweak VEV serving merely as a normalization factor, v k f stand for the projections of the SM Higgs VEV onto the underlying-theory doublets H k that couple to the fermionic bilinears at the d = 4 level and c k f cover all remaining constant O(1) numerical factors (Clebsches, symmetry coefficients and so on). It is clear that if D 4 , the set of indices corresponding to the doublets relevant at the d = 4 level, is small (i.e., if the number of such doublets is less than 4) the fits of Y f 's in terms 10 of Y i 's may be quite non-trivial, {F, F C } strongly constrained and, thus, the theory predictive.

III. D=5 PLANCK-SCALE FLAVOUR EFFECTS
In reality, however, the renormalizable-level fits may be incomplete because the effective Yukawa matrices may be affected by physics at the Planck-scale M P l which, in the effective theory picture, may enter the game by means of non-renormalizable d > 4 operators. For instance, the presence of the d = 5 operators of the Yukawa type [class O 2 in the list (2)] inflicts additional shifts in the relevant matching conditions between the effective (running) lowenergy Yukawa couplings Y f and the underlying GUTtheory couplings in the form where ∆Y f may be schematically written as (again, no summation over f and no generation indices) Here the meaning of v k f is the same like above, D 5 is the set of relevant doublet indices that is summed over (note, however, that D 5 is not 11 necessarily the same as D 4 in (15)), S 5 indexes the SM singlets with GUT symmetry breaking VEVs, κ kl are the coefficients of the relevant d = 5 operators as in (2) and, as before, c kl f cover all the remaining numerical factors.
Comparing (15) with (17) and assuming that the coefficients c kl f and κ kl are at most O(1) one may expect that the Planck-induced contributions ∆Y f to the full effective Yukawa couplingỸ f in (16) should typically come with an additional M G /M P l ∼ O(10 −2 ) factor. 10 The fundamental doublet projections v i f are, in principle, calculable functions of the scalar potential parameters and, usually, turn out to be correlated among themselves. An extreme example of this is the situation in the minimal SUSY SO(10) model [33] which was eventually discarded [34,35] just due to such correlations. 11 At d = 5 there may be contributions in (17) from scalar multiplets containing SM doublets that would not be present in (15) and vice versa. As an example consider the 126 and 210 scalars in SO(10) -the latter irrep can, indeed, couple to the fermionic bilinears only at the d = 5 level.
Although such corrections may naively appear to be small, they may make the matrices diagonalising the completeỸ f 's (which we shall from now on denote by {F ,F C }) significantly different from those obtained at the d = 4 level (denoted by {F, F C }). The point is that the entries of the Y f matrices (and, most importantly, their eigenvalues) are often smaller than O(M G /M P l ) ∼ O(10 −2 ); hence, even an O(10 −2 ) correction can be enough to change the shape of the "small" part of the Yukawa matrix completely (we further elaborate on the formal aspects of this in Appendix B).
On the other hand, even though the two sets of rotation matrices {F, F C } and {F ,F C } may look dramatically different, the resulting partial proton decay widths may still be rather similar since the relevant formulae (3)-(7) depend only on their specific products and, as we shall see, in some cases there may be reasons to expect significant cancellations of such d > 4 effects.
A. Robustness of the renormalizable-level p-decay estimates with respect to d > 4 effects At first glance, it may seem rather hopeless to attempt to say anything general enough to be interesting about the possible differences of the two sets of matrices {F, F C } and {F ,F C } -i) either one fits the effective Yukawa matrices in a renormalizable model and then has little or no grip onto a typically yet larger set of the higher-order operators, or ii) one includes nonrenormalizable operators into the game right away (because it may be necessary to do so otherwise no consistent parameter-space points may be found at all, see e.g. [30,[36][37][38]) and then, naturally, never asks about the renormalizable case because it makes little sense.
The main scope of this work is to argue that the situation corresponding to the case 12 i) above may be slightly more subtle and, in fact, under some circumstances, one may say something sensible about the robustness of the p-decay estimates based on the renormalizable-level Yukawa sector fits even without a detailed knowledge of the structure of the higher-dimensional contributions therein.

First look: The problem in full generality
Let us start with assuming for the moment an ideal world in which we have enough computing power to generate all (with perhaps some given granularity in practice) fits of the effective Yukawa matricesỸ f in terms of the underlying renormalizable-level Yukawa couplings subject to sum-rules dictated by the unified model under scrutiny and perhaps even more power to repeat the same exercise for the more complicated d > 4 case.
The former, in other words, amounts to getting first all possible structures of Y k 's in (15) associated to all attainable configurations of v l 's which, in the d = 4 case, yield theỸ f 's that (after the necessary renormalization group evolution to our energies) encode the desired spectra of the SM fermions together with their mixing in the charged-current interactions (aka the CKM and PMNS matrices). With these at hand one would then easily derive the complete set of the possible {F, F C } matrices which shall be eventually used to estimate the proton lifetime and, in particular, the associated theoretical uncertainty corresponding to the fact that the low energy data can not pinpoint the "true" solution among all these possibilities. In the second step one repeats the same exercise with just a little bit more of freedom due to the presence of the extra couplings corresponding to the d > 4 operators and derives all possible {F ,F C } associated to these "extended" fits together with the relevant proton lifetime estimates.
Given this it is immediately clear that: • The set F of all thus obtained "realistic" {F, F C } is a subset of the setF of all "realistic" {F ,F C }; • Without any specific constraint on the size of the higher-order operators the set of the "realistic" {F ,F C } may be so large that one effectively looses any grip on the vector leptoquark interactions and, subsequently, the theoretical uncertainties of the proton lifetime estimates within such scenarios rocket (and, hence, exhibit the behaviour depicted by the blue bars in FIG. 1).
The point we will try to make is that, in some cases, even a simple extra assumption such as an additional O(10 −2 ) suppression associated to each subsequent step on the effective operator ladder, in conjunction with specific features of the renormalizable-level fits such as their symmetry in the generation space, may be enough to cor-relateF to F to such a degree that the proton lifetime estimates obtained within the humble d = 4 approach may actually represent a very good approximation to the "true" (i.e., full theory) predictions.

The trick: Small perturbations and continuity
Needless to say, the programme sketched in the previous part (i.e., obtaining the sets F andF -both complete -and comparing the spans of the associated proton lifetime estimates in order to asses the robustness of those based only on F) is intractable 13 . However, for small 13 Besides intractability it does not even make sense to do that because with the completeF at hand nobody would care about F anymore.
|∆Y f | 10 −2 (for all f 's) the task to learn something aboutF can be accomplished even without embarking on its full determination by assuming continuity in the change of the fitted renormalizable-level Yukawa couplings as functions of the size of the non-renormalizable contributions.
Technically, what we have in mind is that the fullF set obtained upon fitting the complete non-renormalizable structure Y f + ∆Y f with no assumptions made on Y f and ∆Y f (besides the smallness of the latter, see (16)), will be essentially 14 the same as the setF obtained by fitting a slight variation of the original formula, namely is assumed to run over the set of all good fits of the renormalizable-level case only.
The point is that due to continuity assumption the choice of specific Y f 's instead of a fully general Y f inflicts only a small change in the structure from which F would have been derived. Indeed, this change can be modelled by a mere reshuffling of the set of the Planck-scale-induced corrections that would have to be summed over anyway to get the completeF (indeed, and ∆Y f → ∆Y f replacements qualitatively correspond to a mere choice of a specific reference element in the set of all possible d > 4 contributions 15 ! Needless to say, this leads to an enormous simplification of the general problem described in Sect. III A 1 and, as such, it represents the central point of this study (and, in fact, the very key to its practical feasibility). Hence, we may simplify our life in mapping theF set by reformulating the general exercise described in Sect. III A 1 into a much more tractable one of employing only the specific shapes of the renormalizable-level Yukawa couplings corresponding to the renormalizable-level fits whilst keeping fully general (i.e., unspecified but small) only the d > 4 contributions. Thus obtainedF should be essentially the same asF.
In conclusion, what one should do in practice is to take all possible renormalizable-level fits Y f of the effective Yukawa structure of the given model; then, for each such Y f , construct the sums , take all cases in whichỸ f 's happen to give the right 14 One may argue that in this way we are mapping theF set corresponding to 2ε rather than the original ε size of the d > 4 perturbations. This is true but, at the same time, we do not really care because in the semi-qualitative discussion to follow this makes no significant difference. 15 One may wonder if the setF is not subject to an extra restriction compared toF due to the fact that Y f andỸ f = Y f + ∆Y f must have the same generalized eigenvalues in order to yield correct fermion masses. However, the Lemma 3 in Appendix B suggests that arranging for the correct generalized eigenvalues (up to the order of O(ε 2 )) does not restrict the set {F ,F C } at all.
SM fermion masses and mixings (as Y f do by construction) and, eventually, derive and save the corresponding {F ,F C }'s. These will be, subsequently, used as inputs of a refined p-decay analysis in the models of interest.

Further comments
There are perhaps few more comments worth making here: First, it may well be the case that the complete set of all possible renormalizable-level fits that the trick above relies on may not be fully available as its determination represents a formidable task on its own. To this end, in what follows we shall do what we can, i.e., we shall take a look onto just a specific (and small) set of popular and simple enough scenarios and, within these, confine ourselves to all available sets of Y f 's that may be found in the corresponding literature. In this sense, the results presented in the next section may not be completely general; nevertheless, they are not useless as they admit to estimate the robustness of at least the existing proton lifetime calculations based on the renormalizablelevel flavour fits.
Second, the entries of the individual ∆Y f 's of eq. (17) may be in general further restricted due to the extra symmetries of the underlying effective operators in the generation space. Out of the restrictions of this kind, the case with both Y f and ∆Y f symmetric for certain f 's is of most interest since, in such a case, also the effective Yukawa couplings of Eq. (16) inherit this symmetry. Consequently, the set {F ,F C } has the simplified structure described in Section II C and the total proton lifetime uncertainties are contained within the red bars in FIG. 1. Remarkably, this is the case for the SO(10) GUT featuring the 10-and 126-dimensional Yukawa active scalars with the first-stage gauge symmetry breaking driven by the 54-dimensional scalar 16 (see, e.g., [9] for a recent study including the Yukawa sector fits). Such situation is, however, rather exceptional. For instance, if the SO(10) gauge group was broken by 45 instead, ∆Y f 's would contain also an antisymmetric part due to the presence of the 120-dimensional representation in the product 45 ⊗ 10. As already mentioned in Section II C, the situation is similar in the case of the simplest SU (5) unifications, where the d = 5 non-renormalizable operators destroy the symmetry of the up-type quark Yukawa matrix (see, e.g., [30]). Thus, in what follows, we decide not to impose any extra generation-space symmetries onto the ∆Y f 's of eq. (17).
Third, the ∆Y f 's of eq. (17) may be, in principle, further correlated across different flavours (i.e., f 's) due to their common origin from a potentially limited set of available d > 4 effective operators. Such correlations are, however, strongly model-dependent; therefore, we choose to ignore such nuances in the current analysis. This means that our results may be viewed as corresponding to the most pessimistic situation and, in reality, the uncertainties of the proton decay estimates within specific scenarios may be smaller. If, on the other hand, the partial proton decay widths turn out to exhibit a certain degree of robustness with respect to the uncorrelated O(M G /M P l ) perturbations, the same behaviour should be reflected also in the real, i.e., more constrained case.
In this respect, the approach of imposing no extra constraints onto the shapes of ∆Y f 's (besides their smallness) is perhaps the only strategy which can, on one hand, reflect the specifics of the underlying renormalizable-level fits and, at the same time, save thus obtained results from any further model-dependent assumptions.

The numerical approach
Let us now describe the technical aspects of the numerical analysis of formulae (18).
Since both Y f andỸ f are assumed to yield the same physical masses (see Section III A 2), one finds The perturbation ∆Y f can be, hence, expressed in terms of the varied rotation matricesF andF C . Note that it is more convenient to search through the space of the unitary matricesF andF C instead of the perturbations ∆Y f since then the constraints regarding the CKM and PMNS matrices (12)-(13) can be easily implemented. Consequently, our strategy is to exploit the set of all possible shapes ofF andF C satisfying (12)- (13) and to check only subsequently whether the resulting perturbation of the Yukawa matrix is as small as required, i.e., whether is satisfied for f = u, d, e. 17 17 We constrain the matrix N in (11) only by the PMNS matrix relation (13) since Yν is usually computed from different Yukawa couplings according to some type of seesaw mechanism, and the corresponding constraints on ∆Yν would be more complicated and model dependent. This follows our strategy to consider the "worst case" scenario; in reality, the true uncertainties in the corresponding proton lifetime may be more constrained. Moreover, since the channels with neutrinos in the final state are always incoherently summed over in (3)- (7), we do not expect that further constraining N would have any significant effect.
Finally, let us comment on the choice of the renormalizable fits that serve as Y f 's in (18). In case of the SO(10) unifications, exact shapes of the fitted Yukawa matrices are available in the literature [21][22][23][24][25]; these served as inputs for Y f 's in our numerical analysis.
On the other hand, to our best knowledge, neither for the SU (5) nor the flipped SU (5) models any reasonably exhaustive classification of working fits of Yukawa sectors of the minimal potentially realistic and renormalizable scenarios is available in the literature. Nevertheless, it may still be interesting to check how much the set {F ,F C } varies from {F , F C } corresponding to a set of essentially random choices of Y f 's which, however, are still assumed to respects at least the basic symmetry properties inherent to the minimal models 18 , see Section II C. In case of the flipped SU (5), symmetric Y d was always assumed for the starting point while in case of the ordinary SU (5) both symmmetric and nonsymmetric versions of Y u were checked since, in realistic models, the latter options is usually realized as explained in Section II C.

Simplest SO(10) GUTs.
Given their relatively rigid Yukawa structure, the SO(10) GUTs provide an ideal setting for us here since a decent number of renormalizable-level Yukawa fits available in the literature can serve as a starting point for our numerical analysis. In total, 8 different Yukawa sector fits for non-supersymmetric SO(10) models available in [21][22][23][24][25] was studied, including both the cases when 10 H ⊕ 126 H 'Yukawa-active' Higgs fields have been taken into account, hence, the renormalizable-level mass matrices were symmetric, and when also the antisymmetric contribution due to the presence of 120 H was added; no significant qualitative differences were observed. As an example, the scan over the space of possible Yukawa matrix perturbations based on the fit obtained in [22] considering 10 H ⊕126 H Higgs sector and the normal neutrino mass hierarchy is presented in plots in FIGs. 2 -4. When computing the partial proton decay widths, M 1 = M 2 ∼ 0.3 × 10 16 GeV was fixed for which the current SK bound on the proton lifetime in the "golden" channel p → π 0 e + is just saturated. The uncertainty in the (inverse) partial proton decay widths or their sums was then plotted against the maximum size of the perturbations (19) |∆Y | ≡ max 18 In practice, the relations (11) were used, i.e., the diagonal part  . The renormalizable-level setting of the Yukawa matrices Y f corresponds to the fit explicitly given in [22] where the 10H ⊕ 126H Higgs content and the normal neutrino mass hierarchy are assumed. Similar behaviour with large spread of the partial proton lifetime values was obtained also for other individual decay channels and also when other fits in [22] served as the initial point.
The spread in the individual (inverse) partial proton decay widths was observed to be large for a wide range of |∆Y | 10 −2 , see FIG. 2 for the example of the "golden" channel p → π 0 e + . This, unfortunately, means that even with a fit to Yukawa sector at hand, robust predictions for the individual decay channels are in general impossible. The same behavior with large uncertainties was observed also if it was summed over all the partial widths with the charged leptons in the final state. On the other hand, the situation became much more favourable when the neutrino channels were summed over as shown in FIG. 3. Let us note that this behaviour can be understood recalling that it is summed over the neutrino species in the final state for the neutrino channels, whereas the production of the τ lepton is kinematically forbidden, hence, there is more room for "rotating away" the proton decay to the unobservable sector in the charged lepton case.
The particular robustness of the decay modes with neutrinos in the final states is subsequently reflected in the robustness of the total proton lifetime in SO(10) unifications, see FIG. 4. This is due to the fact that, in this scenario, these partial widths are never significantly suppressed with respect to those into charged leptons.

Flipped SU(5) unifications.
As explained above, the lack of dedicated fits of the Yukawa structure for this class of unification models lead us to choosing a random starting point when the stability of this scenario with respect to the Planck-scale corrections was examined. Remarkably, the qualitative log 10 |∆Y| The inverse of the sum of the partial p → X + ν decay widths for X = π, K in case of the SO(10) unifications (the same renormalizable-level point as in FIG. 2 was used, and again similar behaviour was observed also for other initial settings based on fits from [22]). Note, however, that the behaviour of the individual neutrino-final-state channels is similar to the situation for the golden channel (see FIG. 2). behaviour was independent of the starting point choice, hence, we believe that also the results regarding the flipped SU(5) unifications are worth presenting here.
As shown in FIG. 5, for |∆Y | only slightly below 10 −2 a significant instability of proton lifetime estimates based on the purely renormalizable structure was revealed. This supports the observation of [7] that the proton decay can be indeed rotated away in the flipped SU (5) scenarios although at the renormalizable level the predictions seem to be rather robust (see the formulas (14)).
On the other hand, since in the simplest flipped SU (5) models 19 the unified gauge group is broken by a scalar representation charged with respect to the U (1) X gauge group, one easily finds that the d = 5 Yukawa-affecting operators like O 2 in (2) are absent. Consequently, the first Planck-induced structures that may generate uncontrolled shifts in the underlying Yukawa couplings emerge only at the d = 6 level, hence, their size is expected to be of the order of 10 −4 . As can be seen in FIG. 5, for such |∆Y | the uncertainty in the total proton decay width becomes reasonably constrained, which means that, in the end, the proton lifetime estimates in the minimal flipped SU (5) scenarios are particularly robust.

SU(5) GUTs.
Similarly as in the case of flipped SU (5) unifications, also for ordinary SU (5) we had to rely on a random renormalizable ansatz for the Yukawa matrices Y f in our numerical analysis. Contrary to the flipped SU (5) case, however, the realistic renormalizable models do not feature any symmetric Yukawa matrices (see Section II C), hence, starting points with U = U C were assumed in general. For different initial settings of this type no common feature was observed -even the total proton decay width could be spread over several orders of magnitude for certain choices of the starting point, although, on the 19 Besides the original works [18][19][20] where the scalar sector is often not considered in detail, we have in mind, e.g., the fully realistic models including also non-zero neutrino masses like [31] or [32] other hand, for the settings with U being close to U C the uncertainty coming from the unknown Planck-scale contributions was much smaller. 20

C. Remarks
First, let us provide a hint on how the results presented in the plots above could be understood analytically. As stated in Lemma 2 in Appendix B, if the first n generalized eigenvalues of a matrix Y f are of the order of O(ε), then an O(ε) correction to Y f changes the n × n upper left blocks of the diagonalization matrices completely. This means that the uncertainty in the rotation matrices (11) and, hence, in the proton lifetime estimates qualitatively changes whenever the size of the perturbations crosses a generalized eigenvalue of Y f (which is proportional to one of the fermion masses). Since the largest entries of the down-quark and charged-lepton Yukawa matrices corresponding to the b and τ masses are around 5 × 10 −3 , for such values of |∆Y | there appears the first "step" in the FIGs. 3,4 (when perturbations with the size above this threshold are added to Y f 's the matrices D, D C , E, E C can be changed completely; for |∆Y | below this threshold only the upper left 2 × 2 corner of these matrices may be significantly varied). The other "step" in plots of FIGs. 3,4 corresponds to crossing the values of the Yukawa matrix entries corresponding to the c and µ masses.
Finally, a comment is worth concerning the SO(10) with the 10 ⊕ 126 ⊕ 54 scalar sector. As mentioned in Section III A 3 this scenario is exceptional since both the renormalizable Y f 's and the d = 5 Planck-induced corrections ∆Y f 's are symmetric, hence, the flavour structure of the partial proton widths with neutrinos in the final state remains fully determined even if the non-renormalizable terms are included. However, even with this extra information at hand, the decay channels with the charged leptons in the final state still exhibit a numerical behaviour similar to the case with general ∆Y f 's (see FIG. 2), i.e., the spread in these partial widths remains rather large.

IV. CONCLUSIONS AND OUTLOOK
In the current study, we have elaborated on the robustness of the gauge-boson-mediated contributions to the proton decay width in the simplest renormalizable unified models based on the SO(10) and SU (5) gauge groups with respect to several types of uncertainties, in particular those due to the presence of the Planck-scale induced operators altering the renormalizable-level Yukawa structure of specific models (such as O 2 in the list (2)). These perturbations, as small as they may seem in comparison to the typically huge effects inflicted by the notorious gauge-kinetic-form-changing d = 5 operators (i.e., O 1 in (2)), are still significant enough to trigger large changes in the mixing matrices governing the relevant baryon and lepton number violating currents and, hence, cripple in principle the credibility of any of the existing renormalizable-level proton lifetime estimates.
Remarkably enough, a thorough numerical analysis reveals vastly different levels of robustness of the relevant proton decay widths across different variants of the simplest SO(10) and SU (5) scenarios. Let us recapitulate the main observations that we managed to make here: Unfortunately, for all scrutinized models, the individual decay channels with charged leptons in the final state were found to be prone to significant destabilisation even for rather small Planck-induced perturbations. Typically, the inflicted theoretical uncertainties prevent, e.g., the "golden" channel p → π 0 e + from discriminating efficiently among different scenarios (even with a fit to the renormalizable Yukawa structure at hand) unless the overall suppression associated to the relevant d = 5 operators happens to be well under the 10 −2 level expected from the simple M G /M P l ratio (see FIG. 2).
Concerning the SO(10) scenarios, our numerical analysis reveals that, for all available renormalizable Yukawa fits within the minimal renormalizable models, the sum of the partial decay widths with neutrinos in the final state, and, consequently also the total proton lifetime turn out to be quite trustable from the flavour structure point of view even if the overall suppression factor associated to the Planck-scale effects is as large as 10 −2 (see FIGs. 3 and 4). This applies, namely, to the minimal potentially realistic scenario with the GUT-scale symmetry breaking triggered by the 45 scalar [39,40] which, besides this feature, exhibits a spectacular level of robustness with respect to the gauge-kinetic effects associated to the O 1 operator in the list (2). Hence, the leading Planck-scale induced theoretical uncertainties plaguing the renormalizable-level proton lifetime estimates in this scenario can be well within the "improvement windows" of the upcoming megaton-scale facilities such as Hyper-K or DUNE. It is also worth mentioning that the alternative scenario in which SO(10) is broken by 54 instead of 45 is yet more robust as far as the flavour strucure is concerned. On the other hand, it suffers from significant theoretical uncertainties in the overall BLNV scale determination which make it somewhat less attractive from the phenomenology point of view.
Concerning the SU (5)-based scenarios, the flipped SU (5) would be our primary choice as (in the minimal variants such as, e.g., [32]) it typically exhibits a higher degree of robustness of the flavour structure governing the proton lifetime calculations due to the generic ab- Let us complete the formulas for the partial proton decay widths (3)-(7) by the definition of the flavour independent prefactors. For the sake of continuity with the previous works we use the parametrization used in [26] based on the chiral Lagrangian: where m p , m η and m K denote the proton, η and kaon mass, respectively, m B is an average baryon mass (m B ≈ m Σ ≈ m Λ ), f π is the pion decay constant, |α|, D and F are the parameters of the chiral Lagrangian, and A L takes into account the renormalization from M Z to 1 GeV. On the other hand, the recent lattice computations [5] predict directly the individual matrix elements without the use of the chiral Lagrangian. These results can be, however, translated to the chiral Lagrangian parametrization (see, e.g., Appendix A of [5]) and the value of the parameter α is then inferred. 21 Let us stress that the possible change in the value of this multiplicative factor given by the improvement of the lattice computations does not affect our results qualitatively, the points in all the plots would be merely shifted in a uniform way.

Appendix B: Matrix diagonalization
Let a complex matrix Y be diagonalized by a biunitary transformation where Y d is a real non-negative diagonal matrix consisting of the so-called generalized eigenvalues of Y . Since the main issue of this paper is the sensitivity of the unitary matrices U, U C to the small perturbations in the matrix Y , let us mention here few mathematical results on this problem. As a first step, however, let us remind the reader about the way in which the matrices U and U C in (B1) are constructed. Since Y † Y is a hermitian matrix, it can be diagonalized as where D is a real non-negative diagonal matrix. The diagonal matrix in (B1) is then defined as Y d = √ D and if its entries are non-zero, then the unitary matrix U C can be defined as Let us note, however, that the matrices U and U C are not defined uniquely and the level of this ambiguity depends on the shape of Y d : i) For non-degenerate and non-zero diagonal entries of Y d , the ambiguity in U in (B2) amounts to U → U P with P being a diagonal unitary matrix. The matrix U C in (B3) is then accordingly transformed as n is allowed where the upper left corner of U n is formed by an n × n unitary block, U jj n = e iφj for j > n and U n ij = 0 otherwise.
then the definition (B3) of U C can not be applied and the relation has to be used instead. The ambiguity in the definition of the rotation matrices than reads 21 Let us note that only the matrix elements of the RL type like π 0 |(ud) R u L |p are relevant for the vector-boson-mediated proton decay, hence, only the parameter α enters the formulas (3)- .
with the same structure of U n , U Cn as in point ii). Here, however, the upper left n×n blocks of U n , U Cn are uncorrelated and the phases of the diagonal entries U jj n , U jj Cn for j > n have to be adjusted in such a way that Y d in (B1) is real and non-negative. When Y is perturbed by an O(ε) amount with ε being a small parameter, one would naively expect also O(ε) changes in the unitary matrices U, U C in (B1). The following statement confirms this expectation under certain assumptions.
Further, let X be an arbitrary complex matrix and let us defineX = U T C XU . Theñ where R d is a real diagonal matrix with R ii d = ReX ii and U = U (1 + εZ) + O(ε 2 ),Ũ C = U C (1 + εZ C ) + O(ε 2 ) (B7) for some antihermitian matrices Z and Z C .
Proof. Let us assume that the generalized eigenvalues of Y + εX and also the corresponding rotation matrices U, U C are changed by the O(ε) values, i.e., the diagonalization of Y + εX follows (B6) withŨ ,Ũ C of the form (B7). ForŨ ,Ũ C to be unitary (up to O(ε 2 ) terms), Z † = −Z, Z † C = −Z C has to be satisfied, hence, Z, Z C are indeed antihermitian. In order to prove the lemma, the matrices R d and Z, Z C will be explicitly constructed.
According to (B2), the matrixŨ is defined by If the shape ofŨ (B7) is plugged in, then the equality of the O(ε) terms yields Multiplying this relation by U † from the left and by U from the right one obtains whereX = U T C XU was defined. Taking the diagonal elements of this matrix relation, one obtains (after dividing by Y jj d )X jj +X jj * = 2R jj d which shows that indeed R jj d = ReX jj as stated in the lemma. On the other hand, the off-diagonal elements of the matrix relation (B8) allow to compute the matrix Z: (B9) Finally, any purely imaginary number can be chosen as the diagonal entries of Z which corresponds to the ambiguity in the definition of the rotation matrices mentioned in point i) above.
The matrix Z C can be constructed analogously when (B4) is taken into account and Finally, one can plug in the shape ofŨ andŨ C into the formula (B6) and the O(ε) part of this relation then yields If the formulas for Z and Z C are plugged in, one can easily check that indeed the off-diagonal elements of the lefthand side are equal to zero. Moreover, if the imaginary part of the diagonal elements is evaluated, one obtains which fixes the (purely imaginary) diagonal entries of Z C .
It is now easy to understand why the above described construction breaks down when Y 11 No information about Z ij for i, j ≤ n can be obtained from (B8) (and analogously, no information about Z ij C for i, j ≤ n is available). Moreover, (B11) simplifies toX ij = R ij d for i, j ≤ n, hence, the upper left n × n block ofX has to be diagonal. This can be ensured thanks to the ambiguity (B5) in the definition of the U and U C matrices in (B1), with U n and U Cn being chosen in such a way that U T CnX U n ij = U T Cn U T C XU U n ij = 0 i = j, i, j ≤ n.
(B12) The shape of the perturbed rotation matrices (B7) given in Lemma 1 can be, hence, still used, however, the ambiguity (B5) is lifted.
Furthermore, let us consider the setting with O(ε) generalized eigenvalues, more precisely let where Y d and Λ n are diagonal matrices with Y jj d = 0 for j ≤ n and Λ jj n = 0 for j > n. In order to illustrate the effect of an O(ε) perturbation in this case, let us define Y 0 = Y − εU * C Λ n U † which obviously has first n generalized eigenvalues equal to zero and the considerations of the previous paragraph can be applied. Y can be then viewed as a perturbation of Y 0 by an O(ε) term which lifts the ambiguity in the definition of U and U C as described in (B12). Similarly, Y + εX can be understood as a different perturbation of Y 0 and clearly, in general, U and U C are fixed in a different way. It is then easy to check the following statement. Lemma 2 Let Y be a complex matrix diagonalized by a biunitary transformation as in (B13), let X be an arbitrary complex matrix and ε a small parameter. Theñ whereR d is a diagonal matrix and the rotation matrices may be written in the form U = U U n (1+εZ)+O(ε 2 )Ũ C = U C U Cn (1+εZ C )+O(ε 2 ) (B15) for some antihermitian matrices Z, Z C and unitary matrices U n , U Cn where U ij n = U ij Cn = 0 for i, j > n, i = j.
In our work we are interested in perturbations preserving the generalized eigenvalues of the original matrix. The above results on the shape of the diagonalization matrices can be used also in this case due to the following simple observation.

Lemma 3
For any complex matrix X there exists a complex matrix X such that the generalized eigenvalues of Y + εX differ from the generalized eigenvalues of Y by at most O(ε 2 ) terms and, at the same time, Y + εX and Y + εX are diagonalized by the same biunitary transformation.
Proof. Looking at the relation (B14) it is enough to define X = X −Ũ *