Generalized Parton Distributions from Lattice QCD with Asymmetric Momentum Transfer: Unpolarized Quarks

Traditionally, lattice QCD computations of generalized parton distributions (GPDs) have been carried out in a symmetric frame, where the transferred momentum is symmetrically distributed between the incoming and outgoing hadrons. However, such frames are inconvenient since they require a separate calculation for each value of the momentum transfer, increasing significantly the computational cost. In this work, by focusing on the quasi-distribution approach, we lay the foundation for faster and more effective lattice QCD calculations of GPDs exploiting asymmetric frames, with freedom in the transferred momentum distribution. An important ingredient of our approach is the Lorentz covariant parameterization of the matrix elements in terms of Lorentz-invariant amplitudes, which allows one to relate matrix elements in different frames. We also use this amplitude approach to propose a new definition of quasi-GPDs that is frame-independent and, more importantly, may lead to smaller power corrections in the matching relations to the light-cone GPDs. We demonstrate the efficacy of the formalism through numerical calculations using one ensemble of $N_f$=2+1+1 twisted mass fermions with a clover improvement. The value of the light-quark masses lead to a pion mass of about 260 MeV. Concentrating on the proton, and limiting ourselves to a vanishing longitudinal momentum transfer to the target, we extract the invariant amplitudes from matrix element calculations in both the symmetric and asymmetric frame, and obtain results for the twist-2 light-cone GPDs for unpolarized quarks, that is, $H$ and $E$.


I. INTRODUCTION
Parton distribution functions (PDFs), which are measurable in processes like inclusive deep-inelastic lepton-nucleon scattering, are key objects containing information about the quark and gluon structure of strongly interacting systems [1]. They provide 1D images of hadrons by describing how the partons are distributed as a function of the momentum fraction x they carry of the hadron's momentum. PDFs are defined through matrix elements of bi-local quark or gluon operators, with the parton fields having a light-like separation and the operators evaluated for the same initial and final hadron state. Generalized parton distributions (GPDs) are generalizations of the concept of PDFs in that the light-like parton operators are computed for different initial and final states [2][3][4]. GPDs therefore depend on the longitudinal momentum transfer ξ and the invariant momentum transfer t to the target, in addition to their dependence on the parton momentum fraction x. While this makes them complicated multi-variable functions, the information encoded in GPDs is much richer than for PDFs. In particular, they provide 3D images of hadrons [5][6][7][8], give access to the angular momenta of partons [3], and have a relation to pressure and shear forces inside hadrons [9][10][11]. The physics of GPDs has been discussed in various review articles [12][13][14][15][16][17][18][19].
Experimental information on GPDs can be obtained from hard exclusive scattering processes such as deep virtual Compton scattering [2-4, 20, 21] and hard exclusive meson production [22][23][24]. But extracting GPDs from such reactions in a model-independent manner is very complicated, mainly because in the observable quantities, like the Compton form factors, the momentum fraction x is integrated over -see [25] for a recent discussion and detailed analysis of this issue. It is therefore very desirable to obtain information on GPDs from first principles in lattice QCD. However, lattice QCD calculations of light-cone correlation functions like PDFs and GPDs are challenging due to the time dependence of those objects. As a result, for a long time lattice QCD calculations were limited to the lowest Mellin moments of the GPDs [26][27][28][29][30], with simulations at the physical point available only in recent years [31][32][33][34][35][36][37][38][39][40][41][42][43][44]. Despite the progress, their dependence on x remained elusive.
Applications of these new developments in the case of GPDs are still somewhat sparse in comparison. Nevertheless, we have seen results for matching [119][120][121][122][123], model studies [124][125][126][127][128] and, in particular, the first pioneering lattice QCD calculations for the pion [129] and the nucleon [130][131][132][133][134][135]. These results are very encouraging, demonstrating explicitly that GPDs can be obtained on the lattice. But it is rather clear that the full mapping of GPDs with respect to their variables, in particular the momentum transfer t and the skewness ξ, is very challenging and computationally much more demanding than for PDFs. Among the reasons for this is that, so far, they have been computed in symmetric frames of reference, where the momentum transfer is equally split between the source and the sink. Consequently, every value of the momentum transfer is obtained from a separate and costly calculation. Here, for the first time, we consider asymmetric frames for the computation of GPDs. As we will demonstrate below, this allows for more efficient calculations, since different momentum transfers can be obtained in a single calculation.
The paper is organized as follows. In Section II, we discuss the kinematics of the symmetric and asymmetric frames for our study and how the two frames can be related through a Lorentz transformation. For both a spin-0 and spin-1 2 target, we introduce and discuss the main theoretical tool of Lorentz-invariant amplitudes in terms of which matrix elements that define GPDs can be parameterized. Based on those amplitudes we also propose a new, frame-independent definition of quasi-GPDs which, in comparison to previously used quasi-GPDs, may converge faster to their respective light-cone GPDs. Section III specializes on the Euclidean case and provides decompositions of lattice-calculable matrix elements in terms of these amplitudes, and our lattice setup. Numerical results are shown in Section IV, where we make a detailed comparison of the symmetric and the asymmetric frames at various stages, in coordinate space and in momentum space. For the proton and ξ=0 we show, in particular, numerical results for the invariant amplitudes and the twist-2 light-cone GPDs H and E. Section V concludes and discusses future prospects.

