Systematic Improvement of x -dependent Unpolarized Nucleon Generalized Parton Distribution in Lattice-QCD Calculation

We present a first study of the effects of renormalization-group resummation (RGR) and leading-renormalon resummation (LRR) on the systematic errors of the unpolarized isovector nucleon generalized parton distribution in the framework of large-momentum effective theory (LaMET). This work is done using lattice gauge ensembles generated by the MILC collaboration, consisting of 2+1+1 flavors of highly improved staggered quarks with a physical pion mass at lattice spacing a ≈ 0 . 09 fm and a box width L ≈ 5 . 76 fm. We present results for the nucleon H and E GPDs with average boost momentum P z ≈ 2 GeV at momentum transfers Q 2 = [0 , 0 . 97] GeV 2 at skewness ξ = 0 as well as Q 2 ∈ 0 . 23 GeV 2 at ξ = 0 . 1, renormalized in the MS scheme at scale µ = 2 . 0 GeV, with two-and one-loop matching, respectively. We demonstrate that the simultaneous application of RGR and LRR significantly reduces the systematic errors in renormalized matrix elements and distributions for both the zero and nonzero skewness GPDs, and that it is necessary to include both RGR and LRR at higher orders in the matching and renormalization processes.


I. INTRODUCTION
An open question in the theory of quantum chromodynamics (QCD) is how the fundamental degrees of freedom, quarks and gluons, comprise the more massive hadrons.The quarks and gluons (known collectively as partons) contribute to a hadron's mass and spin but cannot be studied in isolation due to confinement.Thus, knowledge of the internal structure of a hadron is highly valued.Great effort has been focused on the study of parton distribution functions (PDFs), which describe the distribution of a hadron's longitudinal momentum among its constituents, and much has been learned about hadronic structure from these studies (see Ref. [1] for a review from Snowmass 2021).However, the PDF only paints a one-dimensional picture of the hadron, since it is dependent solely on longitudinal momentum.Generalized parton distributions (GPDs) contain more information about the hadron, including spin structure, form factors and how the longitudinal momentum of the parton depends on the distance from the center of the hadron.The unpolarized GPD is comprised of two functions commonly denoted H and E, defined in terms of matrix elements on the lightcone as where L(−z/2, z/2) is a link along the lightcone, Q µ = (p ′′ − p ′ ) µ is the momentum transfer, and ξ = p ′′+ −p ′+ p ′′+ +p ′+ is the skewness.In the limit Q 2 → 0 and ξ → 0, the H GPD reduces to the PDF.The E GPD is inaccessible in this limit, since it is multiplied by the momentum transfer vector.GPDs can be probed experimentally by processes such as deeply-virtual Compton scattering (DVCS) [2,3] or deeply-virtual meson production (DVMP) [4], and their study will be an important experimental program at the future Electron-Ion Collider (EIC) [5][6][7][8][9].Lattice QCD involves converting the QCD path integral from continuous Minkowski spacetime to discrete Euclidean spacetime, making field-theory calculations amenable to supercomputers.It can provide early in-sight into GPD functions complementary to experimental programs.The computation of the Bjorken-x dependence of parton distributions can be studied in the framework of lattice QCD using one of several recent methods: the "hadronic-tensor approach" [10][11][12][13][14][15], the Comptonamplitude approach (or "OPE without OPE") [16][17][18][19][20][21][22][23][24][25][26][27][28], the "current-current correlator" method [23,[29][30][31][32][33][34][35], or large-momentum effective theory (LaMET) [36][37][38], which is our focus in this paper.
The method of LaMET begins with the study of spatially separated, equal time, matrix elements of boosted hadrons computed directly on the lattice: where Γ = γ t , γ t γ 5 , γ t γ ⊥ for the unpolarized, helicity and transversity GPDs, respectively.W (−z/2, z/2) is a lattice link from the coordinate (0, 0, 0, −z/2) to (0, 0, 0, z/2), since we may assume without loss of generality that the average momentum ( ⃗ p ′ + ⃗ p ′′ )/2 is along the z axis.The bare matrix elements are then renormalized and Fourier transformed to momentum space to obtain the quasi-GPD.The final step is to match the quasi-GPD to the lightcone to obtain the GPD.GPDs have been studied in LaMET in the Breit-frame setting on the lattice.The GPD on the lattice was first studied in the case of the pion in Ref. [39] and carried out at physical pion mass [40,41] by MSULat in the zero-skewness limit.The nucleon unpolarized and helicity GPDs were studied in Refs.[42][43][44] and the transversity ones in Ref. [45].Recently, the ETMC and ANL/BNL groups have computed bare matrix elements in asymmetric frames [46] to help reduce the computational cost of the lattice calculation.
Since the aforementioned numerical studies of GPDs, developments in the framework of LaMET include renormalization-group resummation (RGR) [47] and leading-renormalon resummation (LRR) [48].RGR is designed to resum the logarithms that arise from the differing intrinsic physical scale and final renormalization scale of the parton.The method is to set the energy scale such that the logarithmic terms vanish and then evolve to the desired scale with the renormalization group.This process can be applied both to the renormalization of the bare matrix elements as well as the perturbative matching.LRR is designed to resum the divergence arising from the infrared renormalon (IRR) which plagues perturbation series [49], and whose effect is more pronounced with the application of RGR alone.The first application of LRR was to the pion PDF in Ref. [48], which showed that LRR in combination with RGR results in greatly reduced systematic uncertainties in the final x-dependent PDF.The ANL/BNL collaboration also applied LRR (and RGR) to their LaMET calculation of the nucleon transversity PDF in Ref. [50] to better control the systematic errors.The field of LaMET has matured to the point at which such systematic uncertainties become an important issue.The methods of RGR and LRR have not yet been applied to GPDs; doing so can lead to more precise calculation of tomography from lattice QCD in the future.
The purpose of this paper is to make the first application of the RGR and LRR improvements to the calculation of the unpolarized nucleon isovector GPD at different skewness, ξ, and squared momentum transfer, Q 2 , in the Breit frame.We use clover valence fermions at physical quark mass with a lattice spacing of a ≈ 0.09 fm and box length L = 64a ≈ 5.76 fm with QCD vacuum composed of N f = 2 + 1 + 1 flavors of highly-improved staggered quarks [51], generated by the MILC collaboration [52][53][54] with one step of hypercubic smearing [55] applied to the gauge links to reduce discretization effects.The valence fermion parameters are tuned so as to produce a physical pion mass (m π ≈ 130 MeV).The same mixed-action setup used in this calculation was previously studied in Refs.[56][57][58][59][60][61][62][63][64][65][66][67][68][69][70] and found to be free of exceptional configurations which can cause the Dirac matrix to be ill-conditioned or the correlation functions to be anomalously large.From a total of 1960 lattice configurations, we use the 501,760 measurements of the bare nucleon matrix elements of Eq. 2 with average boost momentum P z = 10× 2π L ≈ 2.2 GeV in Ref. [42].More information on the bare matrix elements such as the sourcesink separation, the momentum smearing and momentum transfer can be found in Ref. [42] and its supplemental material.The ground-state nucleon bare matrix elements are extracted by simultaneously fitting multiple source-sink separations with skewness values of ξ = 0 and ξ = 0.1.For each skewness value, we have momentum transfer Q 2 ∈ {0.0, 0.19, 0.39, 0.77, 0.97} GeV 2 and Q 2 = 0.23 GeV 2 respectively.This paper is laid out as follows.In Sec.II, we describe the methodology of RGR and LRR as well as the outline of our calculation of the GPDs at zero skewness from the bare matrix elements.We also show results for zeroskewness GPDs for both zero and nonzero momentum transfer, demonstrating the improvements afforded by both RGR and LRR as well as matching at both next-toleading-order (NLO) and next-to-next-to-leading-order (NNLO).In Sec.III, we show nonzero-skewness GPDs at NLO.We conclude in Sec.IV.

II. ZERO-SKEWNESS GPDS AT NLO AND NNLO
In this section we present the zero-skewness (ξ = 0) unpolarized isovector nucleon GPD at both zero (H GPD only) and nonzero (H and E GPDs) momentum transfer Q 2 .When both ξ = 0 and Q 2 = 0, the unpolarized GPD reduces to the unpolarized PDF.The renormalization procedure and the transformation to momentum space are also described in this section, since the same methods are used for all values of momentum transfer and skewness.We describe the lightcone matching for the case of zero skewness and postpone the discussion of nonzero skewness matching to Sec.III.
We begin with the renormalization of the bare matrix elements.We perform the renormalization in the hybrid-ratio scheme [71], in which the bare matrix elements are renormalized in the ratio scheme up to distances z s = 3a ≈ 0.27 fm with our lattice spacing, and at large distances the linear divergence and renormalon divergence are removed by an exponential term.The ratio scheme involves dividing the bare matrix element at nonzero boost momentum by those at zero boost momentum at fixed z.The fully renormalized matrix element (for both H and E) is given by (3) where we have used bare unpolarized pion matrix elements at zero boost-momentum, h B π (z, P z = 0) [72] for the ratio scheme at z < z s .At Q 2 = 0, we normalize the matrix elements to 1 at z = 0.The terms δm and m 0 are, respectively, the linear divergence and the renormalon divergence.The linear divergence is due to the self-energy of the Wilson line in the bare matrix element, and the renormalon divergence arises from the fact that the perturbation series used to calculate δm is not convergent to all orders [48,49,71,72].We determine the linear divergence by following the same procedure as in Ref. [71] by fitting the zero-momentum pion matrix elements to the exponential decay Be −δm z in the interval z = [0.54,1.53] fm, as shown in the left-most panel of Fig. 1, where B and δm are fitting parameters.This same procedure was performed with the same data in our previous work [73] in which we find δm = 0.668 (10) GeV.
While the computation of the linear divergence would seem to be subjective, it is compensated for by the calculation of the renormalon divergence such that their sum, δm + m 0 , is constant in a fixed scheme [71].The renormalon divergence is determined by demanding that the short-distance physics (z ≲ 0.3 fm) agrees with the theoretical predictions of the operator-product expansion (OPE).The functions that appear in the OPE (and describe the short distance physics) are known as Wilson coefficients, which we denote by C 0 (z, µ), where z is the Wilson length, and µ is the energy scale.For a Wilson length z and renormalization scale µ, which is the final desired energy scale for the lightcone PDF renormalized in the modified minimal-subtraction (MS) scheme, the unpolarized Wilson coefficients are at NLO [74] and at NNLO [75] where l(z, µ) = ln z 2 µ 2 e 2γ E /4 , γ E is the Euler-Mascheroni constant, α s (µ) is the strong coupling at energy scale µ, C F is the quadratic Casimir for the fundamental representation of SU(3) and n f is the number of fermion flavors.
The determination of the renormalon divergence can be improved with the additions of renormalization-group resummation (RGR) [47,48] and leading-renormalon resummation (LRR) [48].The difference between the intrinsic physical scale and the final renormalization scale results in logarithmic terms that require resummation.This is achieved by setting the renormalization scale such that the logarithmic terms vanish and then evolving to the desired scale using the renormalization-group equation: where γ(µ) is the anomalous dimension, which has been calculated up to three loops [76].The energy scale at which the logarithms vanish is µ = 2e −γ E /z ≡ z −1 as can be seen in Eqs. 4 and 5. Thus, we can improve the computation of the Wilson coefficient with RGR giving where k = 1 for NLO, k = 2 for NNLO and β(α) is the QCD beta-function.For brevity, we define The Wilson coefficients are perturbation series which can contain a renormalon divergence [77].We account for this using the LRR method, in which the Wilson coefficient is modified according to Eq. 14 of Ref. [48] C where r i are the coefficients of the renormalon series in α s and C PV (z, µ) is the contribution of the renormalon to the Wilson coefficients after a Borel transformation originally derived in Refs.[78,79].Explicit definitions can be found in Eqs. 12 and 13, respectively, of Ref. [48] We can then combine the RGR and LRR improvements into a single high-quality Wilson coefficient The Wilson coefficients with different improvements yield different central values and uncertainties for the renormalon divergence, m 0 .We use the same procedure to compute the renormalon divergence as Ref. [48] in which is fitted to m 0 z +c for multiple sets of z values.We interpolate the matrix elements h B π (z, P z = 0) as in our previous work [73] and determine m 0 (z) with the inputs of {z − 0.02 fm, z, z + 0.02 fm} to a maximum of z = 0.2 fm.A plot of m 0 to different orders and with different improvements as a function of fitting range is shown in the middle and right panels of Fig. 1.We seek a plateau in the values of m 0 across different fitting ranges which signals a stable and reliable measurement of the renormalon divergence and select the corresponding value as the measurement of m 0 .The results with the smallest errors as well as clear plateaux are those for which RGR and LRR are applied simultaneously.Having determined both the linear divergence and the renormalon divergence, we now have fully renormalized matrix elements in the hybrid-ratio scheme (using Eq. 3).
To obtain the quasi-distribution, we first extrapolate the renormalized matrix elements to infinite distance with a view to performing a Fourier transform.The extrapolation model is inspired by the small-x physics we expect to see in the PDF [71,80,81], which is itself governed by the large-distance behavior of the renormalized matrix elements: where A, m and d are fitting parameters.The data used to fit the extrapolation must be at sufficiently large z that we can realistically model the large-distance behavior.
We then Fourier transform to momentum space to obtain the quasi-GPDs with the convention where F is either H or E corresponding to the respective GPD functions.By extrapolating the renormalized matrix elements to infinite distance, we remove unphysical oscillations from the quasi-GPDs that would otherwise occur in the Fourier transform.The final stage in the calculation is the perturbative matching to align the ultraviolet (UV) behavior of the quasi-GPD with the lightcone.The matching formula is where K is the matching kernel.Once again, this formula applies to both the quasi-H and quasi-E GPDs.For zero skewness, ξ = 0, the kernel has been calculated up to NNLO in the hybrid-ratio scheme for unpolarized GPDs in Refs.[47,75,82].For nonzero skewness, the kernel has been computed up to NLO for unpolarized GPDs (as well as helicity and transversity GPDs) [83], and we discuss it in more detail in Sec.III.The RGR process applied to the matching is designed to resum logarithmic terms that occur in the matching kernel.The philosophy is the same as that of the determination of the renormalon divergence with RGR in that we set an energy scale such that the logarithms vanish and then evolve to the final desired energy scale.This time, the anomalous dimension is the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation where P(z) is the DGLAP kernel, which has been calculated up to three loops [84].The formula is applicable to both H and E GPDs.We use the same algorithm for RGR matching as in Ref. [47].However, this formula is only applicable to zero-skewness GPDs.At nonzero skewness, a different evolution formula is required for |x| < ξ.The corresponding formula is the Efremov-Radyushkin-Brodsky-Lepage (ERBL) equation [85][86][87][88] and |x| < ξ is known as the ERBL region.In this x range there are two distinct scales that emerge, which cannot be eliminated simultaneously by the choice of a single energy scale.A more sophisticated technique must be developed for this case in the future.
We begin by looking into the special case of the nucleon unpolarized GPD at Q 2 = 0 and ξ = 0, which is equivalent to the PDF.We use the four different methods of computing the renormalon divergence m 0 at NLO and again at NNLO with the renormalized matrix elements h R H (z, ξ = 0, Q 2 = 0).Our notation for the different schemes is "(N)NLO×RGR" for the RGR improvement only, "(N)NLO+LRR" for the LRR improvement only, and ((N)NLO+LRR)×RGR for both the RGR and LRR improvements. 1 We show the real and imaginary parts of the matrix elements for (N)NLO in the top (bottom) of Fig. 2. The (N)NLO, (N)NLO×RGR, (N)NLO+LRR and ((N)NLO+LRR)×RGR matrix elements are plotted in blue, red, green and purple, respectively.Except for (N)NLO, the data points are offset slightly from their true z values to allow for readability.The plots contain both statistical error bars and combined statistical and systematic error bars from scale variation.In the case of the renormalized matrix elements, the systematic errors are computed by scale variation as was used in Ref. [48].When RGR is applied to the Wilson coefficients, we vary the initial energy scale used in the RGR process, c ′ ×z −1 , before evolving to the final desired one.The central value corresponds to c ′ = 1.0; the upper and lower error bars are derived by varying c ′ from 0.75 to 1.5.The range c ′ ∈ [0.75, 1.5] corresponds to a change of approximately 15% on either side of α s (µ = 2.0 GeV).This creates two additional curves with the maximum (minimum) value corresponding to the upper (lower) systematic error.The systematic errors are asymmetric, since the strong-coupling dependence on energy scale is nonlinear.When RGR is not applied to the Wilson coefficients, the systematic errors are determined by computing the renormalon divergence at energy scales 0.8 GeV and 2.8 GeV with 2.0 GeV being the central value.These scale variations yield different measurements of the renormalon divergence and the upper-and lower-values are interpreted as the upper and lower systematic errors, respectively.The RGR and LRR improvements to the Wilson coefficients give different central values and uncertainties in the renormalon divergence, resulting in different systematic errors in the renormalized matrix elements.
Examining the four NLO schemes in the top row of Fig. 2, we can see that the relative systematic errors are reduced by approximately 15% to 35% from NLO to NLO+LRR.The reduction from NLO×RGR to (NLO+LRR)×RGR, however, is approximately 70% to 90% showing that leading-renormalon resummation has a much greater effect when used in combination with RGR.This is to be expected, since the Wilson coeffi-cients used to compute the renormalon divergence m 0 are series expansions in the strong coupling α s , and the renormalon divergence does not emerge until we expand the series to a power n in the strong coupling where n ∼ 1/α s (µ) [49,89].At our smallest energy scale used at fixed order, µ = 0.8 GeV, α s (µ) ≈ 0.5, and the renormalon divergence will not emerge unless we expand beyond quadratic terms in the strong coupling; however, this does not mean that the renormalon divergence is irrelevant.We can see that there is an increase up to fifteenfold in the absolute systematic errors from NLO to NLO×RGR, since the latter does not account for the renormalon divergence.When we compute the Wilson coefficients (and hence the renormalon divergence) at NLO×RGR, we set the initial energy scale to µ = z −1 .At small z, this is a large energy scale, which results in a small α s , meaning that the renormalon divergence does not emerge at NLO×RGR in the series expansion.The opposite occurs at large z and, hence, the renormalon divergence can emerge at NLO×RGR.This divergence is passed on to the calculation of the renormalon divergence, resulting in large systematic errors, particularly at large z, where the renormalon divergence occurs sooner in the series expansion.For this reason, there is a significant difference between NLO and NLO×RGR.This reasoning also applies at NNLO, in fact, to a greater extent, as can be seen in the bottom half of Fig. 2.
With the above eight sets of renormalized matrix elements, we then construct the quasi-distributions.First, we take each set of the real and imaginary renormalized matrix elements at large Wilson-line displacement and extrapolate them to infinite distance using Eq.11.Here, we select the range z ∈ [8a, 15a] = [0.72,1.35] fm for both the real and imaginary parts.For all schemes as well as both real and imaginary parts, the χ 2 /dof values are less than 1, which indicates the extrapolation formula is good model for unpolarized GPD matrix elements.We construct a renormalized matrix element as a full function of z by making a piecewise function.At small z, we interpolate the lattice data, and at large z, we use the extrapolation model with the best-fit parameters.
We then Fourier transform our full function into momentum space using Eq. 12 to obtain the quasi-PDF and finally match to the lightcone using Eq. 13.When RGR is not included in the calculation, the matching is performed at fixed order; that is, K is evaluated at a fixed energy scale µ.When we include RGR, we perform matching at the energy scale µ = 2xP z , which removes the large logarithms in the kernel, and then evolve to the desired scale with the DGLAP formula in Eq. 14.The DGLAP equation begins to break down for |x| ≲ 0.2, since the strong coupling α s (µ = 2xP z ) becomes nonperturbative in this region.Hence, we do not plot the unpolarized lightcone GPD data within this region and shade it in light gray.In addition, the LaMET expansion breaks down for small-and large-|x| as in the matching formula in Eq. 13; we approximate the region where these corrections become greater than or equal to one and shade these regions in dark gray.
The systematic errors for the unpolarized PDFs are computed by renormalizing the bare matrix elements with the upper and lower values of the renormalon divergence given by varying the scale.We then perform the large-distance extrapolation, Fourier transformation and matching on the matrix elements.When RGR matching is used, we set the initial scale to µ = c ′ × 2xP z with the central value corresponding to c ′ = 1.0 and the systematic error bands coming from c ′ = 0.75 and 1.5 as in the determination of the renormalon divergence.This gives us a central value for the PDF as well as two other values which correspond to the upper and lower systematic errors.
In Fig. 3 we show the lightcone unpolarized GPDs in the "PDF limit" (Q 2 = 0 and ξ = 0) with statistical errors (inner error bands) and combined statistical and systematic errors (outer error bands).Since we have computed the unpolarized GPD, the regions x > 0 and x < 0 correspond to the combinations ξ) ("antiquark region"), respectively.The top (bottom) row shows the PDFs at (N)NLO.The left column shows no modifications and LRR only, the right column shows the RGR modification only and both RGR and LRR.We plot the (N)NLO, (N)NLO×RGR, (N)NLO+LRR and ((N)NLO+LRR)×RGR PDFs in blue, red, green and purple, respectively.In Eq. 13, there are corrections to the lightcone GPD that are suppressed with P z but grow at finite P z as x → 0 or |x| → 1. We, therefore, shade in the regions at small and large |x|, where the LaMET calculation breaks down.
Examining the x-dependent GPDs, we consider first the four NLO schemes (top row of Fig. 3).The statistical errors are more or less constant across the four of them, since the bare matrix elements are the same.It is clear that the systematic errors are at their minimum when both the LRR and RGR improvements are applied simultaneously.Indeed, much of the behavior of the systematic errors we see in Fig. 2 for the renormalized matrix elements also occurs in the PDFs.Examining the large-x region, we see that the four schemes become compatible with zero as x → 1 within one to two sigma.Across the quark region as a whole, we see that the central values across the four schemes are in general agreement for x ≳ 0.3 and the main difference is the variation in systematic errors.We anticipate this result from the fact that the renormalized matrix elements differ in the renormalon divergence and its uncertainty; the m 0 parameter in the four schemes are all compatible, but the error bars differ a great deal from one scheme to another.The antiquark region, given the fluctuations across the different schemes, is compatible with zero.Larger boost momenta will be required to improve the antiquark signal, as was demonstrated in Refs.[90][91][92][93].Hereafter, our focus will be on the quark region.
Turning next to the four NNLO results in the bot-tom half of Fig. 3, we see once again that the smallest systematic errors occur with (NNLO+LRR)×RGR.The systematic errors are approximately the same for NNLO and NNLO+LRR, as we would expect from the renormalized matrix elements in Fig. 2 having similar systematic errors for the two schemes.In the quark region, there is, in fact, little difference between the central values at (NLO+LRR)×RGR and (NNLO+LRR)×RGR except in the endpoint regions; however, going to higher order reduces the systematic errors.We note that while the (NNLO+LRR)×RGR scheme has the smallest systematic errors, the central value remains consistent with both NNLO and NNLO+LRR.The central values of the NNLO×RGR results differ from the other three NNLO results due to the enhancement of the renormalon divergence when RGR is applied on its own.Our results have shown that much of the advantage due to renormalization-group resummation and leadingrenormalon resummation comes from a reduction in the systematic errors computed from scale variation.The improved systematic errors with these schemes show that the benefits are transferable across different LaMET calculations, since their effects were first demonstrated in the case of the pion PDF [47,48] and pion DA [94].Given the significant differences in our (N)NLO×RGR and ((N)NLO+LRR)×RGR results and errors, we have shown that the renormalon divergence is a source of systematic errors that cannot be ignored as an esoteric phenomenon.

