Semi-leptonic three-body proton decay modes from light-cone sum rules

Using light-cone sum rule techniques, we estimate the form factors which parametrise the hadronic matrix elements that are relevant for semi-leptonic three-body proton decays. The obtained form factors allow us to determine the differential rate for the decay of a proton ($p$) into a positron ($e^+$), a neutral pion ($\pi^0$) and a graviton ($G$), which is the leading proton decay channel in the effective theory of gravitons and Standard Model particles (GRSMEFT). The sensitivity of existing and next-generation neutrino experiments in detecting the $p \to e^+ \pi^0 G$ signature is studied and the phenomenological implications of our computations for constraints on the effective mass scale that suppresses the relevant baryon-number violating GRSMEFT operator are discussed.


Introduction
Proton decay constitutes one of the most sensitive probes of high-scale physics beyond the Standard Model (BSM). Most of the existing proton decay searches have focused on two-body decay channels, since transitions such as p → e + π 0 that involve a positron (e + ) and a neutral pion (π 0 ) are generically predicted by theories of grand unification (GUTs). In a broader context, studying baryon-number violation can be interesting in the light of supersymmetric theories, baryogenesis and theories of quantum gravity, where the global symmetries of the Standard Model (SM) are expected to be broken at some level -for a review on baryon-number violation in various BSM models see [1]. The null results provided by these experiments (cf. [2][3][4] for comprehensive summaries of the available experimental results) together with the observation that many higher-dimensional operators violating baryon number (B) by one unit induce multi-body proton decay modes, make proton decay processes with more complicated final states interesting search targets for existing and next-generation neutrino experiments like Super-Kamiokande (SK) [5], Hyper-Kamiokande (HK) [6], JUNO [7] and DUNE [8].
While for all two-body proton decays into anti-leptons and pseudoscalar mesons, lattice QCD (LQCD) techniques nowadays allow direct computations of the relevant hadronic matrix elements within uncertainties of (10 − 15)% [9][10][11][12][13][14][15][16], LQCD calculations of threebody decay modes do not exist at present although the formalism and methodologies are in principle known [17]. Model estimates of three-body final-state proton decay rates are therefore available only for selected modes [18] or rely on naive dimensional analysis and phase-space arguments [3,4]. In our recent article [19] we have shown that by employing light-cone sum rules (LCSRs) [20][21][22][23][24][25] it is possible to reproduce the LQCD results for the hadronic matrix elements relevant for GUT-like proton decay. The goal of this work is to apply the LCSRs formalism developed in the latter publication to the case of semileptonic three-body proton decay processes. In particular, we will describe in detail the calculation of all form factors needed to compute the differential decay rate for the process p → e + π 0 G with G denoting a graviton. This decay mode is the leading proton decay channel in the effective theory that describes the interactions of gravitons and SM particles aka GRSMEFT [26,27]. In a companion paper [28] we will analyse a broad range of possible laboratory probes of the GRSMEFT, showing that proton decay measurements set the nominally strongest bound on the effective mass scale that suppresses the GRSMEFT interactions. While in this work the focus lies on obtaining predictions for p → e + π 0 G, the provided analytic expressions and numerical results can be used to obtain the differential decay rates of other proton decay modes with a single pion in the final state.
This article is structured as follows. In Section 2 we introduce the relevant baryonnumber violating GRSMEFT interactions and provide a suitable representation for the decay amplitude of p → e + π 0 G in terms of a leptonic and a hadronic part. In Section 3 we outline the calculation of the relevant hadronic matrix elements using LCSR techniques, while the structure of the resulting LCSRs is discussed in Section 4. We turn to the numerical evaluation of the LCSRs in Section 5, providing predictions and uncertainty estimates for the proton-to-pion form factors in the physical region. The differential decay rate of p → e + π 0 G is computed in Section 6. This section also contains a discussion of the sensitivity of existing and next-generation neutrino experiments to the studied proton decay signature. Our conclusions and an outlook are presented in Section 7. Supplementary material is relegated to a number of appendices.