A. Symmetric and asymmetric frames
The initial and final momentum of the hadron are frame-dependent, and the most widely used frame of reference to calculate GPDs is the symmetric frame. In such a case, the momentum transfer is symmetrically distributed between the incoming (p i ) and the outgoing hadrons (p f ) (see, e.g., Eqs. (48) - (49)). In this section, it is convenient to define the momenta in terms of the average momentum P = 1 2 (p i + p f ) and the momentum transfer ∆ = p f − p i , (1) The above expressions are general for any frame, but will differ numerically in each frame. An alternative setup to the symmetric one, is an asymmetric frame, where the momentum transfer is not shared between the incoming and outgoing hadrons but is rather entirely applied to the incoming hadron (see, e.g., Eqs. (50) - (51)). This frame is of interest for this work and will be used throughout this paper. For completeness, we remind the reader that the energies of the initial and final states, E i/f = m 2 + ( p i/f ) 2 (where m is the mass of the hadron), are also different in the two frames.
While computations in the symmetric frame have been extensively used in model calculations [124,125], due to the nice symmetry properties of the correlators, they are notoriously difficult to calculate in lattice QCD mainly due to the computational cost to extract a range of values for the momentum transfer. More specifically, the information on the momentum transfer is present in both the initial and final states. As a consequence, every value of ∆ requires a separate calculation. Such a constraint places severe limitations on GPD calculations in terms of the range of values of the momentum transfer that can be accessed, and, consequently, skewness. So the question arises whether it is meaningful to calculate GPDs in asymmetric frames, which can be computationally less expensive. As we will see below, one of the approaches to handle calculations in asymmetric frames is to relate the setup of the symmetric frame to the asymmetric one via an appropriate Lorentz transformation. For instance, a Lorentz transformation along the z-direction is not optimal for this purpose because it requires a spatial operator distance (say z = (0, 0 ⊥ , z 3 = 0)) to develop a non-zero temporal component (that is z → (z 0 = 0, 0 ⊥ , z 3 )) which is problematic for lattice QCD calculations. However, any Lorentz transformation transverse to the z-direction leaves the operator distances unchanged. This means that if one begins with a spatial operator distance in one frame, its counterpart in the other frame of reference remains spatial. Such a transformation is called "transverse boost". We illustrate this point below by focusing on the simplest case of a transverse boost in the x-direction and zero skewness. Note that this method can be generalized for any general transverse boost and for an arbitrary value of skewness. Let us begin with relating the incoming state in the two frames, p s i = (E s i , −∆ 1,s /2, 0, P 3 ) and p a i = (E a i , −∆ 1,a , 0, P 3 ). Lorentz transformation provides p s = Λ LT p a , This gives, and, Now, we relate the outgoing state in the two frames, p s f = (E s f , ∆ 1,s /2, 0, P 3 ) and p a f = (E a f , 0, 0, P 3 ). (Note that the energies of the incoming and outgoing states are different in the asymmetric frame.) Following the steps outlined above, we find, and, element ..γ 0 .. in the symmetric frame can be written as a linear combination of matrix elements of different operators ..(γ 0 + γ 1 ).. in the asymmetric frame, This equation essentially shows how the 0 th component of a 4-vector changes under a Lorentz transformation. This implies that a transverse boost that uniquely fixes (β, γ) (Eqs. (10) and (11)) allows for an exact calculation of quasi-GPDs in the symmetric frame through matrix elements of the asymmetric frame. However, Eq. (15) also makes it clear that a quasi-GPD defined through the operator γ 0 is intrinsically Lorentz non-invariant. In the limit of a large momentum, though, we recover implying that the contribution from the matrix element ..γ 1 .. is essentially a power correction at finite values of momentum P 3 .
We now illustrate the Lorentz non-invariance of the (above) historic definition of quasi-GPD through an altogether different approach and then motivate a new definition for quasi-GPDs that is frame independent, and more importantly, may potentially reduce power corrections. We call this the second approach or the amplitude approach. As a first step, we build a Lorentz-covariant decomposition of the matrix element in Eq. (14) in terms of the available vectors (P µ , z µ , ∆ µ ), where m is the mass of the target. Here, A i 's are the Lorentz-invariant (and, thus, frame-independent) amplitudes whose arguments are functions of Lorentz scalars 1 . The light-cone GPD H in both symmetric and asymmetric frames is defined from the correlator where z µ = (0, z − , 0 ⊥ ). Note that ∆ + /P + = z · ∆/z · P , so the above GPD is the same in both frames as long as z · P , z · ∆ and ∆ 2 are held to be the same. The light-cone GPD H in the momentum space is defined, where the skewness parameter ξ = −∆ + /(2P + ). In the literature, the light-cone GPD has also been defined in the symmetric frame as [3,20] H for a lightlike vector n µ ∝ (1, 0, 0, −1). In the symmetric frame, the average momentum P has its dominant component along the light-cone direction that is anti-collinear to n. The above expression allows us to generalize H as which is independent of the orientation of z µ and equivalent to Eq. (18) in the coordinate system where z µ = (0, z − , 0 ⊥ ). Therefore, H is Lorentz invariant as long as the scalars z · P, z · ∆, ∆ 2 are fixed, and the H GPD in the momentum space is the Fourier transform by integrating along a fixed direction in the (z · P, z · ∆) plane with z · ∆ = −2ξ(z · P ), i.e., (Note that x is the Fourier conjugate of z · P .) Now, we turn to the quasi-GPD H which in the coordinate space is connected to the light-cone GPD H through the matching formula [122], whereC is the short-distance matching coefficient that can be calculated perturbatively [121][122][123], and µ is the renormalization scale in the MS scheme. (We will revisit the derivation of the matching equation towards the end of this section.) At leading order in α s , the above formula indicates that H collapses to H in the light-cone limit z 2 → 0, Therefore, a natural candidate for a frame independent quasi-GPD is the generalization of Lorentz-invariant H to include z 2 = 0, i.e., Note that this result in the forward limit agrees with the quasi-PDF definition using the γ 0 matrix element [57]. Since both sides of Eq. (24) are Lorentz invariant (recall also Eq. (21)), at finite z 2 the difference between H and H are frame-independent subleading power corrections in A i 's. Correspondingly, the quasi-GPD H is defined as where the measure d(z · P ) = −P 3 dz 3 with fixed P 3 .
A direct implication of Eq. (23) is that the skewness of the GPD H, ξ = −∆ + /(2P + ), is equal to the quasi-skewness ξ = −∆ 3 /(2P 3 ) of the corresponding GPD H, as they both are given by −z · ∆/(2z · P ). To better understand this, let us recall the derivation of the factorization formula for the quasi-GPD. At short z 2 , the matrix element in Eq. (14) has an operator product expansion (OPE) [121], where C n are Wilson coefficients, F n is a special polynomial series, O µ0µ1...µn are the conformal operators, and B n,m (µ) are perturbative coefficient functions that diagonalizes the anomalous dimension matrix of the operators that mix with O µ0µ1...µn . The conformal operator is defined as [136,137] n µ0 n µ1 . . . n µn O µ0µ1...µn = where C 3/2 n,m x 2m is a Gegenbauer polynomial in x. The off-forward matrix element of O µ0µ1...µn is the Gegenbauer moment, which, according to Lorentz covariance, can be parameterized as where φ n,m (t) are frame-independent form factors.

