Light-cone sum rules for proton decay

We estimate the form factors that parametrise the hadronic matrix elements of proton-to-pion transitions with the help of light-cone sum rules. These form factors are relevant for semi-leptonic proton decay channels induced by baryon-number violating dimension-six operators, as typically studied in the context of grand unified theories. We calculate the form factors in a kinematical regime where the momentum transfer from the proton to the pion is space-like and extrapolate our final results to the regime that is relevant for proton decay. In this way, we obtain estimates for the form factors that show agreement with the state-of-the-art calculations in lattice QCD, if systematic uncertainties are taken into account. Our work is a first step towards calculating more involved proton decay channels where lattice QCD results are not available at present.


Introduction
Proton decay would be a clear indication of physics beyond the Standard Model (SM) if it was measured. While the global symmetries of the SM forbid proton decay, these symmetries arise accidentally. In fact, considering baryon-number violation at the perturbative level is well motivated by theories of grand unification (GUTs) [1,2], supersymmetric theories [3][4][5], models of baryogenesis [6,7] and more generally in theories of quantum gravity, where the global symmetries of the SM are expected to be broken at some level [8,9]. In the case of GUTs, baryon number is typically violated by tree-level interactions, and proton decay is mediated by the massive gauge fields of the spontaneously broken unified gauge group [10][11][12][13]. On the experimental side, searches for simple proton decay channels, such as decays into a pseudoscalar meson and an anti-lepton, provide very strong constraints on the proton lifetime of τ p 10 34 years [14]. These bounds allow one to probe theories that predict proton decay up to extremely high energy scales, thereby putting severe constraints on the scale of unification.
Obtaining such constraints requires a theoretical prediction for the widths of semileptonic proton decay channels, which in turn relies on the knowledge of the hadronic matrix element of the underlying proton-to-meson transition. Water-Cherenkov experiments such as Super-Kamiokande [15,16] attempt to measure the decay products of protons that are approximately at rest, and thus the relevant energy scale for the hadronic transition is given by the proton mass. A perturbative description of the relevant hadronic matrix elements in QCD is not possible at this energy scale because of large radiative corrections due to the exchange of soft gluons. A prediction for the hadronic matrix elements by other means is therefore required to probe baryon-number violating new physics with the help of experimental data from proton decay searches.
Early attempts to compute the hadronic matrix elements of semi-leptonic proton decays date back to the '80s and employed non-relativistic quark models, often based on the approximate SU (6) flavour-spin symmetry of the partons [17][18][19][20][21], bag models which allow for relativistic partons [22][23][24][25][26][27], or QCD sum rules [28]. Also an effective chiral theory was proposed in the articles [29][30][31][32], which can be used to derive relations among the various two-body decay widths but still contains a priori unknown low-energy constants. As a result the latter approach cannot predict the absolute value of the proton decay width without further input. The methods mentioned above have also been applied to estimate these low-energy constants, in which case the final predictions for the hadronic matrix elements suffer from additional systematic uncertainties due the approximate nature of the effective chiral theory. Moreover, the results of these model calculations differ by up to an order of magnitude from each other (see Table VI in [33] for a summary and comparison). On the other hand, lattice QCD (LQCD) groups have by now achieved to directly compute the needed hadronic matrix elements within uncertainties of (10 − 15)% [33][34][35][36][37][38][39][40]. These results cover all two-body decays into pseudoscalar mesons and light anti-leptons, which are relevant for GUTs.
Experimental bounds are today available for a broad range of baryon-number violating processes [41][42][43]. In particular, inclusive proton (p) decay searches for processes like p → π 0 + + X that involve a neutral pion (π 0 ) and an anti-lepton ( + ) might be of interest if baryon-number violation does not become manifest in a simple two-body decay. For example, the case where X is a graviton may provide relevant constraints on theories where baryon-number violation occurs in connection with gravity such as in the effective theory of gravitons and SM particles called GRSMEFT [44,45]. LQCD results are not available for processes of this kind, which raises the question: how can one obtain estimates of the proton lifetime in such cases?
In this work, we establish a method that allows to estimate the hadronic matrix elements that enter processes of the type p → π 0 + (+ X). As a proof-of-principle we apply our general approach to the simple two-body case p → π 0 e + with e + a positron, leaving the application to three-body proton decay processes such as transitions involving an additional graviton for future work. In fact, studying the simple decay mode p → π 0 e + allows us to make a thorough comparison with the latest LQCD results [39]. In this way we are not only able to validate our method but can also assess the systematic uncertainties that plague our estimates. Our method employs the techniques of light-cone sum rules (LCSRs) in QCD. In particular, we perform an operator product expansion (OPE) on the light-cone, which allows us to factorise the hard scattering from the soft interactions. While the hard-scattering kernel can be computed perturbatively in QCD, the soft contributions are parametrised in terms of condensates and distribution amplitudes (DAs) of the final-state pion, which enter the LCSRs as input. The light-cone expansion works well if the momentum transfer q from the proton to the pion is large in magnitude and space-like, i.e. q 2 < 0. We therefore cannot directly compute the hadronic matrix elements at the physical point of the two-body decay kinematics, where q 2 is fixed and equal to the square of the positron mass. However, we are able to find values in the space-like regime at q 2 −0.5 GeV 2 , which are close enough to the physical regime to provide an estimate of the hadronic matrix elements at the physical point by means of suitable extrapolations. Albeit our approach does not achieve the same level of accuracy as the state-of-the-art LQCD calculation [39], we believe that the results of this work are promising, because the obtained precision is better than the methods that have been developed in the '80s to estimate proton decay rates. Furthermore, the LCSR approach developed by us in this article should be able to at least provide order-of-magnitude estimates for hadronic matrix elements that enter certain three-body proton decay processes. Such decays could be phenomenologically relevant (see for instance [42]) but only model estimations exist for selected modes [46], making three-body final-state proton decay processes an interesting target of future LQCD studies [47].
Our work is organised is as follows. In Section 2 we put the discussion of hadronic matrix elements for the p → π 0 e + decay on a more systematic footing. In particular, we use all of the dimension-six operators in the so-called SM effective field theory (SMEFT) that are relevant for this decay, so that our analysis is model-independent. These operators are typically generated by baryon-number violating new physics that can be integrated out below a certain (large) energy scale. After that we decompose the hadronic matrix elements into form factors. These form factors enter a correlation function that is computed with the help of LCSR techniques in Section 3, which allows us to derive the LCSRs for the form factors relevant for proton decay in GUTs. In Section 4 we turn to the numerical evaluation of the LCSRs and compute the form factors in the regime of virtual momentum transfer. Eventually, we compare our findings to the results of the latest LQCD computation [39] and discuss in detail the uncertainties that enter our final estimates. We conclude and present an outlook in Section 5. Technical details are relegated to three appendices.

