Generalized Parton Distributions from Lattice QCD with Asymmetric Momentum Transfer: Axial-vector case

Recently, we made significant advancements in improving the computational efficiency of lattice QCD calculations for Generalized Parton Distributions (GPDs). This progress was achieved by adopting calculations of matrix elements in asymmetric frames, deviating from the computationally-expensive symmetric frame typically used, and allowing freedom in the choice for the distribution of the momentum transfer between the initial and final states. A crucial aspect of this approach involves the adoption of a Lorentz covariant parameterization for the matrix elements, introducing Lorentz-invariant amplitudes. This approach also allows us to propose an alternative definition of quasi-GPDs, ensuring frame independence and potentially reduce power corrections in matching to light-cone GPDs. In our previous work, we presented lattice QCD results for twist-2 unpolarized GPDs ($H$ and $E$) of quarks obtained from calculations performed in asymmetric frames at zero skewness. Building upon this work, we now introduce a novel Lorentz covariant parameterization for the axial-vector matrix elements. We employ this parameterization to compute the axial-vector GPD $\widetilde{H}$ at zero skewness, using an $N_f=2+1+1$ ensemble of twisted mass fermions with clover improvement. The light-quark masses employed in our calculations correspond to a pion mass of approximately 260 MeV.


I. INTRODUCTION
Parton distribution functions (PDFs) play a crucial role in understanding the quark and gluon structure of strongly interacting systems [1].These functions, measurable in processes such as inclusive deep-inelastic lepton-nucleon scattering, provide valuable insights into the distribution of partons within hadrons as a function of their momentum fraction, denoted as x.PDFs are defined through matrix elements of bi-local operators, where the parton fields are separated by a light-like interval, and the operators are evaluated for the same initial and final hadron states.Generalized parton distributions (GPDs) extend the concept of PDFs by considering light-like parton operators computed for different initial and final states [2][3][4].GPDs introduce additional dependencies on the longitudinal momentum transfer (ξ) and the invariant momentum transfer (t) to the target, in addition to the parton momentum fraction (x).While this multi-variable nature makes GPDs more complex, they offer a wealth of information beyond PDFs.In particular, GPDs provide three-dimensional images of hadrons [5][6][7][8], enable access to the angular momenta of partons [3], and offer insights into the pressure and shear forces within hadrons [9][10][11].Recently, it has been discovered that GPDs exhibit chiral and trace anomaly poles, which provide insights into phenomena such as mass generation in QCD, chiral symmetry breaking, and confinement [12][13][14].Understanding these imprints can offer valuable insights into fundamental aspects of QCD.We also refer the reader to several other review articles that extensively discuss the physics of GPDs [15][16][17][18][19][20][21][22].
In our recent publication of Ref. [85] 1 , we achieved significant advancements in enhancing the computational efficiency of lattice QCD calculations for off-forward matrix elements.The work was also extended to calculate the Mellin moments of the unpolarized GPDs [89], including high moments enabling a physical picture of quark distribution in the transverse plane.In these calculations, we employed a unique approach using asymmetric frames, which differ from the more commonly used symmetric frames.In this approach, the entire momentum transfer ∆ is applied to the initial state (source) of the nucleon.This choice not only reduces the computational cost but also offers the advantage of covering a broader range in t ≡ ∆ 2 , enabling us to effectively map the GPDs across a larger t-space.In our previous work, our focus was on unpolarized quark GPDs (H and E) at zero skewness.We introduced a novel Lorentz-covariant parameterization for the vector matrix element in terms of Lorentz-invariant amplitudes.This also allowed us to establish connections between matrix elements from any two kinematic frames.Additionally, we employed this amplitude-based approach to propose a frame-independent definition of quasi-GPDs and demonstrated that these definitions can potentially result in reduced power corrections in the matching relations to light-cone GPDs.In this work, we extend the amplitude-based approach to compute the axial-vector GPD H at zero skewness.For a comprehensive discussion on the inaccessibility of E at ξ = 0, we refer to Sec.II, where we discuss the intricacies and reasons behind this limitation.
The paper is structured as follows.In Section II, we begin by presenting the definitions of axial-vector light-cone and quasi-GPDs.We then shift our focus towards discussing the Lorentz-covariant decomposition of axial-vector matrix elements in terms of the Lorentz-invariant amplitudes.Furthermore, we establish the relations between these amplitudes and the GPDs H and E. Based on these amplitudes, we propose a few potential candidates for a new, frame-independent definition of quasi-GPDs under the constraints of finite boost momentum.We thoroughly explore the interpretations of these new definitions, carefully examining the subtleties involved while also addressing the important issue of uniqueness/non-uniqueness in their formulation.In Section III, we provide the Euclidean decompositions of lattice-calculable matrix elements in terms of these amplitudes.We also outline our lattice setup for the calculations in position space.Section IV is dedicated to our numerical results, accompanied by a detailed comparison between the symmetric and asymmetric frames at different stages, both in coordinate space and momentum space.Notably, we provide numerical results for the invariant amplitudes and the twist-2 light-cone GPD H, specifically for ξ = 0. Finally, in Section V, we conclude our findings and discuss potential future prospects for further research and exploration in this field.

II. STRATEGY OF FRAME TRANSFORMATION
Computing GPDs in the symmetric frame presents significant challenges in lattice QCD.Extracting a range of momentum transfers requires separate calculations for each ∆, severely limiting the accessible momentum transfer range.This prompts the question of calculating GPDs in computationally advantageous asymmetric frames.One approach, as outlined in our previous work [85], establishes a connection between the symmetric and asymmetric frames through a suitable Lorentz transformation.Employing a Lorentz transformation along the z-direction does not work since a spatial operator distance will receive a nonzero temporal component, which cannot be dealt with in lattice-QCD calculations.In contrast, transverse Lorentz transformations ("transverse boosts") preserve the spatial operator distance.In our second approach, we have developed a Lorentz covariant formalism that allows calculations in any frame.By parameterizing the relevant matrix element using Lorentz-invariant (frame-independent) amplitudes, we establish connections between different frames.In the following sections, we will explore this approach and its implications for computing axial-vector GPDs in asymmetric frames.