C. Spin-1/2 particles
In this section, we turn our attention to spin-1/2 particles, such as the proton. As in the case of spin-0 particles, it is crucial to derive a Lorentz-covariant decomposition of the vector matrix element for spin-1/2 particles. It turns out that constraints from parity alone are sufficient to write down the general structure of the vector matrix element. One ends up in finding eight linearly-independent Dirac structures multiplied by eight Lorentz-invariant (frame-independent) amplitudes, where with a summation implied for repeated indices. Also, we use the compact notation A i ≡ A i (z · P, z · ∆, ∆ 2 , z 2 ). The steps involved in the derivation of Eq. (35) are outlined in Appendix A. This derivation parallels the steps presented in Ref. [138]. (See also Ref. [139] where this matrix element was parameterized in momentum space and for a straight Wilson line.) Note that the amplitudes A 1 , A 2 , and A 3 are analogous to the spin-0 case. We also note that one can choose to work with a basis of different parametrization other than Eq. (35). However, the number of amplitudes will remain the same and, hence, one would always require eight independent lattice matrix elements to disentangle the amplitudes. Therefore, there is no obvious gain in computational cost if a different parametrization is used.
For spin-1/2 particles, there are two (vector) light-cone GPDs H and E defined through [7], After using µ = + in Eq. (35), we can perform a change of basis of the resulting expression to map the A i 's onto the GPDs in Eq. (36), where the arguments of the A i 's have no dependence on z 2 . We can make the above expressions formally Lorentz invariant as (see Sec. II B), We emphasize that one can arrive at Eqs. (39) and (40) by contracting both sides of Eq. (35) with z µ (where z µ is an arbitrary light-like vector) and by ensuring that z 2 = 0 (recall Eq. (20) and Eq. (21)). This implies that light-cone GPDs are frame-independent as long as the Lorentz scalars such as (z · P s/a , z · ∆ s/a , (∆ s/a ) 2 ) are taken to be the same in the two frames.
We now turn to quasi-GPDs. As emphasized in Sec. II B, the essence of the matching equation is the equivalence of the quasi-GPDs and the light-cone GPDs at the leading order. Therefore, a natural way to define the quasi-GPDs H and E is through a Lorentz-invariant generalization of the light-cone definitions in Eqs. (39) and (40) to z 2 = 0, i.e., where now the arguments of the A i 's have a non-zero dependence on z 2 . We expect that the definitions in Eqs. (41)-(42) may have a faster convergence to the light-cone GPDs at the leading order, although such a statement needs a rigorous justification from the theory side. 2 Furthermore, these definitions of quasi-GPDs differ from their (respective) light-cone GPDs by frame-independent power corrections beyond the leading order.
Historically, quasi-GPDs have been defined through the γ 0 operator as, After using µ = 0 in Eq. (35), we can perform a change of basis of the resulting expression to map the A i 's onto the quasi-GPDs in Eq. (43). The relations in the symmetric frame read, On the other hand, the relations in the asymmetric frame read, Several comments are in order: For finite values of the momentum, the above expressions contain additional amplitudes that are not present in the light-cone expressions. Thus, contrary to their forward limit where arguments are made in favor of γ 0 because of reduced amplitudes [57], here we find that additional amplitudes are found in γ 0 in the off-forward limit. (Note that the different definitions of quasi-GPDs preserve the norm; see also Ref. [125].) Second, the intrinsic Lorentz non-invariance associated with the historical definitions of quasi-GPDs (formally) implies that the basis vectors (γ 0 , iσ 0∆ s/a ) do not form a complete basis for a spatially-separated bi-local operator for finite values of momentum. Therefore, for defining quasi-GPDs in a Lorentz-invariant manner, one needs to have a different basis other than just (γ 0 , iσ 0∆ s/a ). From this perspective, one can infer that our Lorentz-invariant definition of quasi-GPDs is essentially a redefinition of quasi-GPDs in terms of a suitable linear combination of operators (which turns out to be γ 1/2 ) that reduces the additional amplitudes present in the historic definitions. This may potentially provide a faster convergence to the light-cone GPD at the leading-order, an argument that requires further theoretical investigation. In Sec. IV 2, we will study numerically three definitions of quasi-GPDs: -Definition in the symmetric frame via the γ 0 operator (H s -Definition in the asymmetric frame via the γ 0 operator (H a 0 (A i ; z) , E a 0 (A i ; z)), Eqs. (46)- (47), As previously stated, the three definitions are not equivalent; they differ in terms of the amplitudes that contribute and power corrections. Thus, it is interesting to numerically compare the convergence of these definitions and also get an idea about the relative size of power corrections.