Phenomenological parametrisation
Integrating out heavy new physics that violates baryon number typically generates the following dimension-six SMEFT operators, where C is the charge conjugation matrix, T denotes the transpose of the Dirac index, the symbols P Γ,Γ denote the left-and right-chiral projectors P L and P R such that Γ (Γ ) denotes the chirality of the first (second) fermion bilinear, abc is the fully antisymmetric Levi-Civita tensor and a, b, c are colour indices. In the following analysis we restrict ourselves to the first generation of up-and down-type quarks u and d and the electron e. We consider all possible chirality combinations of the interaction (2.1) in order to provide a model-independent analysis. Note that the Wilson coefficients c ΓΓ , which encode short-distance physics, have a mass dimension of −2.
The transition matrix element of the proton decay p → π 0 e + induced by an insertion of an operator entering (2.1) can be factorised into a hadronic and leptonic part (up to electroweak corrections), Here u p (p p ) denotes the spinor of the proton with momentum p p andv c e (q) is the charge conjugate anti-spinor of the electron with momentum q ≡ p p − p π . The main goal of the following analysis is to calculate the hadronic matrix element H ΓΓ (p p , q) of the p → π 0 transition, where all the quark fields are evaluated at the space-time point x = 0. For an on-shell proton the above matrix element can be decomposed into two form factors as follows, with m p = 938 MeV the proton mass. Notice that the form factors are related due to parity, which is conserved in QCD. Specifically, one has with n = 0, 1. In this work we calculate the combinations ΓΓ = RR, LR explicitly, which covers all chirality combinations due to the above relations. The starting point for evaluating the form factors W n ΓΓ (q 2 ) in LCSRs is the correlation function where T denotes time ordering and the current η (η ≡ η † γ 0 ) is a combination of three quark fields that interpolates the proton, 0|η(0)|p(p p ) = m p λ p u p (p p ) . (2.7) Here λ p denotes the couplings strength of the current η to the physical proton state. The strongly-interacting parts of the dimension-six operators (2.1) are represented by To obtain a parametrisation of the hadronic matrix elements H ΓΓ (p p , q) we insert a complete set of intermediate states that have the same quantum numbers as the proton into (2.6) and isolate the pole contribution of the proton to obtain the hadronic representation of the correlation function: with > 0 and infinitesimal, σ pq ≡ σ µν p µ q ν with σ µν ≡ i/2 (γ µ γ ν − γ ν γ µ ) and the ellipsis denotes contributions from heavier states, i.e. excited states and the continuum. The four independent Dirac structures in (2.9) can be used to derive LCSRs for the form factors W n ΓΓ (q 2 ) or combinations of them. The corresponding scalar functions Π had,α ΓΓ depend only on the square of the proton momentum p 2 p and on the square of the momentum transfer Q 2 ≡ −q 2 . They are conveniently parametrised in terms of dispersion integrals, where α = S, P, Q, T and we have introduced the spectral densities Separating the ground-state contribution from the contribution of heavy states denoted by ρ cont,α ΓΓ (s, Q 2 ), the four spectral densities appearing in (2.9) can be cast into the form and m π = 135 MeV is the pion mass. We stress that the relations (2.13) only hold onshell, i.e. if s = m 2 p . This is however guaranteed by the δ s − m 2 p function appearing in (2.12). Under the assumption of a global quark-hadron duality [48] (see also [49] for a review) the contributions of heavy states can be approximated by (2.14) where ρ QCD,α ΓΓ (s, Q 2 ) are the spectral densities in QCD and we will explain how to compute them in the next section. The approximation (2.14) is expected to work well for a sufficiently large continuum threshold s 0 , which is a free parameter and has to be determined within the LCSR calculation. A more detailed discussion on how to fix s 0 is provided in Section 4, but ideally it is chosen low enough to cover even the lightest excitation, which is the Roper resonance with a mass of 1.44 GeV.

