Unpolarized and helicity generalized parton distributions of the proton within lattice QCD

We present the first calculation of the $x$-dependence of the proton generalized parton distributions (GPDs) within lattice QCD. Results are obtained for the isovector unpolarized and helicity GPDs. We compute the appropriate matrix elements of fast-moving protons coupled to non-local operators containing a Wilson line. We present results for proton momenta $0.83,\,1.25,\,1.67$ GeV, and momentum transfer squared $0.69,\,1.38$ GeV$^2$. These combinations include cases with zero and nonzero skewness. The calculation is performed using one ensemble of two degenerate mass light, a strange and a charm quark of maximally twisted mass fermions with a clover term. The lattice results are matched to the light-cone GPDs using one-loop perturbation theory within the framework of large momentum effective theory. The final GPDs are given in the $\overline{\rm MS}$ scheme at a scale of 2 GeV.

We present the first calculation of the x-dependence of the proton generalized parton distributions (GPDs) within lattice QCD. Results are obtained for the isovector unpolarized and helicity GPDs. We compute the appropriate matrix elements of fast-moving protons coupled to non-local operators containing a Wilson line. We present results for proton momenta 0.83, 1.25, 1.67 GeV, and momentum transfer squared 0.69, 1.38 GeV 2 . These combinations include cases with zero and nonzero skewness. The calculation is performed using one ensemble of two degenerate mass light, a strange and a charm quark of maximally twisted mass fermions with a clover term. The lattice results are matched to the light-cone GPDs using one-loop perturbation theory within the framework of large momentum effective theory. The final GPDs are given in the MS scheme at a scale of 2 GeV.
Introduction. Quantum Chromodynamics (QCD) is the fundamental theory describing the strong interactions among quarks and gluons (partons). The strong force is responsible for binding partons into hadrons, such as the proton, that makes the bulk of the visible matter in the universe. Studying how the properties of protons emerge from the underlying constituents and their interactions has been an important experimental and theoretical endeavor since the mid-20 th century. These studies led to the realization that high-energy scattering processes can be factorized into perturbative and non-perturbative parts. The latter includes information about the parton structure of the proton [1]. This resulted in the introduction of a complete set of key quantities, namely the parton distribution functions (PDFs) [1], generalized parton distributions (GPDs) [2][3][4], and transverse momentum dependent distributions (TMDs) [5,6]. These describe the non-perturbative dynamics of the proton, and in general hadrons, in terms of their constituent quarks and gluons [7].
The inner structure of hadrons is studied at major experimental facilities using a variety of high-energy processes, and serves as input for collider experiments such as the LHC at CERN. Most of the available information concerns form factors (FFs), measured using elastic scattering, and PDFs, using inclusive and semi-inclusive deep inelastic scattering, as well as Drell-Yan processes. GPDs are accessed through deeply virtual Compton scattering (DVCS) and deeply virtual meson production (DVMP) [8]. To date, GPDs and TMDs are basically unknown.
There are two unpolarized GPDs, H q (x, ξ, t) and E q (x, ξ, t), and two helicity GPDs, H q (x, ξ, t) and E q (x, ξ, t). The superscript q refers to a given quark flavor, and here we study the isovector combination u − d.
GPDs are functions not only of the longitudinal momentum fraction x (0 ≤ x ≤ 1) carried by the partons, but also of the skewness ξ ≡ −∆ + /2P + and the momentum transfer squared, t ≡ ∆ 2 . ∆ + and P + are the plus component of the momentum transfer and the average proton momentum, respectively. Two kinematical regions arise based on the values of ξ and x: the so-called Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) region [9][10][11][12] defined for x > |ξ|, and the Efremov-Radyushkin-Brodsky-Lepage (ERBL) [13,14] region for x < |ξ|. Physical content can be attributed to each region [15] using light-cone coordinates and the light-cone gauge. In the positive-(negative-) x DGLAP region, the GPDs correspond to the amplitude of removing a quark (antiquark) of momentum k from the hadron, and then inserting it back with momentum k + ∆. In the ERBL region, the GPD is the amplitude for removing a quark-antiquark pair with momentum −∆.
While GPDs are multidimensional objects, they lead to simpler quantities when certain limits are taken, or when integrating over selected variables. For example, the forward limit of the unpolarized case, ∆ = 0, gives the quark, f 1 (x) = H q (x, 0, 0), and antiquark PDFs, f 1 (x) = −H q (−x, 0, 0). Equivalently, in the helicity case one has g 1 (x) = H q (x, 0, 0) and g 1 (x) = H q (−x, 0, 0). Integrating over x for nonzero ∆, GPDs give the usual FFs. In the case of the unpolarized GPDs, for example, one obtains the Dirac (F 1 ) and Pauli (F 2 ) FFs, Similar equations can be written for the polarized GPDs. Taking integrals of GPDs over x leads to a tower of Mellin moments that also have a physical interpretation. The integrals can be decomposed into three generalized FFs (A q (t), B q (t), C q (t)), that in the forward limit can be related to the total angular momentum of quarks using Ji's sum rule [2], The connection of GPDs with other quantities demonstrates the information they encode, in both coordinate and momentum spaces. Despite their importance, it is very difficult to extract them experimentally, even though data are available since the early 2000's. These data are limited, covering a small kinematic region, and are indirectly related to GPDs through the Compton FFs. This poses limitations in their extraction, and the fact that more than one independent measurements are needed to disentangle them [16][17][18][19].
Nevertheless, the interest in GPDs is renewed due to the advances both on the experimental and the theoretical side, as well as the expertise gained from recent studies of PDFs. It is, thus, of utmost importance to have ab initio computations of GPDs, that will help map them over different regions of x, ξ, and t. Lattice QCD is the only known formulation that allows a quantitative study of QCD directly using its Lagrangian. Lattice QCD is based on a discretization of Euclidean spacetime and relies on large-scale simulations.
Since parton distributions are light-cone correlation functions [20], it is not straightforward to calculate them using the Euclidean lattice formulation of QCD. The large momentum effective theory (LaMET) proposed by Ji [21] provides a promising theoretical framework to extract light-cone quantities using matrix elements computed in lattice QCD. Within LaMET, one can access light-cone quantities via matrix elements of boosted hadrons coupled with non-local spatial operators, which are calculable on the lattice, and yield what is referred to as quasi-distributions. For large enough momenta, quasidistributions can be related to the physical distributions via a matching procedure [22,23]. The first investigations led to encouraging results on the determination of PDFs [24,25]. Since then, the method has been advanced and attracted a lot of attention, see e.g. Refs. , and revitalized other approaches [64][65][66][67][68][69][70], as well as gave rise to the development and investigation of new ones  (for recent reviews, see Refs. [92,93]). Calculations of the unpolarized, helicity, and transversity PDFs using simulations of twisted mass fermions with physical values of the pion mass were carried out in Refs. [38,42,47]. Recently, a preliminary study of nucleon GPDs was also presented, demonstrating the applicability of the quasidistribution methodology to GPDs [94]. The quasi-GPDs approach has also been studied using the scalar diquark spectator model [95,96].
Extracting GPDs using lattice QCD. For the calculation of GPDs, we define quasi-distributions with boosted proton states and introduce momentum transfer (denoted Q in Euclidean spacetime) between the initial and final states. The matrix element of interest is given by where |N (P i ) (|N (P f ) ) is the initial (final) state labeled by its momentum, and t = −Q 2 . For simplicity, we drop the index q, since in this work we only consider isovector quantities. It is essential that the proton is boosted with large momentum in a given direction, so that the matching of the quasi-GPDs to the light-cone GPDs is valid. The boost is in the direction of the Wilson line (W (0, z)), P 3 = (P i3 + P f 3 )/2. Quasi-GPDs depend on the quasiskewness, defined as ξ = − 2P3 and equal to the light-cone skewness up to power corrections. The Dirac structure Γ defines the type of GPD, and we employ γ 0 and γ 5 γ 3 for the unpolarized and helicity GPDs, respectively [122].
Another aspect of the calculation is the renormalization, as the divergences with respect to the regulator must be removed prior to applying Eq. (6). We adopt [123] the non-perturbative renormalization scheme of Refs. [29,30], and refined in Ref. [47]. This procedure removes all divergences, including the powerlaw divergence with respect to the ultraviolet cutoff. The renormalization functions, Z Γ , are obtained nonperturbatively by imposing RI-type [97] renormalization conditions, given in Eq. (S7). In a nutshell, the final values of Z Γ are obtained at each value of z separately, at a chosen RI scale (aµ 0 ) 2 . For each value of z at a given µ 0 , we take the chiral limit using a linear fit in m 2 π . As described in the supplement, the available matching equations [98] require that the quasi-GPDs are in the RI scheme. Therefore, we renormalize the matrix elements using the estimates for Z Γ in the RI scheme at a given scale, µ 0 , chosen to be (aµ 0 ) 2 ≈ 1.17. This scale enters the matching kernel, which converts the quasi-GPDs to light-cone GPDs. The latter are always given in the MS scheme at 2 GeV, regardless of the scheme used for quasi-GPDs. Within this work, we explored a few values of the scale within the range (aµ 0 ) 2 ∈ [1 − 5]. We find that the dependence on (aµ 0 ) 2 is within the reported uncertainties.
The renormalized matrix elements are decomposed into the form factors {F H , F E } and {F H , F E }, for the unpolarized and helicity case, respectively. The decomposition is based on continuum parametrizations, which in Euclidean space take the form where Q ≡ P f − P i , and m is the proton mass. O Γ (z) ≡ ψ (z) ΓW (0, z)ψ (0) and Γ ≡ū N (P f , s ) Γ u N (P i , s) with u N the proton spinors.
The matrix elements h Γ (z, P 3 , t, ξ) depend on z, which varies from zero up to the half of the spatial extent L of the lattice. One way to reconstruct the x-dependence of the GPDs is via a standard Fourier transform, e.g., we define the quasi H-GPD as H q : This simple Fourier transform suffers from an ill-defined inverse problem [81]. One alternative reconstruction technique that we adopt here is the Backus-Gilbert (BG) method [99] that leads to a uniquely reconstructed quasidistribution from the available set of matrix elements. More details can be found in the supplement.
The matching formula is available to one-loop level in perturbation theory, for general skewness [98] [124]. In fact, in the limit of ξ → 0, one recovers the matching equations for quasi-PDFs. Furthermore, the matching kernels of H-and E-GPDs are the same [98]. We provide details on the matching in the supplement.
Numerical techniques. For this calculation, we employ an ensemble with two light, a strange and a charm quark (N f = 2 + 1 + 1) using the twisted mass formulation [100,101] with clover improvement [102], generated by the Extended Twisted Mass Collaboration (ETMC) [103]. The ensemble has a spatial (temporal) extent of 3 fm (6 fm) (32 3 ×64), a lattice spacing of 0.093 fm and pion mass of about 260 MeV. For the isovector combination u−d, we need to evaluate only the connected diagram (see Fig. S1).
To increase the signal-to-noise ratio, we use momentum smearing [104], which has been very successful in the calculation of matrix elements of non-local operators with boosted hadrons [27,38,47,103]. We find that momentum smearing decreases the gauge noise of the real (imaginary) part by a factor of 4-5 (2-3) (see, e.g., Fig. S2). To further suppress statistical uncertainties, we apply stout smearing [105] to the links of the operator. The effectiveness of the stout smearing in proton matrix elements was demonstrated in Refs. [106,107].
Ensuring ground-state dominance in h Γ is essential and is controlled by the time separation between the source (initial state) and the sink (final state). This separation, t s , needs to be large in order to suppress excitedstates contributions to the matrix elements. We construct a suitable ratio of two-and three-point functions (see Eq. (S4)), to cancel out unknown overlap factors. Multiple ratios are obtained, for each operator insertion time t ins = 1, . . . , t s − a (assuming the source time is zero). Ground-state dominance is established when the ratio becomes time independent for values of t ins (plateau region) that are far away enough from the source and the sink (see Eq. (S5)). The matrix elements h Γ (z, P 3 , t) are extracted from a constant fit within the plateau region. Here, we choose t s = 1.12 fm [47], and use the sequential method at fixed t s value.
The most common definition of GPDs is in the Breit frame, in which the momentum transfer Q is equally shared between the initial and final states. This has important implications for the computational cost of extracting h Γ (z, P 3 , t, ξ) as compared to the usual FFs. For different momentum transfers, both the source and the sink momenta change, requiring separate inversions for each value of Q. The statistics used for the results presented in this work is given in Tabs. SI-SII. We note that, for the largest value of proton momentum, P 3 = 1.67 GeV, the number of measurements required to reach sufficient accuracy is 112 192.
Results for the matrix elements h Γ . The renormalized matrix elements are decomposed into F H , F E , F H , and F E using Eqs. (4)- (5). To disentangle F H and F E , we use h γ 0 (z, P 3 , t, ξ) projected with the unpolarized pro-jector, P 0 ≡ (1 + γ 0 )/4 and the polarized projector, P κ ≡ (1 + γ 0 )iγ 5 γ κ /4. For the helicity matrix element, h γ 3 γ 5 (z, P 3 , t, ξ), we use the polarized projector, P κ , where both κ = 3 and κ = 3 are necessary to disentangle F H and F E . We note that for zero skewness, only κ = 3 leads to a non-zero matrix element for the axial vector operator, which is related to F H . Thus, for ξ = 0 we cannot access the E-GPD. In fact, the inaccessibility of E-GPD is a general feature due to its vanishing kinematic factor at ξ = 0, and is not related to the choice of the projector.
For the largest momentum, P 3 = 1.67 GeV, we find similar magnitude contributions from both projectors P 0 and P 1 . These matrix elements are combined to solve a system of linear equations to extract F H and F E . Due to its kinematic coefficient, F E has, in general, larger errors than those for F H . We find that the momentum dependence changes based on the values of z, and on the quantity under study. This momentum dependence propagates in a nontrivial way to the final H-and E-GPDs, as one has to reconstruct the quasi-GPDs in momentum space, and then, apply the appropriate matching formula, which depends on the momentum P 3 . The matrix element h γ 5 γ 3 at zero skewness leads directly to F H , as the kinematic factor of E is zero. More details and plots can be found in the supplement.
Results on the GPDs. The P 3 -convergence of the GPDs is of particular interest, as the matching kernel is only known to one-loop level. For H-GPD and H-GPD at ξ = 0, we find that the momentum dependence is small and within the reported uncertainties. Convergence is also observed for E-GPD for the two highest momenta and the region x > 0. We note that the statistical errors on E-GPD are larger than those of the H-GPD, a feature already observed in F E . We refer the Reader to the supplement for more details. and zero skewness are shown in Fig. 1 and Fig. 2 for the unpolarized and helicity GPDs, respectively. For each case, we compare the GPDs with the corresponding PDFs, that is f 1 (x) for the unpolarized, and g 1 (x) for the helicity. We observe that the GPDs are suppressed in magnitude as compared to their respective PDFs for all values of x 0.7. In fact, H-GPD has a steeper slope at small x values. The smaller magnitude of the GPDs is a feature also observed in the standard FFs, which decay with increasing −t. For the large-x region, both distributions decay to zero in the same way. For the antiquark region, we find that the GPDs are compatible with the corresponding PDFs. We note that the statistical uncertainties of GPDs are similar to the PDFs, allowing for such qualitative comparison.
The extraction of the GPDs for ξ = 0 differs from the one for ξ = 0, as a different matching kernel is required. Also, unlike the ξ = 0 case, both helicity GPDs contribute to the matrix element, and therefore a decomposition is required. The comparison between the zero and non-zero skewness is shown in Fig. 3 and Fig. 4, for P 3 = 1.25 GeV. The main feature of the GPDs at ξ = 0 is that an ERBL region (|x| < 1/3 in our case) appears, differentiating it from the DGLAP region (|x| > 1/3). The behavior of the GPDs as a function of t for a fixed x is as expected; increasing −t suppresses the GPDs.
Concluding remarks. We presented first results on the unpolarized and helicity GPDs for the proton, employing the quasi-distribution approach, which has been very successful for the extraction of PDFs within lattice QCD. In the case of GPDs, a non-zero momentum is transferred between boosted initial and final states. The lattice QCD data were renormalized non-perturbatively, and the Backus-Gilbert method was used to extract the x-dependence of quasi-GPDs. Applying matching to the latter within the LaMET approach yielded the light-cone GPDs in the MS scheme at 2 GeV.
The momentum dependence of GPDs for P 3 =  at nonzero skewness, are reassuring, as with increasing −t, the magnitude of GPDs is suppressed. With our calculation, we demonstrate that extracting GPDs with controlled statistical uncertainties is feasible. Their accuracy permits qualitative comparison with their corresponding PDFs.
In the near future, we will investigate systematic uncertainties, as studied for PDFs [47]. The pion mass dependence will also be investigated using an ensemble with quark masses fixed to their physical values. In a follow-up calculation, we will also explore the chiral-odd transversity GPD, for which there are two additional form factors, leading to a more evolved decomposition. This makes the disentanglement of the transversity GPDs more challeng-ing.
The current work demonstrates the feasibility of the quasi-distributions approach for GPDs using computational resources that are within reach. Extracting GPDs from first principles can potentially be combined with future experimental data, as GPDs are at the heart of planned experiments at JLab [108] and the Electron-Ion-Collider (EIC) [109]. Therefore, GPDs are the objects to drive the efforts of the nuclear and hadronic physics communities for the next decades. The work is based on the calculation of proton matrix elements of the nonlocal operator containing a Wilson line in the z direction, W (0, z), that is (S1) An important requirement of the quasi-distribution approach is that the hadron is boosted with a momentum in the same direction as the Wilson line, therefore P = (0, 0, P 3 ). GPDs are multidimensional objects and require momentum transfer between the initial and final states. In Euclidean space, this is defined as Q 2 , which is related to its Minkowski counterpart as t ≡ −Q 2 . An important parameter of GPDs is the skewness, which is proportional to the momentum transfer in the direction of the boost. In the quasi-distribution method, the relevant quantity is the quasi-skewness defined as ξ ≡ − Q3 2P3 .
N ( x, t s ) N ( 0, 0) We calculate the isovector flavor combination, which requires calculation of only the connected diagram shown in Fig. S1. To obtain the ground state of the matrix element, one must calculate two-point and three-point correlation functions, where N α (x)= abc u a α (x) d b T (x)Cγ 5 u c (x) is the interpolating field for the proton, τ is the current insertion time, and Γ αβ and Γ αβ are parity projectors. Without loss of generality, we take the source to be at (0, 0). Γ αβ is the parity plus projector Γ= 1+γ4 2 , and Γ α β = 1 4 (1 + γ 0 )iγ 5 γ j . For nonzero momentum transfer, one must form an optimized ratio in order to cancel the time dependence in the exponentials and the overlaps between the interpolating field and the nucleon states, In the limit (t s − τ ) a and τ a, the ratio of Eq.(S4) becomes time-independent and the ground state matrix element is extracted from a constant fit in the plateau region, that is To improve the overlap with the proton ground state, we construct the proton interpolating field using momentumsmeared quark fields [104], on APE smeared gauge links [110]. The momentum smearing technique was proven to be crucial to suppress statistical uncertainties for matrix elements with boosted hadrons, and in particular for nonlocal operators [27]. In this work, we can reach P 3 = 1.67 GeV at a reasonable computational cost. The momentum smearing function S on a quark field, ψ, reads where α G is the parameter of the Gaussian smearing [111,112], U j is the gauge link in the spatial j-direction. P = (P 1 , P 2 , P 3 ) is the momentum of the proton (either at the source, or at the sink) and ξ is a free parameter that can be tuned so that a maximal overlap with the proton boosted state is achieved. For ξ = 0, Eq. (S6) reduces to the Gaussian smearing function. In our implementation, we keep ξP parallel to the proton momentum at the source and at the sink. Such a constraint requires separate quark propagators for every momentum transfer, because the gauge links are modified every time by a different complex phase. However, this strategy avoids potential problems due to rotational symmetry breaking. It also has the benefit that every correlator entering the ratio of Eq. (S4) is optimized separately. As an example, we show in Fig. S2 the bare matrix elements of the vector operator with and without momentum smearing. For this comparison, we use the unpolarized parity projector, momentum boost P 3 = 0.83 GeV, and momentum transfer t = −1.39 GeV 2 . The individual components of the momentum transfer are given by the vector Q L 2π = (0, −2, 2). The number of measurements is 1616 for t s = 8a, and the data using the momentum smearing correspond to ξ = 0.6 after optimization. As can be seen, the use of momentum smearing significantly decreases the statistical uncertainties for both the real and imaginary parts of the matrix elements. In particular, the statistical accuracy increases by a factor of 4-5 in the real part, and 2-3 in the imaginary part, depending on the value of z.
In Table SI, we summarize the statistics for each value of the nucleon momentum P 3 , the momentum transfer squared t, and skewness ξ. In Table SII, we give the numbers of measurements for the corresponding PDFs, to which the computed GPDs can be compared.