B. Zero-Skewness H and E GPDs at
In this section, we examine our results for both the unpolarized zero-skewness H and E GPDs at nonzero momentum transfer.The range of momentum transfer values used in this calculation is Q 2 ∈ {0.19, 0.39, 0.77, 0.97} GeV 2 .We start this subsection by showing an example of the renormalized matrix elements at the intermediate value Q 2 = 0.39 GeV 2 to demonstrate the effects of LRR and RGR on the calculation.The same procedures are applied to all of our ξ = 0 GPD functions at all momentum transfers.Since we have already studied the effects of NLO, NNLO and the applications of LRR and RGR in Sec.II A, we do not show every case here but restrict ourselves to NLO, NNLO, (NLO+LRR)×RGR and (NNLO+LRR)×RGR.
In Fig. 4, we show the real (left column) and imaginary (right column) renormalized matrix elements for the H and E GPDs, h R H (top row) and h R E (bottom row), at zero skewness and Q 2 = 0.39 GeV use the same fitting range for the large-distance extrapolation as was used in the PDF case (Q 2 = 0).However, at Q 2 = 0.77 and 0.97 GeV 2 , we use the fitting range z ∈ [11a, 15a] = [0.99,1.35] fm, since the h R E matrix elements change sign at larger range for this momentum transfer, and such behaviour cannot be accommodated by the extrapolation model in Eq. 11.
As in the case of the renormalized matrix elements at Q 2 = 0 in Fig. 2, we see a significant decrease in systematic errors from NLO to (NNLO+LRR)×RGR (30% to 70%) and an even greater decrease from NNLO to (NNLO+LRR)×RGR (70% to 90%) in Fig. 4.This is to be expected, since the systematic errors of the renormalized matrix elements are governed by the renormalon divergence, which is itself determined by the Wilson coefficients.The same benefits afforded by RGR and LRR that we see in Fig. 2 should occur at Q 2 = 0.39 GeV 2 , since the same Wilson coefficients are used and improved in the same ways.In addition, the systematic errors increase from NLO to NNLO as in the Q 2 = 0 case as we would expect from the behavior of the renormalon divergence.
In Fig. 5, we show the unpolarized H and E GPDs in the NLO, NNLO, (NLO+LRR)×RGR and (NNLO+LRR)×RGR cases for Q 2 = 0.39 GeV 2 .The inner error bars are statistical and the outer error bars are combined statistical and systematic errors the latter computed in the same way as in the Q 2 = 0 case in Sec.II A. As in the PDF case shown in Fig. 3, the systematics are at a minimum in the (NNLO+LRR)×RGR scheme both for H and E GPDs.The upper and lower systematic errors increase from NLO to NNLO for almost the whole interval x ∈ [0.2, 0.8] which shows that the need to account for both the large logarithms and the renormalon divergence persists across different Q 2 values.Also, the systematic errors decrease by up to 40% from (NLO+LRR)×RGR to (NNLO+LRR)×RGR in the interval x ∈ [0.3, 0.9] both for H and E GPDs.This shows the benefits of going up to two loops in the matching process.Once again, the central values for all four schemes are in general agreement, showing that the main improvement afforded by RGR and LRR is a reduction in systematic errors.It is also evidence for convergence in the matching procedure, since the central values for (NLO+LRR)×RGR and (NNLO+LRR)×RGR are close.It is to be expected that the improved systematic errors persist across Q 2 values since the RGR and LRR improvements are universal and should be applicable in all LaMET calculations.The fact that the systematics increase from NLO to NNLO and decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR, shows again that the handling of systematic uncertainties must keep pace with higher FIG.3: Isovector nucleon lightcone PDFs at NLO (top row) and NNLO (bottom row) without improvement (blue bands), with LRR only (green), with RGR only (red) and with both LRR and RGR (purple) improvements.The dark-gray regions are the x values at which the LaMET calculation breaks down.In addition, when RGR is applied (right column), the matching formula breaks down for |x| ≲ 0.2, so this region is shaded in light gray.
orders in the matching and renormalization processes.
These are the first applications of the RGR and LRR improvements to the LaMET calculation of the unpolarized nucleon GPD as well as the first application of hybrid-ratio renormalization to the same.We plot both the H and E GPDs for Q 2 values from 0.19 to 0.97 GeV 2 (as well as the H GPD for Q 2 = 0) at both (NLO+LRR)×RGR and (NNLO+LRR)×RGR in Fig. 6.The left (right) column corresponds to the unpolarized H (E) GPD.The top (bottom) row corresponds to ((N)NLO+LRR)×RGR.Once again, the inner bands correspond to statistical errors and the outer bands correspond to combined statistical and systematic errors computed as in Sec.II A. We see that both the H and E GPDs decrease with Q 2 , as has been seen in previous calculations of nucleon GPDs [42,44].In all cases, the systematic errors are greatly reduced once again by the simultaneous additions of RGR and LRR.This is more evidence for the universality of the renormalizationgroup resummation and leading-renormalon resummation.
The central values decrease from the quasi-GPD to the (NLO+LRR)×RGR GPD across all Q 2 and again when going from (NLO+LRR)×RGR to (NNLO+LRR)×RGR.This is to be expected as the matching process tends to decrease the GPD value in the mid-to large-x regions and increase the value at small-x.This is due to the probability of a parton carrying a high momentum-fraction decreasing as the hadron approaches the lightcone.However, with the application of RGR in the matching, we cannot reliably study the small-x region, |x| ≲ 0.2.Nevertheless, this first application of the RGR and LRR methods to GPDs at nonzero momentum transfer is a step toward precision GPDs from lattice QCD.