A. Definitions of GPDs
To begin, let us revisit the definition of light-cone quark GPDs for a spin-1/2 hadron.In position space, GPDs characterize non-local quark field matrix elements, which are defined as follows: where Γ represents a gamma matrix.The gauge invariance of this correlator is ensured by the presence of the Wilson line In Eq. ( 2), the parameter g represents the strong coupling constant, while A + denotes the plus-component of the gluon field on the light cone.The initial (final) hadronic state in Eq. ( 1) is described by its 4-momentum p i (p f ) and helicity λ (λ ′ ).We introduce the following kinematic variables: the average 4-momentum of the hadrons P , the (aforementioned) 4-momentum transfer ∆, the skewness ξ (which is defined for hadrons with a large light-cone plusmomentum and represents the longitudinal momentum transfer to the hadron), and the (aforementioned) invariant squared 4-momentum transfer t, We use the definitions in Eq. ( 3) in both the symmetric and asymmetric frames.At twist-2, the correlator with Γ = γ + γ 5 in Eq. ( 1) can be characterized by two distinct axial-vector GPDs, H and E. In position space, the expression is given by [16] For the expression corresponding to Eq. ( 1) in momentum space, the Fourier transform is taken with respect to P • z while keeping P + fixed, leading to Now, let us redirect our focus to quasi-GPDs, which are defined in position space through the equal-time correlator [61] Here, the Wilson line is given by For Γ = γ 3 γ 5 , one finds with the quasi-GPDs H 3 (z 3 , ξ, t; P 3 ) and E 3 (z 3 , ξ, t; P 3 ).Eq. ( 8) is the quasi-GPD counterpart of Eq. ( 4).For the expression corresponding to Eq. ( 8) in momentum space, we perform a Fourier transform with respect to P • z while keeping P 3 fixed, yielding In Ref. [90], it was proposed to define the quasi-counterpart of the light-cone GPDs H and E using Γ = γ 3 γ 5 , as presented in Eq. ( 8).The rationale for selecting γ 3 γ 5 instead of, for instance, γ 0 γ 5 is the absence of mixing with other operators under renormalization, where this mixing is regarded as a lattice artifact caused by chiral symmetry breaking [90].Furthermore, Ref. [72] argues that, for this definition, it becomes necessary to substitute γ + γ 5 /P + with γ 3 γ 5 /P 0 in the prefactor of H to ensure consistency with the forward limit.The definition in Eq. ( 9) also produces the correct local limit when integrated with respect to x.