Renormalization
We employ an RI-type renormalization prescription, using the momentum source method [113,114] that offers high statistical accuracy. The appropriate conditions for the renormalization functions of the nonlocal operator, Z Γ , and the quark field, Z q , are Tr (S(p)) −1 S Born (p) Note that Eq. (S7) is applied at each value of z separately. V(p, z) (S(p)) is the amputated vertex function of the operator (fermion propagator) and S Born (p) is the tree-level of the propagator.
We use five N f = 4 ensembles as given in Tab. SIII, which have been produced for the calculation of the renormalization functions at the same β value as the N f = 2 + 1 + 1 ensemble used for the extraction of the matrix elements of Eq. (S1). The RI renormalization scale µ 0 , defined in Eq. (S7), is chosen appropriately to have suppressed discretization effects, as explained in Ref. [114]. We employ several values that have the same spatial components, that is p = (p 0 , p 1 , p 1 , p 1 ), so that the ratio p 4 (p 2 ) 2 is less than 0.35, as suggested in Ref. [115]. In this work, we use different values of µ 0 ((a µ 0 ) 2 ∈ [1, 5]) to check the dependence of the matching formalism on µ 0 . For each µ 0 value, we apply a chiral extrapolation using the fit to extract the mass-independent Z RI Γ,0 (z, µ 0 ). As mentioned in the main text, the matching kernel of Ref. [98] requires that the quasi-GPDs are renormalized in the RI scheme. For consistency, we use the same for ξ = 0 and use the renormalization functions defined on a single RI renormalization scale, (aµ 0 ) 2 ≈ 1.17. This scale also enters the matching equations. We find negligible dependence when varying µ 0 .