Preliminaries
Considering the case of one generation of fermions, baryon-number violation ( / B) is induced by only a single dimension-eight operator in the GRSMEFT [26,27]. We write this operator in the following way where u and d denote the up-and down-quark field, e is the electron field, P R projects on right-handed fields, abc is the fully antisymmetric Levi-Civita tensor with a, b, c denoting colour indices, C is the charge conjugation matrix, T denotes the transpose with respect to Dirac indices and σ µν ≡ i (γ µ γ ν − γ ν γ µ ) /2 with γ µ the usual Dirac matrices. Furthermore, C µνρσ represents the Weyl tensor which is the traceless part of the Riemann tensor R µνρσ . It takes the form with g µν the metric tensor, R µν the Ricci tensor, R the Ricci scalar and the brackets denote index antisymmetrisation, i.e.
Notice that the Wilson coefficient c / B entering (2.1) carries mass dimension −4. The amplitude for the decay p(p p ) → e + (p e )π 0 (p π )G(p) can be written as where κ ≡ 2/M P with the reduced Planck mass M P ≡ 1/ √ 8πG N 2.435 · 10 18 GeV, and G N 6.709 · 10 −39 GeV −2 [29] is the gravitational or Newton constant. In addition, u p (p p ) denotes the spinor of the proton with four-momentum p p ,v c e (p e ) is the charge conjugate anti-spinor of the electron with momentum p e and ε * µρ (p, λ) denotes the conjugate of the polarisation tensor of the graviton with four-momentum p and polarisation λ. The variable q ≡ p p − p π = p + p e denotes the four-momentum transfer from the proton to the neutral pion, and enters A p → e + π 0 G through the hadronic tensor H µν (p p , q). Notice that in order to obtain (2.3) the gauge of the graviton is chosen such that the polarisation tensor is transverse and traceless, i.e. the following terms can be omitted in the decomposition of the decay amplitude (2.3)

Hadronic form factors
In what follows we discuss the necessary steps to calculate the hadronic tensor H µν (p p , q) with the help of LCSRs in QCD. We employ the notation and the conventions introduced in our earlier article [19]. The starting point for the sum rules is the correlation function where T denotes time ordering and the current η is a combination of three quark fields that interpolates the proton 0|η(0)|p(p p ) = m p λ p u p (p p ) .
Here m p 938 MeV is the proton mass and λ p denotes the coupling strength of the current η to the physical proton state. The strongly-interacting part of the dimension-eight operator (2.1) is encoded by By following the standard procedure the hadronic representation of the sum rules can be cast into the form with > 0 and infinitesimal and the ellipsis denotes contributions from heavier states, i.e. excited states and the continuum. The hadronic tensor H µν (p p , q) that characterises the p → π 0 transition can be parameterised by four independent form factors w n with n = 1, 2, 3, 4 in the following way Here the proton spinor is understood to be on-shell, µνρσ is the fully antisymmetric Levi-Civita tensor with 0123 = 1 and we have introduced the abbreviations σ µp ≡ σ µν p ν , σ pq ≡ σ µν p µ q ν and µνpq ≡ µνρσ p ρ q σ . We add that after making use of the Dirac equation and algebraic identities the result (3.5) matches the decomposition provided in [30,31]. The hadronic representation of the correlation function (3.4) can be written as where µνρp ≡ µνρσ p σ . The eight Dirac structures in (3.6) can be used to derive LCSRs for the four form factors w n or combinations of them. The corresponding scalar functions Π had α with α = S, A 1 , A 2 , A 3 , A 4 , T 1 , T 2 , T 3 depend only on the square p 2 p of the proton fourmomentum and on the square Q 2 ≡ −q 2 of the four-momentum transfer between the proton and the neutral pion. They can be expressed as dispersive integrals as follows are spectral densities. In this way the ground-state contribution can be separated from the contribution due to heavier states denoted by ρ cont α (s, Q 2 ): On-shell, i.e. for s = m 2 p , the functions W α (Q 2 ) ≡ W α (m 2 p , Q 2 ) take the following form with m π 135 MeV denoting the mass of the neutral pion.