B. Parameterization of an axial-vector matrix element
Now, let us discuss the Lorentz-covariant decomposition of the axial-vector matrix elements, specifically Eq. ( 1) with Γ = γ µ γ 5 , for spin-1/2 particles in position space.By incorporating parity constraints, we establish that the axial-vector matrix element can be expressed as a combination of eight distinct Dirac structures, each multiplied by a corresponding Lorentz-invariant amplitude.The choice of basis for the amplitudes is not unique, and here we employ where ϵ µP z∆ = ϵ µαβγ P α z β ∆ γ .We note that the above equation holds for a general value of z and has a smooth z → 0 limit.The amplitudes A i are frame-independent, while the basis vectors are generally frame-dependent.Note also that the basis vectors in Eq. (10) do not contain factors of z2 .(For example, such factors can occur when working with an orthogonal set of basis vectors.)For the amplitudes, we adopt the concise notation A i ≡ A i (z • P, z • ∆, ∆ 2 , z 2 ) for brevity.The procedure of deriving these results closely follows the steps outlined in Ref. [91], with a similar treatment found in Ref. [92] where the matrix element was parameterized in momentum space using a straight Wilson line.It is worth noting that the number of amplitudes is the same as for the vector current (Γ = γ µ ) discussed in Ref. [85].
While it is possible to work with alternative sets of basis vectors, the number of independent amplitudes will remain unchanged, requiring eight independent lattice matrix elements to disentangle all the amplitudes.Furthermore, in Appendix B, we present a comprehensive discussion on the symmetry properties of the amplitudes implied by Hermiticity and the time-reversal transformation.In particular, the relations in Eq. (B3) imply that the amplitudes A 3 , A 4 , and A 8 are odd in ξ.This symmetry behavior, plus the requirement of a well-defined forward limit of the matrix element in Eq. ( 10), leads us to conclude that these three amplitudes vanish for ξ = 0.In the analysis of the lattice data presented in this work, we first kept those amplitudes as nonzero and indeed found them numerically to be compatible with zero.More discussion can be found in Sec.IV.
We will now establish connections between the light-cone GPDs and the amplitudes.As certain quantities depend on the kinematic frame, it becomes crucial to differentiate and specify the relevant frame.To achieve this distinction, we employ superscripts s and a to denote the symmetric and asymmetric frames, respectively.We note that in the symmetric frame, the momentum transfer is equally distributed between the initial and final states, while any other distribution is considered asymmetric.After substituting µ = + in Eq. ( 10), we can apply a basis transformation to relate the A i 's in the resulting expression to the GPDs in Eq. ( 4): where the A i 's are evaluated at z 2 = 0. We emphasize that in the aforementioned equations, we have expressed the kinematic variables multiplying the amplitudes using Lorentz-invariant scalars.It is crucial to note that this particular re-writing is unique.Furthermore, as evident now, these equations exhibit the property of Lorentz invariance, which guarantees their validity and applicability across different reference frames.It is worth noting that due to the prefactor ∆ + , E drops out of the parametrization, Eq. ( 4), at ξ = 0. Upon first look the expression for E, Eq. ( 12), seems only valid for ξ ̸ = 0.However, based on symmetry arguments (as discussed before and given in Appendix B), A 3 is odd in ξ and vanishes at ξ = 0.One can therefore reliably determine the zero-skewness limit of A 3 /ξ (and of E) by calculating A 3 for nonzero ξ and extrapolating the r.h.s. of Eq. ( 12) to ξ = 0. Note that a recent work has unveiled the possibility of obtaining a glimpse into E at ξ = 0 by studying a specific twist-3 GPD [93].Now, we shift our focus to quasi-GPDs.As highlighted in Ref. [85], one plausible approach to define the quasi-GPDs is by starting from the Lorentz-invariant light-cone definitions in Eqs.(11) and (12) to incorporate z 2 ̸ = 0. Throughout our discussions, we will refer to this definition as the Lorentz-invariant (LI) quasi-GPD: where now the A i are evaluated at z 2 ̸ = 0 2 .In simple words, this definition of the quasi-GPD is based on the same functional form in terms of the A i as the light-cone GPD (see Eq. (11) and Eq. ( 12)).One finds that this definition is given by an operator which combines (γ 0 , γ 1 , γ 2 )γ 5 , rather than the conventional operator γ 3 γ 5 (further explained in the subsequent paragraph and sections).Because of this operator structure, a different matching coefficient is required compared to the one used for γ 3 γ 5 .Specifically, one can disregard the contributions from the operators γ 1 γ 5 and γ 2 γ 5 as they are relatively suppressed by a factor of 1/(P 3 ) 2 .In our numerical results, we chose to implement the matching for γ 3 γ 5 (see Ref. [71] and also Eq. ( 81)) due to the unavailability of the matching kernel for γ 0 γ 5 in the literature.Consequently, this necessitates a new calculation for both the matching and the renormalization of the γ 0 γ 5 operator, along with addressing its mixing effects.Note also that the difference between results with matching for γ 0 γ 5 and γ 3 γ 5 is finite and is expected numerically to be very small.
We now turn to the set of quasi-GPDs already introduced in Eq. (8).By setting µ = 3 in Eq. ( 10), we can perform a change of basis to transform the resulting expression and establish a mapping between the A i and the quasi-GPDs defined in Eq. (8).The relations are as follows: Similar to the light-cone case, it is important to note that there are no divergences arising from the terms A 3 /ξ and A 4 /ξ as ξ → 0 in Eq. ( 16).We iterate that the reason for the well-behaved ξ → 0 limit is the fact that these amplitudes are odd in ξ and vanish for ξ → 0 (see also Appendix B).Moreover, if one intends to calculate the value of E at ξ = 0 using Eq. ( 16), extrapolation from nonzero ξ values becomes necessary.Now, let us examine the frame (in)dependence of Eqs. ( 15)-( 16).It is noteworthy that these equations hold true in both symmetric and asymmetric frames, which is why we have refrained from using explicit "s/a" superscripts to denote the GPDs.The kinematical prefactor of the amplitudes can (again) be uniquely expressed through Lorentz scalars: In our previous study [85], we emphasized the frame-dependence of the conventional definitions of unpolarized quasi-GPDs that employ γ 0 .In contrast, we observe that the helicity quasi-GPDs defined through γ 3 γ 5 are the same in the symmetric and non-symmetric frames.This can be understood since the two frames are connected through a transverse boost which preserves the 3-component.By applying the same reasoning, it can be inferred that quasi-GPDs defined using γ 0 γ 5 will exhibit frame dependence.Furthermore, given the ability to reformulate the traditional definition involving γ 3 γ 5 in a Lorentz-invariant manner, it is clear that it emerges as an additional contender for a Lorentzinvariant definition.Consequently, this example explicitly demonstrates the lack of uniqueness in Lorentz-invariant definitions for quasi-GPDs.The distinction between the Lorentz-invariant definitions presented in Eqs. ( 13)-( 14) and Eqs. ( 17)-( 18) can be attributed to terms proportional to z 2 associated with the amplitudes A 4 and A 7 .
We repeat that the two sets of quasi-GPDs discussed above are not equivalent, as they differ in the contributing amplitudes and both explicit and implicit power corrections.Henceforth, when referring to power corrections, we specifically denote corrections that are proportional to z 2 .The additional amplitudes in Eqs. ( 17)-( 18) can be interpreted as contamination arising from explicit power corrections, which could potentially be suppressed by considering higher values of the momentum.One may therefore speculate that Eqs. ( 13)-( 14) converge faster compared to Eqs. ( 17)- (18).However, it is essential to acknowledge that the amplitudes themselves also contain implicit power corrections, so the above statement should be examined case by case.(The presence of additional amplitudes could potentially mitigate the implicit power corrections inherent in the amplitudes stated in Eqs. ( 13)-( 14).)Ultimately, the actual convergence of the different quasi-GPD definitions is determined by the underlying non-perturbative dynamics.Therefore, it is important to perform numerical comparisons to assess the convergence behavior of these definitions and gain insights into the relative magnitude of power corrections in each case.
We conclude this section by briefly discussing the symmetry properties of H (and E) in position space, as these properties play a crucial role in leveraging symmetries to improve statistical precision in lattice calculations.The Hermiticity constraint provides the symmetries of GPDs under the transformation P 3 → −P 3 for a fixed value of z 3 .Similarly, for a fixed P 3 and with ∆ → −∆, the Hermiticity constraint unveils the symmetries of GPDs under the transformation z 3 → −z 3 .Notably, we find that the real part of H satisfies H(−P 3 ) = + H(P 3 ) (and likewise for E).Furthermore, we observe H(−z 3 ) = + H(z 3 ) (and likewise for E).(The imaginary parts of the GPDs satisfy the same constraints as their real parts, with the exception of a negative sign.)Finally, we would like to reiterate that a comprehensive analysis of the symmetries at the level of amplitudes is provided in Appendix B.