III. NONZERO-SKEWNESS GPDS
In this section, we show the results for GPDs evaluated at nonzero skewness.While the LRR method is directly transferable to ξ ̸ = 0, the RGR matching is not.In x space, the GPD is often broken down into two regions: the DGLAP region for |x| > ξ and the ERBL region for |x| < ξ.While the DGLAP evolution in Eq. 14 is applicable to the corresponding region, a different scaling formula is required in the ERBL region, and there is the additional issue of two different intrinsic scales, which cannot be eliminated simultaneously by the judicious selection of a single initial energy.For this reason, we only examine the effects of RGR on the renormalized matrix elements and confine our attention to the NLO and NLO+LRR GPDs in momentum space.
We start by looking at the renormalized matrix elements for both h R H and h R E at Q 2 = 0.23 GeV ξ = 0.1 with statistical errors (inner bars) and combined statistical and systematic errors (outer bars) for all NLO four schemes in Fig. 7 with the outer systematic error bars computed as detailed in the previous sections.The Wilson coefficients for the nonzero skewness are the same as those in zero skewness case.Although we cannot yet apply RGR matching at nonzero skewness, we can see that the systematic errors in the renormalized matrix elements follow the same pattern as in the zero-skewness cases in Figs. 2 and 4.This is evidence that the same improvements in the matching process adjusted for nonzero skewness should be equally effective as at ξ = 0.The matching kernel of the nonzero-skewness GPDs, K ξ , differs from the zero-skewness one used in the ,ξ,μ) ,ξ,μ) ,ξ,μ) Eq. ( 13) due to the fact that skewness parameter ξ encapsulates the change in the struck hadron's longitudinal momentum.To date, the ξ ̸ = 0 matching kernel, K ξ , has only been computed up to NLO for unpolarized GPDs in the hybrid-ratio scheme in Refs.[83,[95][96][97]; however, the kinematic setup of the kernels in Refs.[95][96][97] in the ERBL region are incomplete.For this work, we adopt the ξ ̸ = 0 matching kernel, K ξ , from Ref. [83]: Note that we modified the kernel to convention in which there is an extra factor of 1/|y| in the integrand, whereas Ref. [83] absorbed this factor into the kernel itself.The nonzero skewness matching kernel K ξ contains singularities at y = 0, y = x and |y| = ξ.It is invariant under ξ → −ξ and recovers the NLO zero-skewness kernel when taking limit ξ → 0. The LRR modification to the match-ing kernel is the same at both zero and nonzero skewness [48].The LRR matching modification is derived from the LRR modification to the Wilson coefficients.
Since the Wilson coefficients are the same for zero and nonzero skewness, the same modification to the matching kernel is applicable.
The final unpolarized H and E GPDs are shown in FIG.7: Real (left column) and imaginary (right column) renormalized h R H (top row) and h R E (bottom row) matrix elements at Q 2 = 0.23 GeV 2 and ξ = 0.1.We show data with NLO (blue), NLO+LRR (red), NLO×RGR (green) and (NLO+LRR)×RGR (purple) improvements.The solid error bars are statistical and the dashed error bars are combined statistical and systematic, the latter arising from the scale variation.Except for NLO (real and imaginary for both h R H and h R E ), the data points shown in the plots have been offset from their exact z value to allow for readability.
Fig. 8 for ξ = 0.1 at NLO and NLO+LRR; we plot vertical dashed lines at x = ±ξ.As in the zero-skewness case, there is little change between the two aforementioned schemes in central values or error bars.This is expected from the fact that the renormalon divergence has a lesser effect at fixed order (NLO) than it does when RGR is included (NLO×RGR).In addition, the GPD suffers a discontinuity at x = ±ξ due to the corresponding singularities in the matching kernel.One difference between our nonzero-skewness H GPD and those in Ref. [44] is that our H GPD does not plateau in the ERBL region.The unpolarized H GPD at ξ = 0.3 in Fig. 3 of Ref. [44] is approximately flat in the region |x| < ξ = 0.3, whereas our H GPD at ξ = 0.1 increases in the region |x| < ξ = 0.1.Our ERBL region lies within the x-range where the LaMET expansion breaks down.For this reason, we should perhaps not expect our calculation to have the same qualitative behavior as that of Ref. [44].The effect of LRR on the x-dependent GPD without RGR is similar to the corresponding effects at zero skewness (Fig. 3); we, therefore, anticipate that the improvements we see for the unpolarized GPDs at zero skewness will also manifest at nonzero skewness once the methods have been adapted for the latter.Because we are ultimately interested in the x dependence, this is an auspicious indication of the benefits of RGR and LRR at ξ ̸ = 0.