Reconstruction of x-dependence
The renormalized matrix elements F G , where G = H, E, are related to the quasi-distributions G q by a Fourier transform: (S10) Inverting this expression relates the quasi-GPDs to the matrix elements. However, the inverse equation involves a Fourier transform over a continuum of lengths of the Wilson line, up to infinity, while the lattice provides only a discrete set of determinations of F G , for integer values of z/a up to roughly half of the lattice extent in the boost direction, L/2a. Thus, the inversion of Eq. (S10) poses a mathematically ill-defined problem, as argued and discussed in detail in Ref. [81]. The inverse problem originates from incomplete information, i.e. attempting to reconstruct a continuous distribution from a finite number of input data points. As such, its solution necessarily requires making additional assumptions that provide the missing information. These assumptions should be as mild as possible and preferably model-independent -else, the reconstructed distribution may be biased.
One of the approaches proposed in Ref. [81] is to use the Backus-Gilbert (BG) method [99]. The model-independent criterion used in the BG procedure, to choose from among the infinitely many possible solutions to the inverse problem, is that the variance of the solution with respect to the statistical variation of the input data should be minimal. The reconstruction proceeds separately for each value of the momentum fraction x. In practice, we separate the exponential of the Fourier transform into its cosine and sine parts, related to the real and imaginary parts of the matrix elements, respectively. We define a vector a K (x), where K denotes either the cosine or sine kernel, of dimension d equal to the number of available input matrix elements, i.e. d = z max /a + 1, where z max /a is the maximum length of the Wilson line (in lattice units) used to determine the quasi-distribution. The BG procedure consists in finding the vectors a K (x) for both kernels according to the variance minimization criterion. The vector a K (x) is an approximate inverse of the cosine/sine kernel function K(x), that is: and K(x ) z/a = cos(x P 3 z) or K(x ) z/a = sin(x P 3 z) are elements of a d-dimensional vector of discrete kernel values corresponding to available integer values of z/a. The function ∆(x − x ) is, thus, an approximation to the Dirac delta function δ(x − x ), with the quality of this approximation depending, in practice, on the achievable dimension d at given simulation parameters.
The vectors a K (x) are found from optimization conditions resulting from the BG criterion. We refer to Ref. [81] for their explicit form and here, we just summarize the final result. We define a d × d-dimensional matrix M K (x), with matrix elements where x c is the maximum value of x for which the quasi-distribution is taken to be non-zero (i.e. its reconstruction proceeds for x ∈ [0, x c ]) and the parameter ρ regularizes the matrix M K . This regularization, proposed by Tikhonov [116], was put up as one possible way of making M K invertible [81,117,118]). The value of ρ determines the resolution of the method and should be taken as rather small, in order to avoid a bias. We use ρ = 10 −3 , which leads to reasonable resolution and is large enough to avoid oscillations in the final distributions related to the presence of small eigenvalues of M K . Additionally, we define a d-dimensional vector u K , with elements The above mentioned optimization conditions lead to: and the final BG-reconstructed quasi-distributions are given by (S15)