III. LATTICE CALCULATION A. Methodology
In this section, we present a synopsis of the methodology for the lattice QCD calculation of proton matrix elements using the axial-vector operator.All expressions here are presented in Euclidean space, where we use lower indices in P and ∆ to avoid confusion in the expressions given previously in Minkowski space.The goal of this calculation is twofold: (a) compare the Lorentz invariant amplitudes extracted from different frames; (b) present results for the H GPD at multiple values of t.We note that our calculation is performed at zero skewness and, thus, E is inaccessible from the matrix elements.The two frames we employ are the symmetric and an asymmetric in which the final state does not contain the momentum transfer In the above equations, a factor of 2π L (L: spatial extent of the lattice) is implied in ⃗ ∆ and P 3 .Note that ⃗ ∆ is the same in both frames, however, differs due to the term containing the energies.In each kinematic frame, we parametrize the lattice matrix elements using the trace where F [γµγ5] is given in Eq. ( 10).We use four parity projectors; the unpolarized, Γ 0 , and the three polarized, Γ k , defined as K is a kinematic factor that has been obtained based on the normalization of the proton state, As we will demonstrate below, the combination of the four projectors and the four directions of the axial-vector operator can disentangle all eight A i for any kinematic setup.
Next, we focus on the ground-state contribution to the matrix elements that we will denote as Π s/a µ (Γ κ ) (µ, κ : 0, 1, 2, 3).We note that the operator γ j γ 5 (j ̸ = 3) has a finite mixing under renormalization for lattice regularizations with chiral symmetry breaking [47,90,94,95].Such mixing is not included in the renormalization analysis in this calculation, as it would require the matrix elements of the tensor operator.However, the effect is found to be small for the twisted mass formulation with a clover term [47,90].The general expressions for Eq. ( 21) in the symmetric frame for zero skewness are Π s 3 (Γ 0 ) = 0 (37) where K simplifies to 2m 2 /(E(E + m)), due to E i = E f ≡ E in the symmetric frame when ξ = 0.It is interesting to observe that A 7 appears only in Eq. (40).In general, we find fourteen nonzero equations, some of which are linearly dependent.For instance, Eqs. ( 26) and ( 27) have the same numerical value besides a multiplicative factor of ∆ 1 and ∆ 2 , respectively.Still, there are eight linearly independent matrix elements that allow one to disentangle all amplitudes A i .Another observation is that the amplitudes A 3 , A 4 , A 8 are decoupled from the other ones.This is an important aspect, as these amplitudes are expected to be zero at ξ = 0 due to theoretical arguments.In Sec.IV, we will comment more about how one can incorporate this information into the analysis.
The trace algebra of Eq. ( 21) in the asymmetric frame of Eq. ( 20) leads to more complicated kinematic coefficients mainly because E i ̸ = E f , as well as the lack of symmetry between p f and p i .At zero skewness, we obtain To summarize our findings, Eqs. ( 25) - (40) and Eqs. ( 41) -( 56) are sufficient to disentangle the A i in the symmetric and asymmetric frame, respectively.This task can be done analytically by inverting the equations, which, however, leads to very complicated general expressions; it is practically more convenient to implement a numerical inversion of the 8 × 8 system for each value of P and ⃗ ∆.Here, we give the expressions for A i using ⃗ ∆ = (∆, 0, 0) as an example.We use a superscript 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 and for the asymmetric frame at ⃗ ∆ = (∆, 0, 0), one obtains The use of the amplitudes A i is a pathway to extracting the quasi-GPDs using lattice data from any kinematic frame.Here, we present two approaches to relate the A i to the quasi-GPDs: (a) the standard γ 3 γ 5 operator (Eq.( 73)); (b) an alternative Lorentz-invariant definition (Eq.( 75)).Our focus is on zero skewness, which only gives access to the H GPD; the kinematic coefficient of E in Eq. ( 8) becomes zero due to the factor ∆ 3 .For the same reason, E does not appear in the parametrization of the matrix elements in the forward limit.In fact, to obtain E-GPD at t = 0, one must parametrize its t dependence, and, similarly, its estimate at zero skewness could be obtained by a fit using ξ ̸ = 0 values.Below, we give the relation between the quasi-GPD of H at zero skewness using the γ 3 γ 5 definition.We note that the standard definition of Eq. ( 73) is Lorentz invariant, and therefore, it is the same in both frames, as discussed in Sec.II B. At zero skewness, one obtains For simplicity, we only show two arguments for H 3 , that is, A i to indicate the frame used in the calculation and z to explicitly show that the relation for quasi-GPDs is given in coordinate space.To keep the expressions compact, we suppress the arguments of the amplitudes A i .It is useful to rewrite Eq. ( 73) in terms of matrix elements in the symmetric frame for the special case ⃗ ∆ = (∆, 0, 0), for which we find As expected, Eq. ( 74) is the usual expression extracted from the matrix elements of the γ 3 γ 5 operator previously used for the helicity GPDs [77].With the A i being frame-invariant, one can use either A s i or A a i in Eq. ( 73); calculating A a i is computationally less costly and, thus, more optimal for lattice QCD calculations.
An alternative approach to extract the light-cone GPDs is through a Lorentz-invariant definition of choice for the quasi-GPDs H and E, as given in Eqs. ( 17) - (18), where H( , by construction.The expression for H at zero skewness simplifies giving For completeness, we provide the expressions of H using matrix elements in each frame.As above, we use as an example the case ⃗ ∆ = (∆, 0, 0) to write H in terms of matrix elements, that is This alternative definition of H can be interpreted as the construction of a new operator that is a combination of γ µ γ 5 with µ = 0, 1, 2, 3, as given in the example of Eqs. ( 76) - (77).We note, however, that in the case of the helicity, the matrix elements γ k γ 5 with k = 0, 1, 2 (k ̸ = 3) have finite mixing in lattice regularization [90], which affects H.In Sec.IV, we will compare the two definitions of H 3 and H and discuss their merits.

B. Computational setup
The proton matrix elements entering Eq. ( 21) (F [γµγ5] ) use a non-local axial-vector operator containing spatiallyseparated quark fields in the ẑ direction.The Wilson line and the momentum boost are also along the ẑ direction.The matrix elements have momentum transfer between the initial and final state, ⃗ ∆ = ⃗ p f − ⃗ p i , and can be written as |N (p i )⟩ and |N (p f )⟩ are the initial (source) and final (sink) states of the proton, while the remaining variables are defined previously.We use momentum smearing [96] to improve the overlap with the proton ground state and suppress gauge noise; Ref. [97] demonstrated that the method is essential for non-local operators.It was also 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 quasi-GPDs calculated in a previous work [77].In addition, we use five steps of stout smearing [98] to the gauge links of the operator with parameter ρ = 0.15, to further suppress gauge noise, as demonstrated in Refs.[99,100].We note that the stout smearing changes both the matrix elements and the renormalization function, but the renormalized matrix elements should remain independent of the stout smearing.Indeed, in Ref. [101] it was examines the effect of the number of stout smearing steps (0, 5, 10, 15, 20) and shows that the renormalized matrix elements are stoutsmearing independent; the test was performed at ⃗ ∆ = 0 and a physical pion mass ensemble.The same conclusions were reached in the case of the gluon PDF [102] calculated using the same ensemble as this work.The matrix element is extracted from the ratio where C 2pt and C 3pt , are the two-and three-point correlation functions.τ is the current insertion time, and t s is the source-sink time separation; the source is taken at zero timeslice.As shown in Table II, we implement all kinematically equivalent momenta that lead to the same value of p 2 i , p 2 f .Thus, to increase statistical accuracy, we average C 2pt for all possible values.We extract the ground-state contribution to F [γµγ5] from R µ by taking a plateau fit with respect to τ in a region of convergence.Here, we indicate the ground state by Π µ (Γ κ ), and their decomposition is given in Eqs. ( 25) - (56).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 a clover term [103].The gluon part of the action is Iwasaki-improved.The volume of the ensemble is 32 3 × 64, and its lattice spacing, a, is 0.093 fm.The quark masses correspond to a pion mass of 260 MeV.The ensemble parameters are given in Table I.Using this ensemble, we obtain the matrix elements at a source-sink time separation of t s = 10a = 0.934 fm, a choice made to effectively manage statistical uncertainties in the matrix elements.The study of excited states via calculations of multiple time separations lies outside the scope of the current project.Details regarding the statistics of the calculation in both the symmetric and asymmetric frames are provided in Table II.In summary, for P 3 = 1.25 GeV we analyze three different values of −t in the symmetric frame.These are complemented by seven values of −t in the asymmetric frame, which significantly enhances computational efficiency.The majority of the values are within the range −t ∈ [0.17 − 1.50] GeV.To examine dependence on the momentum boost, we focus on −t s = 0.69 GeV 2 , where we use three values of P 3 , that is 0.83, 1.25, and 1.67 GeV.Notably, the asymmetric frame offers a computationally advantageous approach, as it enables the acquisition of multiple values of ⃗ ∆ within the same computational cost.To elaborate, while each value of t in the symmetric frame requires a separate calculation, the data production in the asymmetric frame is divided into two groups: one for (±∆ x , 0, 0) and its permutations and another for (±∆ x , ±∆ y , 0) and their permutations.The benefits of acquiring the data summarized in Table II are three-fold, as it allows: -Comparison of results on the three values of P 3 at fixed −t to assess P 3 dependence; -Comparison of the estimates for A i in the two frames using −t s = 0.69 GeV 2 and −t a = 0.65 GeV 2 .
-Extraction of the −t dependence of the GPDs and apply parametrizations.
In comparing A i between frames, it's worth noting that although t s and t a are not precisely identical, they are closely aligned, differing by only 5%, which allows for a meaningful and reliable comparison.

A. Comparison of kinematic frame
In this section, we focus on the setup with P 3 = ±1.25 GeV and ⃗ ∆ = { 2π L (±2, 0, 0), 2π L (0, ±2, 0)}, implemented in both the symmetric and asymmetric frames.This setup gives −t s = 0.69 GeV 2 in the symmetric frame, and −t a = 0.65 GeV 2 in the asymmetric frame.Our first goal is to compare the A i between the two frames, in a similar fashion as our previous work for the unpolarized GPDs [85].Once agreement is established between the A i from the two frames, all data from Table II will be analyzed to extract the −t dependence of GPDs.
Before presenting the matrix elements, we show the ratio of Eq. ( 79) for two representative cases, R 3 (Γ 3 ) and R 2 (Γ 2 ) and P 3 = 1.25 GeV and ⃗ ∆ = 2π L (2, 0, 0).For better clarity, we use the symmetric frame in which the data have less statistical fluctuations.We choose [3a − 7a] for the fit with respect to the insertion time We remind the reader that the data of Figs. 1 -2 are only one of the eight kinematically equivalent cases.These are averaged at the amplitude level according to the symmetry properties of the latter.Below we present selected matrix elements in the two frames.Although direct numerical comparisons of these values may not yield direct physical insights, they do prove invaluable in assessing the signal quality and the extent of symmetry breaking concerning the sign of P 3 and z.In Fig. 3, we show the real and imaginary parts of the bare matrix element Π 3 (Γ 3 ) for the eight combinations of aP 3 = ± 2π L 3 and a ⃗ ∆ = 2π L (±2, 0, 0), 2π L (0, ±2, 0).Similarly, Figs. 4 -6 show Π 0 (Γ 0 ), Π j (Γ j ; ∆ j = 0), and Π j (Γ j ; ∆ j ̸ = 0) (j = 1, 2), respectively.Note that Π j (Γ j ) leads to independent equations for ∆ j = 0 and ∆ j ̸ = 0 (see, e.g., Eq. ( 30)).All the plots presented in this section offer side-byside comparisons between the symmetric and asymmetric frame data.It is essential to underline that the numerical values of matrix elements in these two frames should not be directly compared.This is because the parametrization of a matrix element for a given operator and parity projector differs between frames in terms of the involved A i and their associated kinematic coefficients.For example, Π s 3 (Γ 3 ) contains information on A 2 , A 6 , and A 7 , while Π a 3 (Γ 3 ) decomposes into A 2 , A 3 , A 4 , A 6 , and A 7 , as can be seen in Eqs. ( 40), (56).Comparison of Π 3 (Γ 3 ) from Fig. 3 in the two frames reveals two features: (a) the matrix elements in the symmetric frame are less noisy than in the asymmetric frame; and (b) the variation of the data between the eight different classes of ±P 3 and ±z is smaller in the symmetric frame.The latter is due to the fact that symmetries in the matrix elements with respect to ±P 3 and ±z are only present in the symmetric frame.Nevertheless, the asymmetry in ±P 3 and ±z is found to be small for this kinematic setup in the asymmetric frame.Similar observations hold for Π 0 (Γ 0 ) and Π 1,2 (Γ 1,2 ), which are not shown here.In the light cone limit, the operators γ 0 γ 5 and γ 3 γ 5 are the components of γ + γ 5 , and, thus, lead to the standard helicity GPDs, H and E. This justifies the large magnitude observed in Figs. 3  -4.However, γ 0 γ 5 has finite mixing under renormalization, while γ 3 γ 5 does not [90].The mixing was previously investigated numerically for twisted mass fermions [90,104] and was found that the inclusion of a clover term in the fermion action suppresses the mixing significantly.To extract the amplitudes A i , we utilize the matrix elements obtained from all possible combinations of operators and projectors, as detailed in Section III A. For each A i , we combine the positive and negative values of P 3 , ⃗ ∆, and z according to their respective symmetry properties, as outlined in Eqs.(B3).Upon averaging the data for A i , we proceed to compare their estimated values in the two frames.This comparison serves as a numerical assessment of their agreement, essentially acting as a consistency check for the lattice estimates of A i .The degree of agreement observed between the two frames offers an assessment of systematic effects, such as those stemming from finite lattice spacing, that may affect the results.
, and A a 7 remain consistent, regardless of whether we include A 3 , A 4 , and A 8 in the analysis or not.This consistency provides additional validation for the results obtained.
In Figs. 7 -8 we present a comparison of the amplitudes using our data obtained with P 3 = ±1.25 GeV and ⃗ ∆ = { 2π L (±2, 0, 0), 2π L (0, ±2, 0)} for both frames (−t s = 0.69 GeV 2 and −t a = 0.65 GeV 2 ).This comparison takes into account all eight combinations of ±P 3 and ± ⃗ ∆.Among the amplitudes, we observe that A 5 has the largest magnitude both in the real and the imaginary parts, followed by A 2 .The remaining amplitudes are notably small or negligible, which can be attributed to the small signal for certain matrix elements.Encouragingly, we find a very good agreement between the two frames for each A i up to statistical fluctiations, as expected given their Lorentz-invariant definition.As previously mentioned, the small differences observed may be associated with the approximately 5% discrepancy between t s and t a , as well as potential systematic uncertainties that have yet to be determined.Regarding the amplitudes A 1 , A 6 , and A 7 , they cannot be directly accessed at z = 0 because their associated kinematic coefficients in Eqs. ( 25) -( 56) become zero.Nevertheless, one may perform extrapolations on their z dependence to estimate A i (z = 0).These findings collectively provide valuable insights into the behavior of the amplitudes under various conditions and kinematic setups.An insightful exploration of the amplitudes involves examining their dependence on momentum transfer.To illustrate this, we will focus on the substantial amplitudes, A 2 and A 5 , and analyze their t dependence.The results for these amplitudes are presented in Fig. 9 and Fig. 10, respectively.We use all the data obtained in this work that covers the range −t ∈ [0.17 − 2.77] GeV 2 .Importantly, since the amplitudes are frame invariant, a single function can describe the data from any frame.This feature allows for a direct comparison of the data, ensuring consistency in the analysis of the two frames.Our observations reveal that as −t increases, both the real and imaginary parts of the amplitudes decrease in magnitude.It is noteworthy that, based on our findings, these amplitudes continue to exhibit non-zero values even at −t beyond 2 GeV 2 .However, it is essential to exercise caution in this high-momentum transfer region.The calculations may suffer from systematic uncertainties and higher-twist contamination, rendering this region less reliable for precise conclusions.It is worth noting that the presence of the pion pole discussed in Refs.[12,13] is argued to extend to the GPD E in the isovector channel, where As discussed in Appendix C, the trace of this pole should be evident in A 5 .By comparing the t-dependence of A 5 with, example, A 2 around z = 0, we can infer that the lattice results for A 5 contains a pion pole behavior, and so should its Mellin moments, for example g P .It is worth highlighting that this z = 0 scenario aligns with our previously established anticipations for the pseudo-scalar form factor exhibiting a pole, as discussed in Ref. [106].The significance of this finding is underscored by the considerably steeper rise of A 5 compared to the other amplitudes.For instance, as observed from the plots in this section, in the t range of [−0.65, −0.17], GeV 2 , the rise factor is approximately 2.7 for A 5 , while for A 2 it is around 1.3.Similarly, in the t range of [−2.29, −0.17], GeV 2 , the rise factor is roughly 15 for A 5 , whereas it is about 2.7 for A 2 .These  findings collectively support the notion that the detection of this pole does not necessitate a full calculation of E.
In the near future, we intend to thoroughly investigate this phenomenon in conjunction with the examination of the manifestation of the pole within the z-dependence of the amplitude.Such extension of our analysis aims to provide a deeper understanding of the pole's influence beyond the context of moments to a broader and more holistic understanding of its impact on the physics under consideration.

B. Quasi-GPDs in coordinate space
Our attention now shifts to the quasi-GPDs, which are renormalized in coordinate space using the RI ′ prescription developed and refined in Refs.[90,101,104].We refer the reader to these publications for more details, as well as the previous work of Ref. [85].In Sec.II B, we discussed that the definition of quasi-GPDs is not unique.In this context, we will consider the standard γ 3 γ 5 definition, denoted as H 3 , as well as an alternative definition that is constructed to be Lorentz invariant, termed H.As mentioned above, H 3 also exhibits frame independence, in contrast to the case of the unpolarized GPDs.This frame independence is linked to the fact that the indices of the axial operator align with the direction of the momentum boost, a point we discussed in detail in Section II B. The relations between the quasi-H GPD and the amplitudes are provided in Eq. ( 73) and Eq. ( 75), for H 3 and H, respectively.To enable the comparison of momentum boost dependence, we utilize the data in the symmetric frame at −t s = 0.69 GeV for all cases.Although there is a slight difference in the imaginary part for intermediate values of z at P 3 = 1.25 GeV, overall, the two definitions remain consistent.In terms of the P 3 dependence, we find that, as P 3 increases, the real part approaches zero at smaller values of z.On the other hand, the imaginary part is enhanced for higher P 3 , a feature also observed in PDFs.This trend suggests a strong link between the behavior of quasi-GPDs and PDFs in response to varying value for the momentum boost.To investigate whether the similarity between the two definitions of quasi-H GPD is specific to the particular −t = 0.69 GeV 2 , we extend the analysis to other values of t.We compare the two definitions for all the data in both the symmetric frame (Fig. 12) and the asymmetric frame (Figs. 13 -14).Remarkably, we consistently find agreement between the two definitions across both frames.This agreement is noteworthy because, theoretically, quasi-GPD definitions are not unique, and one might expect variations in the results.The level of agreement observed suggests a robustness in the analysis and interpretation of these GPDs, despite the lack of a unique definition.In Figs. 15 -16, we present the −t dependence of H 3 and H in coordinate space to provide a comprehensive overview.The shapes of these functions exhibit striking similarities, as previously discussed.However, there are notable distinctions worth highlighting.For the real part, at −t = 0.69GeV 2 , H 3 maintains an equidistant relationship from both H 3 (−t = 0.34GeV 2 ) and H 3 (−t = 0.81GeV 2 ).In contrast, H(−t = 0.69GeV 2 ) shifts closer to H 3 (−t = 0.34GeV 2 ).The difference between the pair H 3 (−t = 1.24GeV 2 ), H 3 (−t = 1.38GeV 2 ) and H 3 (−t = 1.52GeV 2 ) is more pronounced compared to the analogous comparison in the H definition.For the imaginary part, H 3 (−t = 0.34GeV 2 ) demonstrates compatibility with H 3 (−t = 0.69GeV 2 ).In the case of H, there is a discernible difference.

C. Light-cone GPDs
The reconstruction of the x dependence of the quasi-GPDs is not unique for a finite number of discrete data, as the standard Fourier transform suffers from the so-called inverse problem3 , which intensifies in the small-x region.In this work, we employ the Backus-Gilbert (BG) reconstruction method [108], which offers a solution to the inverse problem.This approach is based on a model-independent criterion to choose the light-cone reconstructed GPDs from the infinite set of possible solutions to the inverse problem.The criterion employed is that the variance of the solution with respect to the statistical variation of the input data should be minimal.Despite the model independence of the Backus-Gilbert method, its reliable applicability may be limited by the small number of lattice data sets that enter the reconstruction.An approach to assess the consistency and reliability of our reconstruction is to conduct a sensitivity test by varying the number of input data used in the Backus-Gilbert method.Specifically, we test for different values of z max , namely z max = 9a, 11a, 13a.Furthermore, the functions H 3 (z) and H(z) for z > z max are assumed to be zero.This assumption helps establish a condition for the reconstruction and is a practical choice for such cases where data beyond this limit are either unavailable or less reliable.Fig. 17 shows the x dependence of the quasi-GPD for the three values of z max .It is found that, for all z max values, there is compatibility for the H quasi-GPD up to x = 0.7, as well as consistency between z max = 11a and z max = 13a up to x = 1.With that, z max = 11a has been chosen for the quasi-GPD to proceed with the matching to the light-cone GPDs.For the latter, we use the one-loop equations of Ref. [71] at ξ → 0, which we include here for completeness.
x−1 ln y y−1 − 1 y < 0 , In the above equations, q denoted a general quasi-GPD, q is the corresponding light-cone GPD, and f 1 is the matching kernel.Since quasi-GPDs are defined in the RI scheme, the kernel contains the so-called RI counterterms in addition to f 1 .The expressions are lengthy and can be found in Ref. [109].The matching formalism is constructed with the quasi-GPDs being defined in RI scheme at a scale of 1.2 GeV, while the light-cone GPDs are in the MS scheme at a renormalization scale of 2 GeV.We have conducted a series of investigations for quasi-GPDs, which can also be performed for the light-cone cases.In this section, we provide a comparison between the two definitions for H, examining their P 3 dependence where data is available, as well as their behavior with respect to −t.Fig. 18 illustrates the momentum boost dependence for both definitions, H 3 and H, specifically for the symmetric-frame data at −t = 0.69 GeV 2 .We observe some residual dependence between the two, albeit it is important to note that the difference between P 3 = 1.25GeV and P 3 = 1.67GeV is expected to be well within unquantified systematic uncertainties.The behavior of the two definitions for H GPD is remarkably similar, with only minor distinctions between P 3 = 1.25GeV and P 3 = 1.67GeV.It is worth emphasizing that although the two definitions are different, they both possess Lorentz invariance.Consequently, each definition offers a unique function that is applicable in any frame.Turning to Fig. 19, we shift our focus to the −t dependence of H 3 and H at |P 3 | = 1.25 GeV, where we have a substantial amount of data.Numerical values and statistical uncertainties are found to be similar for both definitions.This result is somewhat expected, as the difference between the two definitions is proportionally tied to A 7 , which is found to be very small.The GPDs exhibit a decaying behavior as −t increases, which parallels the behavior of the form factors. Notably, for −t > 1.5 GeV 2 , we observe a negligible dependence on −t, where the GPDs become of similar magnitude.It is important to recognize that this observation is qualitative in nature, as at such values of −t and for P 3 = 1.25 GeV, the lattice results have increased higher-twist contamination.Nevertheless, we have included this data, as it was obtained at no additional computational cost due to the use of an asymmetric kinematic frame.To conclude this discussion, we provide the data of Fig. 19 in a three-dimensional plot to demonstrate both the −t and x dependence of H GPD. For completeness, we show both definitions we explored in this work.

V. SUMMARY AND FUTURE PROSPECTS
This work builds upon recent advancements that enable the extraction of the x dependence of GPDs from lattice QCD matrix elements calculated in any kinematic frame [85].The main motivation is to efficiently compute GPDs across a range of −t values with efficient use of computational resources, a task that proves challenging in the symmetric kinematic frame.Our approach hinges on the decomposition of matrix elements into Lorentz invariant amplitudes, which can then be related to the quasi-GPDs.In this work, we concentrate on the axial-vector case for the proton, providing a detailed framework to extract the helicity GPDs, H and E. To illustrate our methodology, we present a proof-of-concept calculation, where we obtain the matrix element in the standard symmetric frame, as well as an asymmetric frame, in which the momentum transfer is assigned to the initial state of the proton; the parameters are chosen to lead to a very similar value of −t.The calculation is performed at zero skewness, which has the limitation that the E cannot be obtained directly from the lattice data.Our analysis involves a comparison of lattice data for the Lorentz invariant amplitudes, A i , which confirms the theoretical expectations that the amplitudes are frame independent.Such a finding paves the way for a complete asymmetric frame calculation to obtain the light-cone GPD H at multiple values of −t.In our work, we employ two definitions for the quasi-GPDs, namely the standard definition of γ 3 γ 5 , as well as a Lorentz invariant definition that is based on the same functional form in terms of the A i as the light-cone GPDs.Our observations suggest that both definitions yield comparable numerical results.However, it is important to highlight that the constructed Lorentz invariant definition receives contributions from finite mixing of different operators under renormalization due to chiral symmetry breaking [90].
Our final results for the GPDs are presented in the MS scheme at a scale of 2 GeV and are summarized in Figs.19 -20 for various −t values covering the range [0.17, 2.77] GeV 2 .While we have achieved a robust signal, reducing statistical noise further proves to be a formidable challenge for off-forward matrix elements.Considering the fact that systematic uncertainties are still unquantified for GPDs, aiming for high statistical accuracy falls outside the scope of this work.In the future, we anticipate a deeper exploration of systematic uncertainties, including the effects of eliminating excited states, as well as the impact of volume and discretization effects.We also aim to investigate the dependence on the pion mass and other systematic factors tied to the quasi-distribution approach, such as finite momentum boost and limitations of the one-loop formalism.There is also the potential for parameterizing the −t dependence and exploring the impact-parameter space across a broad range of −t values.Another avenue is introducing nonzero skewness, which would facilitate the direct extraction of E from lattice data and allow for extrapolation to ξ = 0.In summary, this study underscores the significant potential for advancing lattice QCD calculations with regard to generalized parton distributions.These advances promise precision calculations that will contribute to developing a framework of global analysis of both current and forthcoming experimental data, in which incorporating lattice data is possible.

5 FIG. 7 .
FIG. 7. Comparison of bare values of A2 and A5 in the symmetric (filled symbols) and asymmetric (open symbols) frame.The real (imaginary) part of each quantity is shown in the left (right) column.The data correspond to |P3| = 1.25 GeV and −t = 0.69 GeV 2 (−t = 0.65 GeV 2 ) for the symmetric (asymmetric) frame.

FIG. 8 .
FIG. 8. Comparison of bare values of z A1, z A6, and z 2 A7 in the symmetric (filled symbols) and asymmetric (open symbols) frame.The notation is the same as Fig. 7. 0

TABLE I .
Parameters of the ensemble used in this work.

TABLE II .
Statistics for the symmetric and asymmetric frame matrix elements are shown.The momentum unit 2π/L is 0.417 GeV.NME, N confs , Nsrc and N total are the number of matrix elements, configurations, source positions per configuration and total statistics, respectively.
FIG.9.The amplitude A2 for all values of −t given in TableII.
FIG.10.The amplitudes A5 for all values of −t given in TableII.
15G.15.The quasi-GPD H3 for several values of the momentum transfer squared, −t in coordinate space.
16G.16.The quasi-GPD H for various values of the momentum transfer squared, −t in coordinate space.