A. Matrix elements parametrization
One of the goals of this work is to calculate in lattice QCD the Lorentz-invariant amplitudes defined in Eq. (35). To this end, we perform a proof-of-concept analysis based on two calculations of the vector matrix elements, as outlined in Sec. II. Two separate calculations are performed, one in the symmetric and one in an asymmetric frame, which allows us to compare the estimates for A i . For self-consistency, in this section we present the setup in Euclidean space, where we use lower indices in P and ∆ to avoid confusion in the equations presented. The notation for the symmetric frame is and for the asymmetric frame, in which all the momentum transfer is assigned to the initial state, is In the above equations, a factor of 2π L (L: spatial extent of the lattice) is included in ∆ 1 , ∆ 2 , and P 3 . As can be seen, the setup corresponds to zero skewness, that is (p i ) 3 = (p f ) 3 = P 3 . Numerical calculations of the matrix elements in the two frames at the same value of the arguments of the A i must be in line with the Lorentz invariance of the A i . Such a numerical confirmation is a highly non-trivial check of the numerical calculations and the underlying equations, which relate the matrix elements and the amplitudes in the two frames. As mentioned previously, the matrix elements are frame-dependent, and, in general, decompose into different sets of A i in the two frames. This is demonstrated in Eqs. (56) - (71) and Eqs. (73) -(88) below 3 . The analysis takes into consideration four matrix elements of the vector operator, that is, γ 0 , γ 1 , γ 2 , and γ 3 . The matrix element of γ 3 is needed only for the extraction of A 2 and A 7 . We note that the operator γ 3 has a finite mixing under renormalization for lattice regularizations with chiral symmetry breaking [32,[140][141][142]. However, for twisted mass fermions, which we use in this work, the mixing is between γ 3 and γ 5 ; the latter has a vanishing physical matrix elements in the forward limit.
To disentangle A 1 -A 8 , we need eight independent matrix elements, which can be obtained using the unpolarized (Γ 0 ) and three polarized (Γ k ) parity projectors defined as The Γ 0 projector corresponds to an unpolarized proton, while the three Γ k to a polarized projector in the k direction. The parity projectors are applied to the right-hand-side of Eq. (35), along with the spinor normalization. Finally, a trace is taken to mimic the procedure of extracting the lattice matrix elements, that is with the following normalization for the spinor sum The trace is performed analytically, and the obtained expressions correspond to the decomposition of the matrix elements, with the ground state denoted by Π µ (Γ κ ). The matrix elements, Π s/a µ (Γ κ ), for each operator γ µ and projector Γ κ combination are given in Eqs. (56) -(71) for the symmetric frame (superscript s), and Eqs. (73) - (88) for the asymmetric case (superscript a). For simplicity of the presentation, we adopt the expressions at zero skewness. The general equations for ξ = 0 can be found in Appendix C.
K is a kinematic factor that depends on the normalization of the proton state, In fact, K takes a simpler form in the symmetric frame, that is 2m The above equations are elegant, which can be attributed to the zero skewness and the simplification of K. A general feature of the set of equations is that some of them express physically-equivalent matrix elements corresponding to momentum transfer along the two different transverse directions. For instance, Eqs. (57) and (58)  The decomposition of Eq. (35) in the asymmetric frame defined in Eqs. (50) -(51) leads to more complicated kinematic coefficients mainly because E i = E f . Also, the kinematic factor K is canceled by the coefficients of A i due to its simple structure, leading to more elegant expressions. The parametrization for each operator and parity projector is The matrix elements above involve more A i compared to the symmetric frame, and, for instance, Eqs. (73) - (75) contain A 1 , A 3 , A 4 , A 5 , A 6 , and A 8 . Consequently, (A 1 , A 5 , A 6 ) and (A 3 , A 4 , A 8 ) are not decoupled in the asymmetric frame, contrary to the symmetric one. In fact, in the symmetric frame, the coefficients of A 3 , A 4 , A 8 in the γ 0 matrix elements vanish due to E i = E f . Similarly, A 4 , A 5 , A 6 , drop out the matrix elements of γ 1 and γ 2 except for the projector Γ 3 . The above simplifications lead to the aforementioned decoupling and to more compact expressions in the symmetric frame. Similarly to the symmetric frame, A 2 and A 7 appear only in the decomposition of Π 3 .
To summarize, Eqs. (56) - (71) and Eqs. (73) - (88), consist the basis of required matrix elements for the calculation of the Lorentz invariant amplitudes, and are used to present numerical results for the A i for the setup of Table I.

B. Decomposition of Lorentz-invariant amplitudes
The Lorentz-invariant amplitudes can be disentangled using the parametrizations given in the previous subsection, which requires the inversions of Eqs. (56) - (71) and Eqs. (73) -(88), in the symmetric and asymmetric frame, respectively. As mentioned previously, the matrix elements of the operators γ 0 , γ 1 , γ 2 are sufficient to disentangle A 1 , A 3 , A 4 , A 5 , A 6 , and A 8 , in both frames. A 1 , A 5 , and A 6 can then be incorporated to the parametrization of the matrix elements of γ 3 to extract A 2 and A 7 .
The analytic inversion of the equations can, potentially, lead to expressions with complicated kinematic coefficients. For simplicity in the presentation, we show the expressions for ∆ = (∆, 0, 0). We use a subscript s and a in the matrix elements to differentiate between the two frames; A i are frame-independent and do not carry such an index. The expressions for the symmetric frame take the form Below, we give A i in the asymmetric frame for ∆ = (∆, 0, 0), which, as anticipated, has more complicated expressions than the symmetric frame.
One of the main motivations of this work is to extract the twist-2 light-cone GPDs for unpolarized quarks, that is H and E. We begin by constructing the quasi GPDs in coordinate space from the γ 0 , H 0 and E 0 , using the Lorentzinvariant amplitudes, A i . With the A i being frame-independent, one can relate H 0 and E 0 to the matrix elements of either frame; this is a powerful relation, as a calculation in the asymmetric frame requires less computational resources (see Sec. III E). With this in mind, and using the definition of Eq. (43) for the quasi-GPDs in the symmetric frame, we adopt Eqs. (44) - (45), which for zero skewness simplify to For simplicity we suppresse the remaining arguments of the quasi-GPDs and the arguments of the amplitudes. It is useful to rewrite Eqs. (105) - (106) in terms of matrix elements in the symmetric frame that, for ∆ = (∆, 0, 0), leads to As expected, Eqs. (107) -(108) are the usual expressions extracted from the matrix elements of operator γ 0 previously used for the unpolarized GPDs [130]. We also find that the z-dependence in the kinematic factors cancels out. While the above observations validate the methodology of relating H 0 and E 0 to A i , using Eqs. (107) - (108) is computationally costly and not optimal for lattice QCD calculations to extract H s 0 (A s i ; z) and E s 0 (A s i ; z). Instead, one may still employ the historically-used symmetric frame definition for the unpolarized quasi-GPDs (Eqs. (105) -(106)), but perform a calculation of the matrix elements in the asymmetric frame for the extraction of the amplitudes A i . We define such a case by H s 0 (A a i ; z) and E s 0 (A a i ; z). Such a possibility is due to the frame invariance of the amplitudes A i . As outlined in Sec. II, the kinematic coefficients between the two frames are related via a Lorentz boost transformation.
An alternative approach is to define the quasi-GPDs in the asymmetric frame using only the matrix elements of γ 0 (Eq. (43)). In such a case, one can relate H a 0 and E a 0 to the A i as Substituting A i in the asymmetric frame for ∆ = (∆, 0, 0) (Eqs. (97) -(104)), we find We note that (H s 0 , E s 0 ) differ from (H a 0 , E a 0 ) due to their Lorentz non-invariant definition. However, in the infinite momentum limit both approach the correct light-cone limit. While working at finite momentum boost, a different functional form in the two frames is found. Here, we compare H 0 and E 0 between the two frames for pedagogical reasons, as exact agreement between them is not anticipated. Theoretically, there is no preference in using H a 0 and E a 0 versus H s 0 and E s 0 . Historically, the latter is employed, and one convenient approach is the extraction of H s 0 (A a i ; z) and E s 0 (A a i ; z) which uses matrix elements in asymmetric frames. Another aspect of this work follows a different approach from the one outlined above. In particular, we propose a Lorentz-invariant quasi-GPDs H and E definition, as given in Eqs. (41) - (42). Such definitions may have faster convergence to light-cone GPDs. However, further theoretical and numerical investigation is required to reach concrete conclusions. The expressions of Eqs. (41) -(42) simplify for zero skewness, that is Being Lorentz invariant, the above definitions are equivalent in both frames, that is, . For completeness, we provide the expressions of H and E using matrix elements in each frame. As above, we use as an example the case ∆ = (∆, 0, 0), which may be written in terms of matrix elements in the symmetric frame or, alternatively in the asymmetric frame We note that the definition of H and E can be interpreted as the construction of a new operator that is a combination of γ 0 , γ 1 and γ 2 . We emphasize that it is important to provide a comparison of the H and E GPDs in the two frames at the same value of t. This requires ∆ a ⊥ = ∆ s ⊥ . Such a relation is P 3 -dependent, but for the values of P 3 employed in this work (see Sec. III E), are numerically similar (t s = −0.69 GeV 2 , t a = −0.64 GeV 2 ).

D. Renormalization and Matching
A schematic structure of the historical definition of a quasi-GPD (say H) is, while, as we will show in Sec. III C, a schematic structure of the Lorentz-invariant definition is, Here, (c 0 , c 1 , c 2 ) are frame-dependent kinematic factors. Specifically, the linear combination c 1 ..γ 1 .. + c 2 ..γ 2 .. is such that it eliminates the additional amplitudes (with respect to the light-cone case) present in c 0 ..γ 0 .. alone, thereby (potentially) projecting the resulting quasi-GPD faster to the light-cone GPD at the leading order, which is Lorentz-invariant. The question that we want to discuss here concerns the strategy to renormalize the linear combination of (γ 0 , γ 1 , γ 2 ) in Eq. (120). Since the UV divergence of the quark bilinear operator is multiplicative and independent of the Dirac Γ matrix [142][143][144], one can just use the renormalization factor for γ 0 to renormalize the quasi-GPD. Besides, since γ 0 and γ 1,2 are free from O(1) operator mixings due to chiral symmetry breaking on the lattice [140,141], this issue is also avoided in the Lorentz-invariant definition. Note that for our numerical results, we will use the matching coefficient from Ref. [121]. It is known that the GPD matching coefficient for the operator γ 0 reduces to that for the corresponding PDF when ξ = 0, even if t = 0 [121,122]. The PDF matching coefficient for γ 0 is for the amplitude A 1 , which is also the only contributing amplitude to the LI definition of the GPD when ξ = 0. Therefore, the matching coefficients for the γ 0 and the LI definitions of the GPDs are equal. We will elaborate this point more, including the general case of ξ = 0, in a forthcoming publication, along with exploring hybrid renormalization [145] for quasi-GPDs, as well as an improved RI/MOM-type [140,146] scheme to eliminate unwanted finite lattice contributions [147].

E. Calculation parameters
We calculate the proton matrix elements of the non-local vector operator containing spatially-separated quark fields, in theẑ direction, connected by a Wilson line. The proton states are momentum-boosted with nonzero momentum transfer between the initial and final state, |N (P i ) and |N (P f ) are the initial (source) and final (sink) state of the proton. The remaining variables are defined previously. We use momentum smearing [148] to improve the overlap with the proton ground state and suppress gauge noise; Ref. [64] demonstrated that the method is essential for non-local operators. More relevant to the present analysis, it was found that the statistical noise is z-dependent and reduces by a factor of 4-5 in the real part, and 2-3 in the imaginary part of the unpolarized GPDs [130]. The matrix element h µ V is extracted from the ratio where C 2pt and C 3pt , are the two-and three-point correlators. τ is the current insertion time and t s is the source-sink time separation; the source is taken at zero timeslice. We extract the ground-state contribution to h µ V from R µ by taking a plateau fit with respect to τ in a region of convergence, indicated by Π j (Γ κ ). For simplicity, the dependence on z, p f , and p i is not shown explicitly in the matrix elements Π j (Γ κ ).
The calculation is performed on a gauge ensemble of N f = 2 + 1 + 1 twisted-mass fermions including clover improvement, and Iwasaki-improved gluons [149]. The quark masses correspond to a pion with mass 260 MeV. The ensemble has a volume of 32 3 × 64 and lattice spacing a = 0.093 fm. Several of the matrix elements beyond the commonly used γ 0 , have small and noisy signal. To keep the statistical noise under control, we use a source-sink separation of t s = 10a = 0.93 fm. The study of excited states via calculations of multiple time separations is beyond the scope of this work. In future precision calculations, we will include excited states studies. In Table I, we give the statistics for the calculation in the symmetric and the asymmetric frame. The Lorentz-invariant amplitudes have definite symmetry with respect to P 3 → −P 3 , z → −z and ∆ → − ∆ (see Appendix B) and, therefore, we combine all data contributing to the same value of momentum transfer squared, t. We remind the reader that for the kinematic setup in the two frames as given in Table I, t s and t a are not the equal, but are sufficiently close to each other for a comparison between the two frames to be meaningful. We emphasize that the asymmetric frame is computationally more efficient, as one can obtain more than one value of ∆ within the same computational cost.  TABLE I. Statistics for the symmetric and asymmetric frame at zero skewness and ts = 10a. NME, N confs , Nsrc and N total is the number of matrix elements, configurations, source positions per configuration and total statistics, respectively.

Matrix elements and Lorentz-invariant amplitudes
In this section, we present selected matrix elements in the two frames. We point out that the matrix elements in the symmetric frame have definite symmetries with respect to P 3 → −P 3 , z → −z, and ∆ → − ∆. However, such symmetries are, in general, not present in the asymmetric frame, which prevents one from taking averages of the matrix elements for ±P 3 , ±z, or ± ∆ before extracting the A i ; the amplitudes have definite symmetries that are outlined in Appendix B. For consistency in the analysis, the averaging over, e.g., ±z is done at the level of A i for both frames.
In Fig. 1, we show the real and imaginary parts of bare Π 0 (Γ 0 ) in both frames for the eight combinations of ±P 3 and ± ∆ given in Table I. Similarly, Fig. 2, Fig. 3 and Fig. 4 show Π 0 (Γ j ), Π j (Γ 3 ) and Π j (Γ κ ) (j, κ = 1, 2, j = κ), respectively, for the P 3 and ∆ combinations that lead to non-vanishing matrix elements. There are several comments and observations for the behavior of the matrix elements. First, Π s µ and Π a µ are not equivalent and, thus, not directly comparable due to their frame dependence. For instance, their frame dependence can be seen in Eqs. (56) -(88) in Sec. III A, where the matrix elements parametrize in different combinations of A i with frame-dependent kinematic coefficients. Also, the numerical values of some of the kinematic factors, e.g., E i , depend on the frame. For example, Π s 0 (Γ 0 ) contains information on A 1 , A 5 , and A 6 , while Π a 0 (Γ 0 ) contains information on A 1 , A 3 , A 4 , A 5 , A 6 , and A 8 . Second, the matrix elements in the symmetric frame have definite symmetries, which are validated in the data shown in Figs. 1 -4. For instance, Eq. (56) has a symmetric real part and an anti-symmetric imaginary part with respect to P 3 → −P 3 at fixed z, which can be traced back to the symmetries of A i (A 1 (−z ·P ) = A 1 (z ·P ), A 5 (−z ·P ) = A 5 (z ·P ), A 6 (−z · P ) = −A 6 (z · P )) combined with the symmetry properties of the kinematic factors. For simplicity, we do not show all arguments of A i . Third, the lack of definite symmetries in the asymmetric frame appears to be a small effect in Π a 0 (Γ 0 ) and Π a 2 (Γ 3 ), and comparable to the current statistical uncertainties. However, the effect in Π a 1 (Γ 2 ) is more significant, particularly when ∆ → − ∆. Finally, some of the matrix elements, e.g., Π a 1 (Γ 2 ), are theoretically nonzero but are suppressed in magnitude. This has implications in the signal of some of the A i , as we will see below.
The next step in the analysis is to decompose the matrix elements into the corresponding Lorentz-invariant amplitudes. The fact that the A i have definite symmetries makes them interesting to isolate and study. This is done for each kinematic setup of Table I (total of eight). The A i from different kinematic setups can be combined according to their symmetries, by symmetrizing or antisymmetrizing with respect to ±P 3 z, based on the findings given in Appendix B. The frame independence of A i is a major advantage of the proposed parametrization, which we observe numerically using our lattice QCD calculation. Such a test should not be understood as a check of the theoretical formulation, but rather a consistency check of the lattice estimates for A i . The extent of agreement in the two frames provides an estimates of systematic effects arising from non-vanishing lattice spacing.
In Figs. 5 -6, we present the bare A i after combining all values of P 3 and ∆. The amplitudes A 2 , A 4 , A 6 , A 7 and A 8 are accompanied with an explicit factor of z n (n = 1, 2) in Eqs. (56) - (88), and therefore, cannot be accessed at z = 0. One may extrapolate their z dependence to obtain A i (z = 0). In the presentation of Figs. 5 -6, we keep the same range in the y-axis for all the amplitudes for a better quantitative comparison.
We find that A 5 has the largest magnitude both in the real and the imaginary parts, followed by A 1 . The remaining amplitudes are found to be very small or negligible, which can be traced back to the small signal for some of the matrix elements, such as Π 1 (Γ 2 ). Overall, we find very good agreement between the two frames for each A i , as expected theoretically. We remind the reader that there is no exact match of the momentum transfer in the two frames (t s = −0.69 GeV 2 , t s = −0.64 GeV 2 ) and small differences may be attributed to the ∼ 7% change in t, as well as unquantified systematic uncertainties. Such a difference between t s and t a is, in general, not an obstacle in

Lorentz invariant and non-invariant quasi-GPDs
In this paragraph, we use the methodology of Sec. III C to calculate the GPDs based on the γ 0 operator (Lorentz non-invariant), H s/a 0 , E s/a 0 , as well as an alternative Lorentz-invariant operator that combines γ 0 , γ 1 , γ 2 , defining H, E. Having the amplitudes A i , one may use them for any definition of quasi-GPDs, as they contain no information on the frame and are interchangeable between the symmetric and the asymmetric frame, as long as one keeps track of the values of P 3 and t that the quasi-GPDs correspond to. The results for P 3 = 1.25 GeV and t s = −0.69 GeV 2 , t a = −0.64 GeV 2 are shown in Figs. 7 -10. In particular, we compare the definitions of H 0 and E 0 , as given in Eqs. (107) -(108) for the symmetric, and Eqs. (111) - (112) for the asymmetric frame. We emphasize that defining H 0 and E 0 through γ 0 is frame-dependent and, thus, H s 0 and H a 0 have a different functional form; similarly for E s 0 E a 0 . Indeed, we find numerically that the real part for both H 0 and E 0 is not in agreement in the two frames (see left plots of Figs. 8, 10); agreement is found in the imaginary part.
Another interesting comparison is for the Lorentz-invariant definitions, H and E, using matrix elements obtained in the symmetric or the asymmetric frame (see right plots of Figs. 8, 10). As expected theoretically, the agreement between the two frames is significantly improved for both H and E. It is natural to also compare H 0 with H and E 0 with E as extracted in each frame. Also in this case, an agreement is not expected at finite P 3 , as the Lorentz-invariant and non-invariant formalism is not the same. However, some similarity is expected because both definitions agree at P 3 → ∞. Such a comparison is shown in Figs. 7, 9, and it is found that, for this kinematic setup, H 0 is fully compatible with H in both frames; in fact, perfect agreement is found in the asymmetric case. An excellent agreement is found between Re[E] and Re[E 0 ] in the asymmetric frame, while in the symmetric frame there is difference. Finally, differences are observed between Im[E] and Im[E 0 ] for both frames. As previously mentioned, these quantities are not expected to be in agreement for finite momentum P 3 . It is also interesting to note that the statistical errors are considerably smaller in E as compared to E 0 . The origin of this behavior is illustrated in Fig. 11, which shows the respective matrix elements separately for the eight equivalent kinematic cases. At least for this choice of P 3 and ∆, the Lorentz-invariant definition improves the statistical quality of the signal, i.e. that these kinematic contaminations introduce additional noise to the extracted quantity. This trend holds for the symmetric frame too. We also note that this effect does not occur, or is strongly limited, in the case of H GPD. Tracing this behavior back even further, the definition of E involves additional matrix elements that subtract the noise present in Π 0 (Γ 1/2 ) (see Fig. 2, particularly the imaginary part). In turn, H 0 is numerically dominated by the less noisy Π 0 (Γ 0 ) (Fig. 1). We remind the reader that, in general, the difference between t a and t s may be responsible for small differences between the quantities An important component of the lattice calculation is the renormalization, which in this work is done in coordinate space using an RI prescription. Since this is a proof-of-concept calculation, we do not focus on the various prescriptions to improve the renormalization, such as the hybrid scheme [145], and reduction of lattice artifacts in the RI estimates [147], or combination of the two. However, we do emphasize that this is an important direction for future systematic studies of GPDs. For compatibility with the matching formalism of Refs. [121,150], we use the standard RI prescription defined on a single renormalization scale, (aµ 0 ) 2 ≈ 1.95, which will also enter the matching equations. We find negligible dependence when varying µ 0 . As discussed in Sec. III D, the appropriate renormalization for H and E is the one of the γ 0 operator, which is valid for both Lorentz-invariant and non-invariant quasi-GPDs. Details on the calculation of the renormalization functions used in this work can be found in Ref. [75]. As an example, in Fig. 12 we compare the bare and renormalized values of quasi-GPDs using the symmetric frame and the two definitions, that is H s/a 0 and H, E s/a 0 and E. The plots demonstrate the challenges related to the renormalization, that is, as z increases the RI prescriptions becomes less reliable. In practice, the value of the renormalization functions increase exponentially due to the linear divergence leading to renormalized functions that do not decay to zero. Such behavior can be seen in H s,R 0 and H R . Re FIG. 12. Comparison of bare and renormalized quasi-GPDs H0, H (left) and E0, E (right) for the real (top) and the imaginary (bottom) part in coordinate space. All cases correspond to the symmetric frame (t = −0.69 GeV 2 ). The renormalized quantities carry a superscript R.

Light-cone GPDs
To extract the light-cone GPDs from the lattice data, one must transform the latter in momentum space to reconstruct their x-dependence. While this is necessary, the Fourier transform from a finite set of quasi-GPD matrix elements is accompanied by the so-called inverse problem 4 , which mainly affects the small-x region. Nevertheless, the moderate-to-large-x region is not sensitive to this inverse problem, thus allowing us to make reliable predictions. In this work, we use the Backus-Gilbert (BG) reconstruction method [152], which uses a model-independent criterion to choose the light-cone reconstructed GPDs from among the infinite set of possible solutions to the inverse problem. The criterion is that the variance of the solution with respect to the statistical variation of the input data should be minimal. We reconstruct the momentum-space distribution by applying BG separately for each value of x. It should be noted that, despite BG being model-independent and better than the naive Fourier transform, there are still limitations due the small number of lattice data sets. In the work presented here, we vary the number of data that enter the reconstruction, that is z max = 7a, 9a, 11a, 13a.
The x dependence of quasi-GPDs for the various definitions is shown in Fig. 13 using the BG reconstruction method with z max = 9a. All definitions for the H quasi-GPD are consistent, with a very mild difference between the definitions of Eq. (105) and Eq. (109) for x < 0.4. Such a difference become negligible after the matching (see, e.g., Fig 18). On the contrary, the E quasi GPD has a noticeable dependence on the definition. More precisely, the results using Eq. (106) and Eq. (110) are in agreement marginally; the agreement improves in the x > 0 region after the matching, as seen in Fig 19. Differences are also observed in Eq. (106) and (Eq. (110)) when compared to the alternative definition of Eq. (114). These differences persist after the matching. Once again, agreement between the various definitions is not expected by construction. The only agreement imposed by theory is the frame independence of  Eq. (114). Indeed, we find that the numerical results are frame independent despite the small difference in the value of t between the two frames.
The final step of the calculation is the application of the matching equations on the x-dependent quasi-GPDs, to connect the lattice data to the light-cone GPDs, as outlined in Sec. III D. We use the one-loop expression of Ref. [121] to relate the quasi-GPDs in the RI scheme at a scale of 1.95 GeV to the light-cone GPDs in the MS scheme at 2 GeV. At zero skewness, the matching coefficient is exactly the same as in the quasi-PDF case [121]. By varying z max , we first investigate the effect of the truncation of the data set entering the reconstruction of the x-dependence. For simplicity, we show the effect in the light-cone GPDs. We find very small dependence between z max = 7a, 9a, 11a, 13a for all quantities calculated in both frames. In Fig. 14, we show the z max dependence of H a 0 and H extracted from the asymmetric frame calculation. Similarly, in Fig. 15, we show E a 0 and E. As can be seen, all values of z max lead to compatible results, with the statistical uncertainties increasing with z max . Hence, we proceed with z max = 9a as the preferred value.
In Sec. IV 2, we compare A i , H 0 , E 0 , H, and E, as extracted from different definitions and frames. It is interesting to present such a comparison in the light-cone GPDs, which is summarized in Figs. 16 -19. In particular, Fig. 16  demonstrates that the Lorentz invariant and non-invariant definitions for the H-GPD lead to the same light-cone GPDs; this holds for both the symmetric and the asymmetric frames. We remind the reader that the two definitions are different and such an agreement is not expected theoretically. The obtained distributions employing the γ 0 operator and symmetric frame definitions of H s 0 and E s 0 coincide with the results of Ref. [130]. Fig. 17 compares E s/a 0 and E as extracted in each frame. Unlike the case of the H-GPD, here we find that the two definitions lead to GPDs that are of similar magnitude and shape, but are not in agreement for most of the x region. Interestingly, the numerical difference between E s 0 and E is more prominent in the symmetric frame. The overall picture in comparing H 0 with H, and E 0 with E is similar to the one for quasi-GPD matrix elements in coordinate space, as presented in Figs

V. SUMMARY AND FUTURE PROSPECTS
Lattice QCD calculations of x-dependent GPDs have so far been defined in the symmetric kinematic frame, which is, however, computationally very expensive to implement. The main complication is that each value of the momentum transfer can only be accessed one at a time, as it appears in both the initial and final states. Furthermore, the symmetric frame requires two separate inversions for two separate momentum smearing at the source and sink. Hence, the current status of GPD calculations is still at the exploratory stage, with a very limited number of values of the momentum transfer, and, consequently, skewness. In this work, we propose a way to resolve this issue via a new parametrization of off-forward matrix elements relevant to GPDs in terms of Lorentz-invariant amplitudes. Specifically, the frame dependence of the matrix elements is absorbed in the kinematic factors of the parametrization, leaving the amplitudes frame independent. Here, we present a lattice QCD calculation of off-forward matrix elements of the non-local vector operators coupled to momentum-boosted proton states. We observe numerically the frame independence of the amplitudes A i , by comparing their estimates as extracted from the symmetric and the asymmetric frame chosen above. Overall, we find very good agreement between the two frames.
A novel aspect of this work is that the applicability of the new parametrization in any frame has major implications in the reduction of the computational cost. Take for instance the fixed-sink sequential inversion approach, and the asymmetric setup used in this work where the momentum transfer is assigned to the initial state, that is p f = P , p i = P − ∆. The computational advantages are two-fold: (a) one can quadruple the number of measurements by adding all permutations of ∆ contributing to the same t; (b) several vectors ∆ may be obtained for a given p f with an overhead of only the contraction cost. The asymmetric frame needs only one inversion corresponding to a single momentum smearing at the source/sink. So there is a factor of two gain in inversion even for a single momentum transfer. Note that for both cases the momentum smearing is optimized for a selected ∆. However, we identify a broad range of values for the momentum smearing parameter in which the signal improvement is close to optimal.
The Lorentz-invariant amplitudes can be related to the quasi-GPDs of H and E in coordinate space. The latter are not uniquely defined and here we focus on three options: (a) definition in the symmetric frame via the γ 0 operator (H s 0 , E s 0 ); (b) definition in the asymmetric frame via the γ 0 operator (H a 0 , E a 0 ); (c) novel Lorentz-invariant definition (H, E). We emphasize once again that the three definitions are not equivalent; they differ by power corrections. The first definition is of particular interest, as it has been used in previous lattice QCD calculations in the symmetric frame. It is still possible to extract the quasi-GPDs in the symmetric frame in a computationally less-costly way by using the Lorentz-invariant amplitudes A i obtained from the asymmetric frame. In such a case, the quasi-GPDs are defined at the value of t corresponding to the asymmetric kinematic setup. The kinematic coefficients of Eqs. (105) -(106) are obtained via a Lorentz transformation for a transverse boost. Another novel aspect in this work is the Lorentz-invariant definition of quasi-GPDs presented in Eqs. (113) - (114). Such a definition should be in agreement between the two frames. We explored this direction in our lattice calculation and our findings confirm such frame independence (see, e.g., right panels of Fig. 8 and Fig. 10).
The proposed parametrization and the introduction of the Lorentz-invariant amplitudes in not limited to the quantities presented in this work. It is a powerful theoretical tool and has a broad range of interesting applicability that extends beyond leading twist. We believe that it has the potential to shape future calculations of GPDs from lattice QCD with a computational cost that is within reach.
Appendix A: Derivation of the Lorentz-invariant parameterization in terms of the Ai amplitudes In this Appendix, we outline the steps used to obtain the Lorentz-invariant amplitudes A i 's that parametrize the position-space matrix element in the vector case. The starting point for a complete parametrization of the vector operator involves considering all possible Dirac bilinears consistent with the Parity constraint (see also next section where we discuss the implications of the Parity constraint), where, A i ≡ A i (z · P, z · ∆, t = ∆ 2 , z 2 ). However, a further reduction in the number of structures is possible as shown in the following: 1. Using the Gordon Identity, Hence, one can drop the term ∝ A 4 .
2. After a multiplication by ∆ µ in Eq. (A2), we find, 3. After a multiplication by z µ in Eq. (A2), we find, Therefore, in the end one is only left with 8 independent structures.
The above expressions remain valid if Parity is applied along with Time-reversal.
So there are two independent Form Factors. The Form Factors are real functions. On the other hand, our decomposition in the local limit reduces to F µ z=0 =ū(p f , λ ) Now, note that for us the A i 's are "complex amplitudes", but for consistency with the local vector operator we must be able to show that the A i 's surviving are real (that is, either the real or its imaginary part survives). Recall that Hermiticity leads to, Then these constraints on these A i 's at z = 0 leads to, This doesn't help fully because we are still left with three A i 's. We must be able to show that we are left with only two A i 's. For this, we turn to Time-reversal transformation to check if it poses any additional constraint on the A i 's. Recall that Time-reversal leads to, Then these constraints on these A i 's at z = 0 leads to, Hence our decomposition is consistent with the local vector current.
Appendix C: Euclidean-space expressions for the traces for any skewness Here we provide the general expressions for the traces for any frame and for any skewness.
• F 0 with unpolarized projector : where the kinematic factor K is defined as, where abcd is a 4-dimensional Levi-Civita tensor.
• F i with unpolarized projector : • F i with polarized projector : • F 3 with unpolarized projector : • F 3 with polarized projector :