LCSR calculation
The basic idea of the LCSRs is to derive a result for Π ΓΓ (p p , q) in QCD while parametrising unknown soft contributions in terms of quantities that can be determined by other means. It can be shown that for large virtualities Q 2 Λ 2 QCD and P 2 p ≡ −p 2 p Λ 2 QCD with Λ QCD 300 MeV the QCD scale, the integrand of the correlator (2.6) can be approximated by an expansion on the light-cone x 2 ∼ 1/Q 2 0 (see [51] and references therein). Schematically, this light-cone expansion takes the form where the Wilson coefficients C k encode the hard scattering process and the objects O k are composite operators of twist k. The matrix elements of these composite operators correspond to the light-cone DAs of the pion which are non-perturbative objects. Performing a Borel transformation with respect to P 2 p then yields an expansion in inverse powers of the two scales that enter our calculation, i.e. it leads to a power expansion in Λ 2 QCD /M 2 and Λ 2 QCD /Q 2 , where M denotes the Borel mass associated to P 2 p cf. (B.4) . In our article we will provide explicit LCSR expressions that include the leading contributions in this expansion, namely the twist-2 and twist-3 DAs. We will however also comment on the possible impact of twist-4 contributions in all cases where such terms could be phenomenologically relevant (cf. Section 4).
In order to carry out the light-cone expansion, we need to choose an explicit form for the proton current η. The most general choice with the appropriate quantum numbers (at lowest order in derivatives and spin) can be written as a linear combination of the following two currents [52]: The current η 1 excites the ground state as well as heavier states, while η 2 almost exclusively excites heavier states [53]. As a result the coupling strength of η 1 cf. (2.7) to the proton state is larger by a factor of about 100 than that of η 2 . Due to its weak coupling to the proton state, the contribution of the current η 2 is expected to be very small in the case at hand, and we therefore choose for simplicity neglecting a possible admixture of η 2 . Notice that our choice of proton current corresponds to the interpolator usually used in LQCD calculations. The expansion of the time-ordered product that occurs in the twist expansion (3.1) is carried out by partially contracting the quark fields, Here we have employed the following basis of gamma matrices used the notationΓ A ≡ C Γ T A C with C = iγ 2 γ 0 and a summation over the index A is implicit. The pairwise contraction of up (down) quark fields is denoted by S ij u (x) S ij d (x) , i, j, k are colour indices and Tr denotes a trace over Dirac matrices. Hereafter we will work in the isospin limit and will therefore drop the flavour index of the contraction.
With the help of (3.5) it is possible to derive the following completeness relation: The contracted fields need to be expanded for light-like distances (including single-gluon emission), which reads [54] where the ellipsis represents terms that lead to contributions of twist higher than three, g s is the QCD coupling constant, we have employed the short-hand notation G ij µν ≡ G A µν T ij A for the gluon field strength tensor with T ij A the SU (3) generators and definedū ≡ 1 − u. We neglect contributions proportional to the quark masses because they are numerically negligible. In the following we consider only single-gluon interactions which is consistent with truncating the expansion (3.1) after the leading-twist contribution [55]. This leads to the two types of one-loop diagrams that are displayed in the top row of Figure 1.
In addition to these leading-order terms factorised contributions of higher twist and multiplicity turn out to be numerically relevant in the case at hand. Such contributions originate from operators with four quark fields or four quark fields and one gluon, such that only one pair of quarks is contracted in the time-ordered product. Part of the respective amplitudes can be approximated by a factorisation into two-or three-particle DAs (of twist two and three) and vacuum condensates of the remaining quark and gluon fields. Such contributions scale with a smaller power of 1/Q 2 than genuine, non-factorisable terms of higher twist and instead are suppressed by powers of 1/M 2 [56,57]. Effectively, these factorised contributions can be taken into account by replacing one of the contractions S ij (x) in the expression (3.4) by the appropriate local terms and condensates as encoded by [58] ∆S ij (x) = − qq 12 Here the ellipsis denotes higher-dimensional condensates and terms with additional gluons, which are neglected in our work because they are numerically small. The parameter m 0 entering (3.8) is associated with the mixed condensate where G · σ ≡ G µν σ µν . The diagrams resulting from the local expansion (3.8) of the contraction are displayed in the middle and bottom row of Figure 1.
The uncontracted quark bilinears in (3.4) form a pion and still need to be expanded around light-like distances to obtain the light-cone DAs. The pion DAs have been extensively studied in the literature (see [59] for a state-of-the-art discussion), and they have definite twist. The only twist-2 pion DA is given by (cf. for instance [60,61]) where q ≡ (u d) T and τ 3 ≡ σ 3 /2 with σ 3 = diag (1, −1) the third Pauli matrix, while f π denotes the pion decay constant given by f π = (130.2 ± 0.8) MeV [62]. 1 The parameters u andū correspond to the momentum fractions of the two quarks that form the pion. The renormalisation scale µ that appears in the twist-2 pion DA φ (2) (u, µ) is set equal to 1 GeV for most of this work. Higher-order contributions to the matrix element (3.10) arise at the twist-4 level. There are two two-particle twist-3 DAs called σ (u, µ). These are defined by [60,61] We also include the only twist-3 three-particle DA called , which depends on the momentum fractions α d , α u and α g of the down quark, up quark and gluon, respectively, as well as on µ. This object is defined as follows [60,61] The normalisation of the twist-3 DAs contains the sum of the up-and down-quark mass, which fixes the values of the parameters µ π and ρ π as well as the quark condensate via the Gell-Mann-Oakes-Renner (GMOR) relation m 2 π −2 (m u + m d ) qq /f 2 π first derived in the article [63]. One obtains (3.14) The DAs can be obtained by a conformal expansion [60]. Explicit formulas for the DAs appearing in our work are provided in Appendix A. Summing up all the contributions of Figure 1, using (3.10) to (3.13) and performing a Fourier integration allows us one derive an analytic expression for the QCD correlation function The explicit expressions for the analytic results of Π QCD,α ΓΓ p 2 p , Q 2 are somewhat lengthy and therefore provided in Appendix C. For the matching with the hadronic representation (2.9), we also introduce QCD spectral densities like it has been done in (2.10) and (2.11) for the hadronic case. The matching conditions for the LCSRs then read Using quark-hadron duality in the form (2.14), we can subtract the unknown contributions of heavy states from the LCSRs. Effectively, this procedure cuts off the spectral integral computed in QCD at the continuum threshold s 0 . Applying a Borel transformation with respect to P 2 p to both sides of the sum rules suppresses heavy contributions exponentially and generically improves the accuracy of the LCSR approach -the accuracy of our LCSRs will be investigated in Section 4. The Borel transformations also remove all terms that are polynomial in P 2 p , which sets all divergent contributions of the dispersion integrals as well as the ultraviolate (UV) divergences of the diagrams in Figure 1 to zero. The Borel transforms that enter our LCSR analysis as well as other calculational details are given in Appendix B. After Borel transformation the matching conditions (3.16) take the form where the expressions for the form factors W α ΓΓ (s 0 , Q 2 ) can be found in (2.13). By an appropriate combination of the four independent relations (3.17) one can derive two LCSRs for each of the two form factors appearing in (2.4). Hereafter we will refer to these combinations as W 0,P ΓΓ (s 0 , Q 2 ), W 0,S+T ΓΓ (s 0 , Q 2 ), W 1,Q ΓΓ (s 0 , Q 2 ) and W 1,T ΓΓ (s 0 , Q 2 ). Notice that the LCSRs (3.17) depend on two unphysical parameters, namely the continuum threshold s 0 and the Borel mass M . However, the results of a LCSR calculation can only be trusted if the final predictions are to a certain extent independent of the exact choice of s 0 and M . In Section 4 we will provide criteria that allow to assess the convergence properties of (3.17), which we will then us to estimate the uncertainties that plague our LCSR results for the form factors W n,α ΓΓ (Q 2 ). The coupling λ p in (3.17) is in principal known from LQCD calculations (see [58] and references therein), but it can as well be extracted from local QCD sum rules of the two-point correlator where the ellipsis denotes the contributions of heavier states. Using λ p from sum rules has the salient advantage that in this way the uncertainties of the form factors due to the input parameters such as the quark condensate qq or m 2 0 are reduced. 2 We therefore choose to fix λ p with the help of sum-rule techniques rather than to take λ p from LQCD computations. The sum rule derived from the structure / p of the correlator (3.18) is typically disregarded due to uncontrollably large radiative corrections as well as large contributions from heavier states [58]. We thus extract λ p from the sum rule for the mass term. The QCD result for the sum rule can be derived by performing a local OPE around x 0 of the time-ordered product in (3.18). If the momentum flow through the correlator is deeply space-like, i.e. −p 2 Λ 2 QCD , this leads to a convergent expansion of local operators with increasing mass dimension [50]. By plugging the OPE into (3.18) one then obtains an expansion of the correlator in terms of condensates. In this work, contributions including condensates up to dimension seven are included but no perturbative QCD corrections.
Using the proton current (3.3) yields [58] The parameterss 0 and M denote the continuum threshold and the Borel mass of the local sum rule (3.19). These parameters can be related to the corresponding parameters of the LCSRs, because the Borel mass is connected to the momentum flow through the proton current. However, we assume for simplicity thats 0 and M are independent parameters and determine them such that the value of λ p does not depend too strongly on the specific choice. Notice finally that the sign of λ p is not fixed by (3.19). More generally, the sign of λ p depends on the (unphysical) phase of the nucleon wave function. The same holds for the overall sign of the form factors W n,α ΓΓ (s 0 , Q 2 ) that are determined from (3.17). The relative sign between W 0 ΓΓ (Q 2 ) and W 1 ΓΓ (Q 2 ) is however fixed by our sum rules. In the following, we will choose a negative sign for the coupling strength of the proton current, i.e. we will employ λ p < 0.