Matching Procedure
Contact between the physical distributions and the quasi-GPDs is established through a perturbative matching procedure. The factorization formula for the Dirac structure Γ takes the form

Results
In this section, we provide more details for the extracted matrix elements and the final GPDs.
In Fig. S3, we show the bare matrix elements for the vector operator, using the projectors P 0 (h γ0,P0 ) and P 1 with κ = 1 (h γ0,P1 ). Note that for the polarized projector, only κ = 1 contributes, as the momentum transfer has zero component in that direction. Focusing on the largest value of the momentum P 3 = 1.67 GeV, one observes that both h γ0,P0 and h γ0,P1 give similar contributions. The decomposition of the renormalized matrix elements leads to F H and F E , shown in Fig. S4. It is interesting to observe that the statistical errors for F E are, in general, larger than those for F H . This effect has its origin in the kinematic coefficients of F E in the decomposition of the matrix element. We find that the momentum dependence changes based on the values of z, and on the quantity under study. This momentum dependence propagates in a nontrivial way to the final H-and E-GPDs, as one has to reconstruct the quasi-GPDs in momentum space, and then, apply the appropriate matching formula, which depends on the momentum P 3 . In Fig. S4 we show the decomposed quantities F H and F E , which have been disentangled using the renormalized matrix elements h γ 0 ,P0 and h γ 0 ,P1 . We find negligible momentum dependence in Re[F H ] for 0 < z < 5, while it has a steeper and compatible slope for the two largest momenta in the region 5 ≤ z ≤ 9. Re[F H ] flattens out for all three momenta for z ≥ 10, which are consistent. Similarly, Im[F H ] is compatible for the largest two momenta for all values of z, while the lowest momentum evinces a clearly slower decay to zero. This might indicate onset of convergence above P 3 = 1.25 GeV. Note, however, that the matching depends on P 3 and thus, more conclusions about convergence can be drawn after applying this procedure. For Re[F E ], the errors are significantly larger than for Re[F H ], as remarked above, and we observe somewhat slower decay at P 3 = 1.25 GeV as compared to P 3 = 1.67 GeV. This may indicate slower convergence in the E-GPD, but may also be a statistical effect, since Im The matrix element h γ 3 γ 5 is shown in the left panel of Fig. S5 for t = −0.69 GeV 2 and ξ = 0. The corresponding F H is shown in the right panel. We note that for zero skewness, the kinematic factor of E is zero, and we only extract H from the lattice QCD data. We observe that both for the real and the imaginary part of F H , there are significant differences between the largest two momenta. Thus, we postpone conclusions about convergence to the discussion of the final GPDs.  We now move on to the discussion of the final GPDs, in particular their convergence in momentum for zero skewness (t = −0.69 GeV 2 ). We compare the unpolarized GPDs for P 3 = 0.83, 1.25, 1.67 GeV in Fig. S8. The H-GPD has negligible P 3 -dependence for every region of x, while the E-GPD exibits convergence between the two largest momenta for x > 0, which is of main interest. We note that the statistical errors on E-GPD are larger than those of the H-GPD, a feature already observed in F E (see Fig. S3). In Fig. S9, we show the momentum dependence of the H-GPD. We observe that the relatively large differences between the renormalized F H for the lowest two momenta and P 3 = 1.67 GeV are compensated by the matching procedure, indicating final convergence within the reported statistical uncertainties. Thus, this conclusion holds for both the unpolarized and the helicity H-GPD.  In Fig. S10 we provide a comparison between the H-and E-GPDs (left panel), and H-and E-GPDs (right panel). In the unpolarized case, H-and E-GPDs are compatible with each other in the quark region. However, the Pauli FF, corresponding to the x-integral of the E-GPD, is considerably larger than the Dirac FF (integral of H-GPD) at this momentum transfer and this is achieved by the larger values of the E-GPD in the antiquark region. For the helicity case, the E-GPD is significantly larger than the H-GPD, which reflects the fact that the axial-vector FF G P is found to be a factor ≈ 3 larger than the G A at this momentum transfer, in a lattice setup with similar parameters [121]. We also note that the integrals of H-,E-, H-and E-GPDs extracted in this work are all compatible with their respective FFs obtained in Ref. [121].