Structure of LCSRs
The derivation of the QCD results for the LCSRs proceeds in full analogy to Section 3 of our earlier work [19] to which we refer the interested reader for all technical details. The analytic expressions for the QCD correlation functions relevant in the context of this article can be found in Appendix B. Rather than repeating the necessary steps to derive them, let us discuss the structure of the Π QCD α functions. A striking feature of the results for the QCD correlation functions is that at the lowest order in the twist expansion. The first non-zero correction to the QCD correlation function Π QCD schematically takes the form Here qq is the quark condensate, q ≡ (u d) T , τ 3 ≡ σ 3 /2 with σ 3 = diag (1, −1) the third Pauli matrix, g s is the QCD coupling constant and we have defined G αβ ≡ αβµν G A µν T A /2 with G A µν the QCD field strength tensor and T A the SU (3) colour generators. The contribution (4.2) corresponds to a three-particle pion distribution amplitude (DA) of twist 4 with the fields evaluated at 0, ux and x. See for instance [32,33] for details. Being of higher twist the correction (4.2) is expected to be small compared to the values predicted for the form factor w 2 by the LCSRs for Π A 1 and Π A 2 -cf. (3.10). This implies that there has to be a cancellation among the contributions of the ground state and that of heavier states in the hadronic sum leading to Π had T 3 0. Similar cancellations are also observed in certain QCD sum rules for the nucleon mass [34,35]. In this case, for a specific choice of the nucleon interpolating current, one of the sum rules starts at higher order in the operator production expansion (OPE), which numerically yields a very small value for the QCD side of the sum rule. In this example, on the hadronic side excitations of the nucleon with even and odd parity contribute with opposite signs leading to a cancellation. In the case of (4.1) the contributions of excited states are not sign-definite, but in principle cancellations may occur if the contributions from heavier states are sizeable. Sum rules that have this feature cannot be used to extract the form factors related to the ground state because the corrections of excited states are just as important as the formally leading ground-state contributions. The LCSR for the correlation function Π T 3 is thus disregarded in our work.
Convergence criteria are now applied to the remaining LCSRs in order to determine the Borel window for each Π α . Notice that compared to the correlator studied for GUT-like proton decay in [19] the hadronic representation of the correlation function (3.6) comprises a larger number of independent Lorentz structures. This feature leads to simpler analytic LCSR expressions for the correlation functions Π α , but it also renders the numerical impact of the dimension-five condensate qg s G · σq with G · σ ≡ G µν σ µν larger than in the GUT case. As a result the LCSRs analysed below will have larger uncertainties than those that have been studied in [19]. In the following, we will use the LCSRs for Π S , Π A 1 , Π A 2 and Π T 1 to extract the form factors w n because they are the most well behaved with regard to the power suppression of higher-dimensional condensates and the dominance of the ground-state contributions. We add that the LCSR for Π A 4 also fulfils the convergence criteria but one would need to combine it with the result for Π T 1 to extract the form factor w 3 and it turns out that the Borel windows of these two LCRSs do not overlap.
In the case of the LCSRs for Π S and Π A 2 we find the window 1.1 GeV M 1.5 GeV with M the Borel mass, while for Π T 1 we obtain 0.7 GeV M 1.1 GeV. In all three cases the Borel analysis has been restricted to 0.6 GeV 2 Q 2 2.5 GeV 2 . The lower limits are obtained by demanding that the mixed condensate qg s G · σq does not account for more than 50% of the total QCD result and as an absolute minimum of the Borel mass we choose 0.7 GeV. The upper limits arise from the requirement that heavier states constitute less than 50% of the total dispersion integrals (3.7) and that the Borel mass should not considerably exceed the continuum threshold s 0 = 1.44 GeV 2 . The latter is chosen as the square of the mass of the Roper resonance [29] which is the lightest excitation in the nucleon spectrum. The Borel transformation ensures that heavier states with a mass m N are exponentially suppressed by a factor of exp −(m 2 N − m 2 p )/M 2 , but the dispersion integral over heavier states, which starts at s 0 , can be modelled as an integral over the QCD result assuming quark-hadron duality [36,37].
In order to determine all four form factors w n one also needs to evaluate Π A 1 , where the power suppression turns out to be less effective than in the other cases. Therefore the use of this LCSR is restricted to Q 2 0.9 GeV 2 because the power suppression becomes more effective for larger values of Q 2 . For Borel masses of 1.1 GeV M 1.5 GeV the relative contribution of the dimension-five condensate qg s G · σq to the LCSR lies between 60% and 100%. The form factor w 2 can therefore only be estimated within systematic uncertainties of the order of 100%. The result of w 2 also enters the prediction for w 3 through the LCSR for Π A 2 but the contribution is suppressed by a kinematical factor of about 0.1 for Q 2 1 GeV 2 cf. (3.10) . As a result the uncertainties plaguing w 2 represent only a subleading part of the uncertainty in w 3 . In fact, it turns out that the differential decay width of p → e + π 0 G receives the dominant contributions from the form factors w 1 and w 3 , meaning that the uncertainty due to w 2 plays only a minor role in the proton decay phenomenology in the GRSMEFT.
If one could resum the expansion of the QCD side to all orders, the dependence on the Borel mass M of the form factors would vanish. Truncating the expansion at some finite order leaves a residual dependence on this parameter, but ideally the results for the form factors w n do not depend too strongly on the exact choice of this scale. Moreover,    observables should not depend on the choice of the effective continuum threshold s 0 which enters the LCSRs because like in [19] we assume quark-hadron duality [36,37]. These two scales are therefore referred to as unphysical parameters in the following.