IV. CONCLUSION AND OUTLOOK
In this paper, we have shown the first application of leading-renormalon resummation and renormalizationgroup resummation to the unpolarized nucleon isovector GPD computed on the lattice in the framework of largemomentum effective theory.We used a lattice spacing a ≈ 0.09 fm with a physical pion mass, N f = 2 + 1 + 1 flavors of highly-improved staggered quarks and an average boost momentum P z ≈ 2.2 GeV with ensembles generated by the MILC collaboration [52][53][54].These matrix elements were renormalized in the hybrid-ratio scheme, applying RGR and LRR.We then extrapolated the renormalized matrix elements to infinite distance and Fourier transformed to momentum space.We report zero-skewness unpolarized nucleon GPDs, H and E, with multiple momentum transfer values Q 2 , which have been matched to two loops as well as improved with both RGR and LRR for the first time.The main advantage of the (NNLO+LRR)×RGR calculation over other schemes is the reduction in systematic errors, since the central values remain compatible between the four schemes as shown in Sec. 2. We also reported GPD functions ξ = 0.1 at a single momentum transfer value of Q 2 = 0.23 GeV 2 in this work.However, only the LRR improvement is applied to matching process up to one loop due to the lack an RGR calculation for nonzero-skewness GPDs to date.The LaMET systematic errors were greatly reduced by the simultaneous application of RGR and LRR in the renormalization and matching processes.For both the renormalized matrix elements and the x-dependent GPDs, the statistical errors remain approximately constant with the RGR and LRR modifications.The improved systematics persist in the determination of the x-dependent GPDs.The fact that systematic errors increase when we go from NLO to NNLO but decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR show that the handling of systematics must keep pace with higherorder expansions in the matching and renormalization processes.In addition, the systematic errors increased when RGR was applied on its own, due to its enhancement of the renormalon divergence.The application of RGR and LRR to multiple Q 2 values at ξ = 0 showed the efficacy of the two processes for nonzero momentum transfer.Finally, we showed that the effects of RGR and LRR on the renormalized matrix elements and GPDs at nonzero skewness are also as promising as those at zero skewness.Future work may involve the modification of the RGR matching to that of nonzero skewness using the ERBL equation in conjunction with the DGLAP equation.In addition, the results could be further improved by performing the LaMET calculation at multiple boost momenta, P z , in order to make an extrapolation P z → ∞ where the parton model is defined.