Numerical analysis
To derive physical predictions from (3.17), we need to find regions where the LCSRs converge sufficiently fast as an expansion in Λ 2 QCD /Q 2 and Λ 2 QCD /M 2 and where the sum rules are to a certain extent insensitive to the choice of the continuum threshold s 0 and the Borel mass M . Therefore, the Borel mass M has to be chosen well above the QCD scale Λ QCD but at the same time well below the mass of the lightest excitation. These conditions are formulated more precisely in the following, and they lead to a set of requirements which are then applied to each of the four LCSRs (3.17) as well as the local sum rule (3.19).
In order to eliminate contributions other than the proton in our sum-rule calculations [64], we use s 0 =s 0 = (1.44 GeV) 2 as a central value in all cases. This value of the continuum threshold s 0 corresponds to the mass of the lightest excited state in the nucleon spectrum, i.e. the Roper resonance. We then vary s 0 (ands 0 ) between (1.4 GeV) 2 and (1.5 GeV) 2 to estimate the uncertainty related to the choice of the continuum threshold. A similar procedure has been adopted in [65,66], and our choice can be further motivated by the observation that for values in this interval, the sum rule (3.19) leads to a good agreement with the LQCD results for λ p (see for instance [34,[67][68][69]).
A lower bound on M is determined by demanding sufficient suppression of higher powers in the OPE. In particular, we require that the contribution of the highest dimensional condensate in each LCSR does not amount to more than approximately 30% of the total QCD result. An upper limit on M is instead obtained by demanding that the ground-state (4.2) We then vary the Borel mass M in this so obtained Borel window to estimate the systematic uncertainty related to the variation of this unphysical parameter.
The physical values of the form factors do not depend on the choice of the continuum threshold s 0 or the Borel mass M . The residual dependence of the form factors extracted from the LCSRs on these parameters originates from the truncation of the expansion in Λ 2 QCD /M 2 at a finite order and the effective description of the a priori unknown contributions of heavy states. Therefore, the predictions of the sum rules are reliable if the dependence on the unphysical parameters is weak, and thus the uncertainties related to the variation of these parameters also quantifies the validity of the predictions.
In order to illustrate the latter statements we show in Figure 2 and  .17) is not present. This in turn leads to a stronger sensitivity on M such that the form factors increase again for larger Borel masses. Also the sensitivity to the continuum threshold is more pronounced if M gets closer to s 0 , as indicated by the widening of the coloured bands in Figures 2 and 3. One can furthermore observe that the sensitivity on the unphysical parameters becomes stronger for larger values of Q 2 . We remark that for Q 2 2 GeV this effect saturates such that the plots on the right-hand side of Figures 2 and 3 represent in a sense worst-case scenarios. The results of the other LCSRs behave similarly to W 0,S+T RR (s 0 , Q 2 ) and W 1,Q LR (s 0 , Q 2 ), so we do not show their dependence on M explicitly.
By considering the momentum range 0.5 GeV 2 ≤ Q 2 ≤ 2.5 GeV 2 , we find the Borel windows 0.7 GeV ≤ M ≤ 1.1 GeV for the LCSRs with α = S, Q, T and ΓΓ = RR and 0.7 GeV ≤ M ≤ 1 GeV for the LCSRs with α = S, Q, T and ΓΓ = LR. The LCSRs for the structure α = P do not meet the above requirements in the Q 2 region of interest, because the contributions of heavy states are large and even dominate the sum rules for certain values of Q 2 . In order to have one common Borel window for each chirality combination, we also choose in the case α = P either 0.7 GeV ≤ M ≤ 1.1 GeV or 0.7 GeV ≤ M ≤ 1 GeV as our Borel window when studying the Q 2 dependence of W 0,P ΓΓ (Q 2 ). In our numerical analysis of the LCSRs we use (m u +m d )/2 = (3.410±0.043) MeV [62] which corresponds to the MS value at 2 GeV. Using the two-loop renormalisation group (RG) running and the one-loop threshold corrections as implemented in RunDec [70,71], we obtain at 1 GeV the value m u + m d = (8.60 ± 0.11) MeV. Employing the GMOR relation this value leads to if the leading-order chiral corrections of [72] are included and uncertainties are added in quadrature. For the non-perturbative parameters defined in (3.14) we then find The parameter m 0 for the mixed condensate as well as the pure-gluon condensate are known from sum-rule estimates evaluated at 1 GeV. We will use the values and uncertainties from [73] which are widely accepted. The relevant numbers read We remark that using (4.5) the local QCD sum rule (3.19) agrees with the LQCD results for λ p within uncertainties (cf. Table I  After having explained how we choose the continuum thresholds and the Borel mass windows and having specified the numerical values of the input parameters, we are in a position to present the results of our LCSR analysis. Our results for the form factors W n,α RR (Q 2 ) and W n,α LR (Q 2 ) are shown in Figure 4 and Figure 5 as coloured lines and bands, respectively. The predictions in the range 0.5 GeV 2 ≤ Q 2 ≤ 2.5 GeV 2 result from a direct evaluation of (3.17). The solid curves correspond to the results obtained for the central values of the unphysical and physical parameters, while the bands reflects the corresponding theoretical uncertainties. The theoretical uncertainties are determined by varying all input parameters independently within their allowed ranges and adding individual uncertainties in quadrature. 3 For Q 2 ≤ 0.5 GeV 2 we instead rely on an extrapolation. Specifically, we consider both a linear and a quadratic fit in Q 2 to the LCSR form factors W n,α ΓΓ (Q 2 ) evaluated in the vicinity of Q 2 = 0.6 GeV 2 , and take the smallest and largest values of the fits at each Q 2 to obtain the displayed uncertainty bands. The results of the quadratic fit to our central form-factor predictions are indicated as dashed lines. For comparison we also show the values of W n ΓΓ (Q 2 ) determined in the recent LQCD study [39]. The numbers given in this article correspond to the MS form factors evaluated at 2 GeV and we use the two-loop RG running (cf. [33,74]) of (2.1) to evolve the form factors down to 1 GeV. The shown single (double) error bars represent the statistical (total) uncertainties of the LQCD predictions. Notice that the LQCD uncertainties at the physical point, i.e. Q 2 0, are dominantly of systematic origin.
As explained before, based on our study of the Borel windows we expect the LCSR prediction for W 0,P ΓΓ (Q 2 ) to be less reliable than the other results because of the large contributions of heavy states. Indeed, comparing the results of W 0,P ΓΓ (Q 2 ) and W 0,S+T ΓΓ (Q 2 ) as shown in Figures 4 and 5, one finds that W 0,S+T ΓΓ (Q 2 ) is closer to the LQCD predictions  than W 0,P ΓΓ (Q 2 ) for both chirality combinations, and that W 0,S+T LR (Q 2 ) itself agrees well with the LQCD calculation within uncertainties. One also observes from Figure 4 that the LCSR predictions for the modulus of W n,α RR (Q 2 ) tend to undershoot the LQCD results. An exhaustive comparison to the shown LQCD results for Q 2 0.5 GeV 2 would require knowledge about the systematic uncertainties of the LQCD calculations for non-zero Q 2 . A full error budget is in Tables 4 and 5  displayed LQCD predictions for Q 2 0.5 GeV 2 . The observed differences between the LCSR and the LQCD results may be related to higher-twist effects. In order to examine this issue, we have calculated the twist-4 corrections to (3.10), which is the only two-particle twist-4 correction [55], using the hadronic input parameters provided in [82]. We find that for Q 2 = 0.5 GeV 2 the relative corrections to the values of the form factors shown in Figure 4 amount to 38% for W 0,P RR , 5% for W 0,S+T RR , 51% for W 1,Q RR and 26% for W 1,T RR . Other twist-4 corrections to the LCSRs arise from additional three-particle DAs (see [60,61] for details), and depending on their size and sign the actual effect of twist-4 corrections may be notably different from the numbers quoted here. Nevertheless, the corrections we have computed are larger for the vectorial structures than for the scalar and tensor structure. This could explain why our LCSR calculation of W n,α RR (Q 2 ) seems to work better for α = S + T, T than for α = P, Q. n Figure 6. Comparison between the physical form factors W n RR (left) and W n LR (right) obtained by our LCSRs (red squares and error bars) and the state-of-the-art LQCD calculation (black dots and error bars) [39]. The shown results correspond to the MS scheme renormalised at 2 GeV. See main text for additional details.
As a comparison, the two-particle twist-4 contributions to the form factor values shown in Figure 5 amount to 22% for W 0,P LR , 7% for W 0,S+T LR , 12% for W 1,Q LR and 31% for W 1,T LR at Q 2 = 0.5 GeV 2 . In this case the tensor structure receives a larger correction than the vectorial structures, but overall the twist-4 corrections seem to be better under control for ΓΓ = LR than for ΓΓ = RR. This may explain why the LCSR predictions for W n,α LR (Q 2 ) are in general in good agreement with the LQCD results. In conclusion, we expect that uncertainties due to higher twist are minor for Q 2 1 GeV 2 , while in the range 0.5 GeV 2 Q 2 1 GeV 2 , twist-4 corrections may in the case ΓΓ = RR account for the differences between our LCSR predictions and the corresponding LQCD results. Notice that on general grounds one would expect that the total uncertainties of the LCSRs become larger for decreasing values of Q 2 , because the power suppression in Λ 2 QCD /Q 2 of the light-cone expansion (3.1) starts to becomes ineffective. In the plots of Figures 4 and 5 this effect is mimicked by our extrapolation procedure that leads to larger total uncertainties for Q 2 0.5 GeV 2 .
The physical form factors W 0 ΓΓ ≡ W 0 ΓΓ (Q 2 0) can be extracted from both the LCSR for W 0,P ΓΓ (Q 2 ) and W 0,S+T ΓΓ (Q 2 ), while in the case of W 1 ΓΓ ≡ W 1 ΓΓ (Q 2 0) one can consider the two independent combinations W 1,Q ΓΓ (Q 2 ) and W 1,T ΓΓ (Q 2 ). Since we believe that the LCSR for W 0,P ΓΓ (Q 2 ) is unreliable, we determine W 0 ΓΓ from the full range of solutions for W 0,S+T ΓΓ (0). The prediction for the form factor W 1 ΓΓ is instead obtained from the extrapolations leading to W 1,Q ΓΓ (0) and W 1,T ΓΓ (0), because in this case the different LCSR estimates result in quite similar numerical predictions (see Figures 4 and 5). At a renormalisation scale of 1 GeV, we obtain in this way the following central values and uncertainties:  Our LCSR predictions have total uncertainties of around (25 − 40)%. In Figure 6 we compare the results (4.6) and (4.7) evolved to 2 GeV to the corresponding LQCD pre-dictions [39]. Notice that two-loop RG effects (see [33,74]) lead to an enhancement of the LCSR results by 8.9% and 9.9%, respectively. From the two panels it is evident that while the LCSR approach does not achieve the (10 − 15)% accuracy of the latest LQCD computations of the form factors W n ΓΓ , the overall agreement between our LCSR predictions and the latest LQCD results is quite compelling.

Conclusions
In our work, we have calculated the hadronic matrix elements of the full set of baryonnumber violating dimension-six SMEFT operators (2.1) using LCSR techniques. These hadronic matrix elements are needed to predict the rates of the main proton decay modes in GUTs, where a proton decays into a pseudoscalar meson and an anti-lepton. Specifically, we have focused on the decay p → π 0 e + , and presented explicit LCSR expressions for the relevant form factors that include the leading contributions in the light-cone expansion, namely the twist-2 and twist-3 pion DAs (cf. Appendix C). We have performed a detailed study of the dependence of the LCSRs on both the unphysical (i.e. the continuum threshold and the Borel mass) and the physical (i.e. the condensates and the pion DAs) parameters, and discussed the possible impact of twist-4 effects. This enabled us 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, i.e. Q 2 = −q 2 > 0, and lies in the range 0.5 GeV 2 ≤ Q 2 ≤ 2.5 GeV 2 . We have then extrapolated our LCSR results to the physical point Q 2 0 by means of both a linear and quadratic fit, including the spread of predictions in our uncertainty estimates. Our analysis indicates that the LCSR for W 0,P ΓΓ (Q 2 ) is not reliable, and we therefore consider only W 0,S+T ΓΓ (Q 2 ), W 1,Q ΓΓ (Q 2 ) and W 1,T ΓΓ (Q 2 ) when determining the final predictions for the physical form factors W n ΓΓ with n = 0, 1 and ΓΓ = RR, LR from the range of different solutions shown in Figures 4 and 5.
Our final results for W n ΓΓ can be found in (4.6) and (4.7), and the LCSR results are compared to the state-of-the-art LQCD predictions [39] in Figure 6. The uncertainties of the LCSR results amount to (25−40)%, while the total accuracy of the LQCD form factors is (10 − 15)%. In view of the inherent systematic uncertainties of LCSRs, it is not clear to which extent possible refinements of our calculations such as including higher-twist contributions or perturbative corrections would allow to increase the precision of (4.6) and (4.7). The observed overall agreement between our and the latest LQCD form factors demonstrates however that LCSRs can be successfully applied to the calculations of proton decay matrix elements, and that such computations can achieve a precision that is better than the methods that have been developed in the '80s to estimate proton decay rates.
This gives us confidence that with the help of the LCSR techniques developed in this article it should be possible to obtain (at least) order-of-magnitude estimates for the hadronic matrix elements that appear in certain three-body proton decays. One possible application is the decay mode p → π 0 e + + G with G denoting a graviton. This channel is expected to be the dominant proton decay mode in the effective theory of gravity coupled to the SM aka GRSMEFT [44,45], since the two-body transition p → e + + G is forbidden by angular momentum conservation. LQCD calculations of three-body proton decay processes at arbitrary kinematics seem to be in reach in the coming years (see [47] for a discussion), but it remains to be seen which accuracy such computations can initially achieve. The calculation of hadronic matrix elements for three-body proton decay modes utilising LCSRs therefore seems to be a worthwhile undertaking, and our work provides the blueprints for such future studies.

Acknowledgments
AH would like to thank Javi Serra for useful discussions on the topic. The analytical calculations in this article were performed with the help of FeynCalc [75][76][77]. Some of the Dirac traces were cross-checked against Tracer [78]. The Feynman diagrams were created with the LaTeX package feynMF [79].

A Pion DAs
We use the following expressions for the pion DAs including terms proportional to the pion mass, which have been derived in [60] (and [55] in the chiral limit) with the help of a conformal expansion. One has where 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 presented in [80], where the moments are obtained by fitting sum rules for the electromagnetic pion form factor to the experimental data of [81]. For the numerical values of the other parameters we rely on the sum rules estimates of [82]: Using the definition f 3π (µ) ≡ f π µ π η 3 (µ) together with (4.4) we then find, where the individual uncertainties are added in quadrature.

B Compendium of analytic formulas
In the following we present a number of useful analytic formulas that enter our computations. As a first step in obtaining the LCSRs (3.17), we have to perform the Fourier transformation, which amounts to solving integrals of the type where the relevant momenta P ≡ q +ūp π and P g ≡ q + αp π with α ≡ α u + uα g arise from combining the exponential factor in (2.6) with those of (3.10) to (3.13). The momentum dependence can be rewritten in terms of Q 2 = −q 2 and P 2 p = −p 2 p using where we employ the notationx ≡ 1 − x for any variable x throughout this section. The UV divergent Fourier integrals are carried out in dimensional regularisation. In this step the poles and the scheme-dependent constants can be dropped in the sumrule calculation as long as we perform a Borel transformation in the end. Only inverse powers of the momenta, arising from finite integrals, and logarithms from the divergent integrals contribute to our sum rules. Moreover, only the imaginary parts of the correlation functions Π QCD,α ΓΓ (s + i , Q 2 ) enter the dispersion integrals (3.17), and as explained in the main text we subtract the heavy contributions for s > s 0 . We then perform a Borel transform, which is defined by for some function F (P 2 p ). In the final step the integration over s is performed. All the steps described above can be translated into certain replacement rules. For the logarithmic terms we find with µ the renormalisation scale and f (u) some function which depends on the momentum fraction u. Here we have introduced where the definition of E n (x) can be found in (3.20). The upper limit of the integration over u satisfies ∆(Q 2 = 0) = 1 = ∆(s 0 → ∞) and ∆ ≤ 1, and it arises because the dispersion integral only has support if s 0 ≥s. For terms involving three-particle DAs one furthermore has where f (u, α u , α g ) now depends on u, α u and α g , and we have eliminated We have furthermore used the abbreviation for the integration measure. The integration boundaries are now modified by the Heaviside step function θ(x), where ∆ g (Q 2 = 0) = 0 = ∆(s 0 → ∞) and ∆ g ≥ 0. For the non-divergent contributions appearing in our LCSRs, we find and for the three-particle integrals:  Dα θ (α − ∆ g ) T (3) (1 − α u − α g , α u , α g ) α 3 e(s g ) Dα θ (α − ∆ g ) T (3) (1 − α u − α g , α u , α g ) α 3 e(s g ) Dα θ (α − ∆ g ) T (3) (1 − α u − α g , α u , α g ) α 2 e(s g ) Recall that ζ = 2u − 1 and notice that we have used the definitions (B.6), (B.7), (B.9) and (B.10) to write the QCD correlation functions in a compact form. We have furthermore suppressed the renormalisation scale dependence of the pion DAs. The analytic expressions for the DAs are collected in Appendix A. Notice that since we have neglected quark-mass effects in (3.7) and (3.8), it would be consistent to set to zero all terms proportional to m 2 π in the formulas (C.1) to (C.8). While these contributions are in fact numerically small, it turns out that they always improve the agreement between the LCSR form factors calculated here and the LQCD form factors computed in [39]. We therefore included the m 2 π terms in the expressions provided above.