Numerical analysis
Using the numerical input and the definitions of the pion DAs of Appendix A we obtain the results for the form factors w n shown in Figure 2. The displayed central values of w n correspond to M = 1.3 GeV for Π S , Π A 1 and Π A 2 , while in the case of Π T 1 we use M = 0.9 GeV. All central predictions employ s 0 = (1.44 GeV) 2 . The total theoretical uncertainties receive contributions from variations of M and s 0 as described above but also from variations of the numerical input parameters within their uncertainties (cf. Appendix A). Each parameter is varied independently while the remaining parameters are kept fixed at their central values. The total uncertainty is then obtained by adding individual uncertainties in quadrature. The values of the form factors w 1 (Q 2 ) and w 4 (Q 2 ) are computed with the help of the LCSRs for Q 2 ≥ 0.6 GeV 2 while for Q 2 ≤ 0.6 GeV 2 the values are obtained by a naive extrapolation. Similarly, the form factors w 2 (Q 2 ) and w 3 (Q 2 ) are predicted by the LCSRs for Q 2 ≥ 0.9 GeV 2 and by the extrapolation for Q 2 ≤ 0.9 GeV 2 . A linear and a quadratic function in Q 2 is taken to extrapolate the form factors which are then fitted to the results of the values obtained from the LCSRs in the vicinity of Q 2 = 0.6 GeV 2 for w 1 (Q 2 ) and w 4 (Q 2 ) and Q 2 = 0.9 GeV 2 for w 2 (Q 2 ) and w 3 (Q 2 ). For a given form factor the quadratic fit is chosen to obtain the central value for w n , while the maximum and minimum of all extrapolations determine the uncertainty band. We remark that the same fitting approach has been successfully used in the LCSR calculation [19] to reproduce the results from LQCD in the case of a GUT-like proton decay. The fit formulas for the form factors w n that we obtain in the physical region, i.e. in the four-momentum range −0.65 GeV 2 −(m p − m π ) 2 ≤ Q 2 ≤ 0, take the following form  Figure 2. Form factors w n (Q 2 ) as a function of Q 2 . The coloured curves and bands correspond to the central values and uncertainties of the LCSRs. In the case of w 1 (Q 2 ) and w 4 (Q 2 ) w 2 (Q 2 ) and w 3 (Q 2 ) the predictions for 0.6 GeV 2 ≤ Q 2 ≤ 1.5 GeV 2 (0.9 GeV 2 ≤ Q 2 ≤ 1.5 GeV 2 ) are obtained by a direct calculation (solid lines), while the predictions for Q 2 ≤ 0.6 GeV 2 (Q 2 ≤ 0.9 GeV 2 ) are obtained by an extrapolation (dashed lines). Consult the main text for further explanations. 3) Here the upper (lower) line in each formula corresponds to the upper (lower) border of the corresponding envelope shown in Figure 2, while the middle line represents the central value of our LCSR form factor prediction. Let us also discuss alternative approaches for extrapolating the LCSR results for the form factors to the physical regime, i.e. to negative Q 2 . A physically well-motivated approach is to model the lowest-lying resonance(s) in the Q 2 -spectrum explicitly by poles in the complex t ≡ −Q 2 plane and capture modifications to the spectrum due to other contributions by an additional series expansion in a new variable z(t) -see for example [38][39][40][41][42] for details on this so-called z-expansion. This parametrisation is often used to take into account resonances that lie below the threshold of the two-particle branch cut which in our case is located at t = t th ≡ (m p + m π ) 2 . However, the lowest-lying resonance is heavier than t th in the case at hand so it is not clear whether the form factors are dominated by a pole contribution close to the upper limit of the allowed three-body kinematics, i.e. t = (m p − m π ) 2 , or by the continuum two-particle contribution. Employing the z-expansion up to the second order including a single pole at the mass m ∆ + = 1.232 GeV [29] of the ∆ + resonance leads to a steeper increase in magnitude of the form factors w n at small Q 2 and in some cases to significantly narrower uncertainty bands. This approach would therefore lead to larger predictions for the form factors w n . In another approach used in [33,43], all contributions but the lowest-lying resonance are modelled by an effective pole at higher mass which is more flexible than a single-pole fit. For the form factors w n this procedure however generates unphysical singularities in the regime of physical Q 2 , rendering it unsuitable. From the above, we conclude that our power-expansion approach (5.1) to (5.4) yields a more conservative estimate of the form factors w n at small Q 2 than the z-expansion including a single pole and is more reliable than the effective two-pole fit. It is therefore preferred with respect to extrapolations featuring explicit poles in the Q 2 -spectrum for the sensitivity studies of semi-leptonic proton decay searches that are discussed in the following section.
We add that the form factors w n are related to the off-shell form factors of the decomposition of the more general matrix element π 0 | abc d α a u β b u γ c |p [30,31], where α, β and γ are Dirac indices. Certain combinations of the form factors w n therefore yield the form factors W k RR with k = 0, 1 that are relevant for GUT-like proton decay. In Appendix C we show that using the results (5.1) to (5.4) allows one to reproduce the physical values of the form factors W k RR calculated in our previous work [19] within uncertainties. This gives us confidence that the naive extrapolation used to obtain the above expressions for w n sufficiently approximates the true scaling in the relevant four-momentum regime.