FIG. 1 :
FIG. 1: Determination of the linear divergence (left-most plot), δm, by fitting the zero-momentum pion matrix element (blue points) to the function Be −δm z (red curve) in the interval z = [0.54,1.53] fm (shaded green).The error bars for the pion matrix elements are included but too small to be visible.The middle (right-most) plots show the renormalon divergence, m0, determined to (N)NLO (solid orange), (N)NLO+LRR (hatched red), (N)NLO×RGR (hatched green) and ((N)NLO+LRR)×RGR (solid blue).The vertical width of each band corresponds to the systematic error determined from scale variation described in Sec.II A.

2 .FIG. 2 :
FIG. 2: Real (left column) and imaginary (right column) renormalized h RH matrix elements at Q 2 = 0 with the top row showing data points of NLO (blue), NLO+LRR (red), NLO×RGR (green) and (NLO+LRR)×RGR (purple) improvements and bottom row with NNLO (blue), NNLO+LRR (red), NNLO×RGR (green) and (NNLO+LRR)×RGR (purple) improvements.The solid error bars are statistical and the dashed error bars are combined statistical and systematic, the latter arising from the scale variation.Except for NLO and NNLO, the data points shown in the plots have been offset from their exact z value to allow for readability.

FIG. 4 :FIG. 5 :
FIG.4: Real (left column) and imaginary (right column) renormalized h R H (top row) and h R E (bottom row) matrix elements of NLO (blue), NNLO (red), (NLO+LRR)×RGR (green) and (NNLO+LRR)×RGR (purple) improvements at Q 2 = 0.39 GeV 2 .The solid error bars are statistical and the dashed error bars are combined statistical and systematic, the latter arising from the scale variation.Except for NLO (real and imaginary for both h R H and h R E ), the data points shown in the plots have been offset from their exact z value to allow for readability.

FIG. 8 :
FIG.8: Lightcone H (left) and E (right) GPDs evaluated at ξ = 0.1 at NLO (blue) and NLO+LRR (green).The inner error bands are statistical and the outer error bands are combined statistical and systematic errors from scale variations.The vertical dashed lines correspond to x = ±ξ.The GPDs suffer a discontinuity at these x values due to the singularity in the matching kernel.
97 GeV 2 97 GeV 2 97 GeV 2 FIG.6:(NLO+LRR)×RGR (left column) and (NNLO+LRR)×RGR (right column) H (top row) and E (bottom row) GPDs at ξ = 0 and variable Q 2 .The Q 2 ∈ {0.0, 0.19, 0.39, 0.77, 0.97} GeV 2 GPDs are plotted in blue, red, green, purple and orange, respectively.In all cases, the inner error bands are statistical and the outer error bands are combined statistical and systematic errors.The systematic errors decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR, but in both cases are very small.The dark-gray regions are the x-values at which the LaMET calculation breaks down.In addition, when RGR is applied, the matching formula breaks down for |x| ≲ 0.2, which is shaded in light gray.Note that the GPDs are suppressed as Q 2 increases.