Proton decay phenomenology
With the help of the expressions (5.1) to (5.4) the p → e + π 0 G decay amplitude (2.3) can now be calculated. One first notices that after making use of the on-shell conditions for the graviton cf. (2.4) the contribution of w 4 vanishes. This feature can be understood by means of the soft pion theorem [44,45]. In fact, in the soft pion limit and recalling that q = p p − p π one finds that all terms but the contribution of w 4 vanish in the hadronic tensor: lim pπ→0 H µν (p p , q) = H µν (p p , p p ) = iw 4 P R σ µν u p (p p ) . (6.1) However, in the limit p π → 0 the pion can be removed from the decay amplitude giving rise to the following relation lim pπ→0 e + Gπ 0 (p π )|L denotes the axial current and the pion field is related to axial current by π 0 ( The second term in (6.2) vanishes unless there are additional poles in the soft pion limit. Such poles occur if the pion is attached to one of the external lines [45] in the p → e + G amplitude, which is formally described by inserting a complete set of intermediate states between the operators in the time-ordered product. The pion however can only couple to the external proton line, so pole contributions arise only when the pion is emitted from the incoming proton. This type of correction thus leads again to the matrix element of the p → e + G transition, which however satisfies e + G|L Since the form factor w 4 itself is non-vanishing it then follows that the associated Lorentz structure (6.1) does not contribute to the proton decay channel p → e + π 0 G at all. By squaring the amplitude, summing over spins and polarisations and calculating the phase space integrals the differential decay width can be computed. Note that the transversality of the graviton see (2.4) has been used to drop unphysical contributions which violates the gauge symmetry of gravity in the weak field limit. The gauge symmetry ensures that negative-norm states cancel out in the sum over polarisations. Therefore the sum has to be constrained to physical polarisations only by employing [46] λ ε * αβ (p, λ) ε γδ (p, λ) = Neglecting the mass of the positron but keeping the mass of the neutral pion, the p → e + π 0 G rate corresponding to the fiducial region of the three-particle phase space defined by an upper cut on the graviton energy can be written as Here we have defined y ≡ 2E G /m p and z ≡ 2E e /m p with E G (E e ) the graviton (positron) energy in the rest frame of the proton, x π ≡ m 2 π /m 2 p , the boundaries for the integral over z are given by and When expressed through the integration variables of (6.6) the scale Q 2 that enters the form factors w n finally takes the following form The GRSMEFT proton decay mode p → e + π 0 G experimentally leads to events that contain a positron, two photons arising from the decay of the neutral pion and missing energy (E miss ) because the graviton escapes the detector undetected. Such a signature has to our knowledge not been searched for directly in experiments that study the possible decays of the proton. As we will show, existing searches that are however sensitive to p → e + π 0 G are the inclusive search p → e + X with X an arbitrary final state and the exclusive search for p → e + π 0 . The total inclusive rate p → e + X can be obtained by employing y fid = 1 − x π in (6.6). Numerically, we find that Λ p = (99 ± 13) MeV . (6.10) Here we have introduced the hadronic parameter Λ p and the normalisation factor 1/(256π 3 ) takes into account the phase-space suppression for a three-body decay. The uncertainty on Λ p is obtained by calculating the minimal and maximal rate that can be achieved by considering all possible combinations of the form factor parameterisations (5.1) to (5.4). Notice that since Λ p appears in (6.10) to the fourth power the LCSR prediction for Γ p → e + π 0 G has an uncertainty of order 50%. The theory uncertainties of Γ p → e + π 0 G are therefore significantly larger than those that plague the GUT predictions for Γ p → e + π 0 obtained in both LQCD [9][10][11][12][13][14][15][16] and LCSRs [19]. Searches for the two-body decay mode p → e + π 0 can also be used to set a bound on the GRSMEFT interaction (2.1), because the cuts that experiments such as SK impose do not fully eliminate the contributions that arise from p → e + π 0 G. The relevant requirements in these experiments are selections on the invariant mass m eπ and the magnitude of the threemomentum p eπ of the e + π 0 system. In the rest frame of the proton these quantities can be expressed in terms of the graviton energy as m eπ = m p (m p − 2E G ) and p eπ = E G . In the latest SK search [47] the definition of the signal region involves the requirements m eπ > 800 MeV and p eπ < 250 MeV, meaning that all p → e + π 0 G events that satisfy the m eπ selection also pass the p eπ cut. To quantify by how much a lower cut m eπ > m cut eπ reduces the observed p → e + π 0 G rate we introduce the acceptance A m cut eπ = Γ p → e + π 0 G fid Γ (p → e + π 0 G) , (6.11) where Γ p → e + π 0 G fid is the fiducial width (6.6) evaluated setting y fid = 1 − m cut eπ /m p 2 and Γ p → e + π 0 G is the total inclusive width given in (6.10). In Figure 3 we show our predictions for A m cut eπ in the range of m cut eπ that is relevant for searches for the two-body proton decay mode p → e + π 0 at existing and next-generation water Cherenkov detectors.
We are now in a position to derive bounds on the Wilson coefficient c / B that multiplies the dimension-eight operator in (2.1). We begin with the inclusive p → e + X decay. The currently best proton lifetime limit from p → e + X is unfortunately more than 40 years old. It reads [48] τ p p → e + X > 0.6 · 10 30 yr , (6.12) and together with (6.10) leads to at the 90% CL. It has been pointed out in [3] that with existing data from water Cherenkov detectors it should be possible to set a significantly better limit on p → e + X compared to the proton lifetime limit reported in (6.12). An estimate of such an improved limit can be obtained from the limit of 1.7·10 32 yr at 90% CL on the p → e + +E miss channel [49] since the latter decay bears close resemblance to the inclusive p → e + X mode. In fact, the authors of [3] estimated that with the available SK data it should be possible to improve (6.12) by around two orders of magnitude. Note that a factor of 100 improvement on τ p (p → e + X) would push the bound (6.13) up to (186 GeV) −4 . In order to determine the SK sensitivity to the p → e + π 0 G signature that derives from the measurement [49], we need to compute the acceptance (6.11) for m cut eπ = 800 MeV. Using the central value of the hadronic parameter Λ p as given in (6.10) we find A (800 MeV) = 2.7 · 10 −3 (1 ± 0.40) . (6.14) The smallness of the acceptance is compensated by the fact that the 90% CL lower limit on the lifetime of the proton in p → e + π 0 is by more than five orders of magnitude better than (6.12) since one has [47] τ p p → e + π 0 > 2.4 · 10 34 yr . (6.15) Combining (6.10) for Λ p = 99 MeV with (6.14) and (6.15) one obtains at the 90% CL. Notice that this bound is very close to the limit that has been quoted above by assuming a factor of 100 improvement of τ p (p → e + X) compared to (6.12) based on the estimate presented in [3]. It is also straightforward to estimate the sensitivity of HK to the Wilson coefficient c / B . Running HK for eight years it should be possible to set the following 90% CL bound τ p p → e + π 0 > 1.0 · 10 35 yr . (6.17) This limit has been obtained in [6] by considering the same signal region for p → e + π 0 as the latest SK search [47]. Consequently, we can again use (6.10), (6.14) and (6.17) to arrive at This bound on the Wilson coefficient of the dimension-eight baryon-number violating operator (2.1) is probably the ultimate limit that can be set with next-generation neutrino detectors, because both JUNO and DUNE are not expected to reach the HK sensitivity to the p → e + π 0 mode cf. [7,8] .

Conclusions and outlook
In our work, we have employed LCSR techniques to calculate the form factors which parametrise the hadronic matrix elements of semi-leptonic three-body proton decays with a pion in the final state. While the presented formalism and the obtained results are general, we have specifically focused on the computation of the differential decay rate for the process p → e + π 0 G. This channel is the dominant proton decay mode in the GRSMEFT, since the two-body transition p → e + G is forbidden by angular momentum conservation. Like in our earlier article [19] our LCSR study includes the leading contributions in the lightcone expansion, namely the twist-2 and twist-3 pion DAs -the explicit expressions can be found in Appendix B -and we have performed a detailed study of the dependence of the obtained LCSRs on both the unphysical (i.e. the Borel mass and the continuum threshold) and the physical (i.e. the condensates and the pion DAs) parameters. In this way we are able to provide results and estimate uncertainties for the form factors in the kinematical regime where the momentum transfer q from the proton to the pion is space-like. We have then extrapolated our LCSR results to the physical regime 0 ≤ q 2 ≤ (m p − m π ) 2 by means of both a linear and quadratic fit, including the spread of predictions in our uncertainty estimates. The resulting uncertainties turn out to be significantly larger than those that plague the hadronic matrix elements that are relevant in the GUT case [9][10][11][12][13][14][15][16]19].
Employing the LCSR results for the form factors we have then studied the sensitivity of existing and next-generation water Cherenkov detectors in looking for the p → e + π 0 G signature. To this purpose we have calculated the rate for p → e + π 0 G differentially in the energies of the final state particles. We have then derived bounds on the amount of dimension-eight baryon-number violation in the GRSMEFT considering both the inclusive search for p → e + X [3,48] and the exclusive search for p → e + π 0 . It turns out that the best constraint arises at the moment from the latest SK search for the two-body decay mode p → e + π 0 [47]. In fact, this search is able to set a 90% CL lower limit of 185 GeV on the effective mass scale that suppresses the relevant baryon-number violating GRSMEFT interactions. HK measurements are expected to be able to push this limit up to 222 GeV. In a companion paper [28] we will analyse a broad range of possible laboratory probes of dimension-six and dimension-eight GRSMEFT operators, showing that the proton decay measurements discussed here set the nominally strongest bound on the effective mass scale that suppresses the GRSMEFT interactions. renormalisation group running and the one-loop threshold corrections as implemented in RunDec [55,56], we obtain at 1 GeV the value m u + m d = (8.60 ± 0.11) MeV. By means of the Gell-Mann-Oakes-Renner (GMOR) relation m 2 π −2 (m u + m d ) qq /f 2 π [57] this value together with the pion decay constant f π = (130.2 ± 0.8) MeV [54] leads to if the leading-order chiral corrections of [58] are included and uncertainties are added in quadrature. In the case of the mixed condensate we employ [59] qg s G · σq = m 2 0 qq , The non-perturbative parameters µ π and ρ π can be fixed via the GMOR relation: This leads to µ π = (1.98 ± 0.05) GeV , ρ π = 0.068 ± 0.002 .
For the pion DAs we use the following expressions which have been derived in [32] (and [60] in the chiral limit) with the help of a conformal expansion. One has whereū ≡ 1 − u and the expansion in terms of the Gegenbauer polynomials C (m) n (ζ) with ζ ≡ 2u − 1 is truncated after n = 4. The hadronic parameters that enter the above definitions depend on the renormalisation scale µ which we set equal to 1 GeV in our numerical analysis.
We adopt the numerical values of the two Gegenbauer moments given in [61], where the moments are obtained by fitting sum rules for the electromagnetic pion form factor to the experimental data of [62]. For the numerical values of the other parameters we rely on the sum rules estimates presented in [63]. Adding all uncertainties in quadrature this leads to together with the definition f 3π (µ) ≡ f π µ π η 3 (µ) and (A.4) are used. In order to obtain the coupling strength λ p we use the QCD sum rule [64] The parameters M ands 0 denote the Borel mass and the continuum threshold of (A.12). The corresponding windows are 0.7 GeV ≤ M ≤ 1 GeV and (1.4 GeV) 2 ≤s 0 ≤ (1.5 GeV) 2 . We furthermore employ [59] α s π G 2 = (0.009 ± 0.009) GeV 4 . (A.14) We add that using the results for the condensates and the variation of unphysical scales specified above, the sum rule (A.12) leads to a good agreement with the LQCD results for λ p (see for instance [9,[65][66][67]).

B Analytic results for LCSRs
Below we provide the analytic expressions for the QCD correlation functions that are relevant for the calculation of the process p → e + π 0 G in the GRSMEFT. The hat on the functionsΠ QCD α indicates that we have subtracted the contributions of heavy states before taking the Borel transform of the QCD results. We obtain where as before q = p p − p π . Other transitions due to operators of the type (C.1) with different chiralities for the quark fields exist and contribute to GUT-like proton decay. But only the matrix element (C.2), where all quarks are right-handed, is related to the decay induced by the GRSMEFT operator (2.1).
The most general decomposition of H αβγ (p p , q) in terms of a set of form factors for an off-shell proton is provided in [30,31] -see in particular (4.64) of the current arXiv version of [31]. Hence, both sets of on-shell form factors w n and W k RR can be related to these off-shell form factors upon using the equations of motion for the proton. In this way it is possible to relate the on-shell form factors for both scenarios of proton decay among each other. We find Numerical results for W k RR (Q 2 ) were computed in our recent work [19], where our findings were also shown to be in agreement with the results of the state-of-the-art LQCD calculation [15] at Q 2 = 0. The form factors w n (Q 2 ) are derived in the same way in this work. One can therefore employ the relations (C.6) and (C.7) to directly assess the validity of first the numerical results presented in Figure 2 for large virtualities (Q 2 Λ 2 QCD with Λ QCD 300 MeV the QCD scale) and second the extrapolations (5.1) to (5.4) for physical momenta (Q 2 0). We find that the relations (C.6) and (C.7) hold numerically within uncertainties in the relevant regime −(m p −m π ) 2 ≤ Q 2 ≤ (2 GeV) 2 , even though the uncertainties of the combinations on the right-hand sides of (C.6) and (C.7) are rather large and the results tend to undershoot the more accurate results for W 0 RR (Q 2 ) and W 1 RR (Q 2 ) obtained by a direct calculation of the corresponding left-hand sides. This feature is illustrated by the two panels in Figure 4. Although the agreement is not perfect, we believe that the shown results validate the LCSR approach employed in this work as well as the extrapolation procedure used to obtain (5.1) to (5.4).