Transversity GPDs of the proton from lattice QCD

We present the first calculation of the $x$-dependence of the isovector transversity generalized parton distributions (GPDs) for the proton within lattice QCD. We compute the matrix elements with non-local operators containing a Wilson line. The calculation implements the Breit symmetric frame. The proton momenta are chosen as $0.83,\,1.25,\,1.67$ GeV, and the values of the momentum transfer squared are $0.69,\,1.02$ GeV$^2$. These combinations include cases with zero and nonzero skewness. The calculation is performed using one ensemble of two degenerate-mass light, a strange and a charm quark of maximally twisted mass fermions with a clover term. The lattice results are renormalized non-perturbatively and finally matched to the light-cone GPDs using one-loop perturbation theory within the framework of large momentum effective theory. The final GPDs are given in the $\overline{\rm MS}$ scheme at a scale of 2 GeV. In addition to the individual GPDs, we form the combination of the transversity GPDs that is related to the transverse spin structure of the proton. Finally, we extract the lowest two moments of GPDs and draw a number of important qualitative conclusions.


I. INTRODUCTION
The current picture on the nucleon structure stems from decades of increasingly precise measurements of form factors (FFs) and parton distribution functions (PDFs), which, in turn, are special cases of more general functions, the generalized parton distributions (GPDs). At a given hard scale Q 2 , GPDs depend on three variables: the longitudinal momentum fraction of the parent nucleon carried by a given parton, x, the square of the four-momentum transferred to the target in a given reaction, t, and on the skewness ξ, which represents the change in the longitudinal momentum fraction induced by the momentum transfer. Physically, GPDs can be seen as correlations between the longitudinal momentum of partons, with a given spin, and their position in the transverse spatial plane of the parent hadron. Together with the transverse-momentum-dependent PDFs, these functions give an overall, three-dimensional, picture of the nucleon, whose comprehension is one of the main goals of the high-energy nuclear physics community.
Although still in their infancy due to constraints in increasing the momentum boost, as well as from systematic effects, such as discretization and volume effects, the lattice QCD calculations in the field of PDFs have advanced enormously. Returning to GPDs, the first perturbative calculation of the matching equations, which relate the distributions with finite momentum boost to the ones with infinite momentum, appeared soon after the original proposal by X. Ji [97,98]. However, it took a few years until the first study of GPDs of the proton within lattice QCD was performed [68]. In that work, the focus was on the chiral-even GPDs for both the unpolarized and helicity case. Here, we extend this work for the chiral-odd GPDs of the proton, following the same methodology as Ref. [68].
The paper is organized as follows. In Sec. II, we outline the general methodology of the GPDs and present the relations between the matrix elements and the GPDs that are needed to disentangle the latter. In Sec. III, we give the lattice details for the isolation of the ground state, the control of statistical uncertainties and the kinematic setup. In separate subsections, we summarize the renormalization procedure, the reconstruction of the x-dependence using the Backus-Gilbert method, as well as the matching formalism. The main results for the matrix elements are presented in Sec. IV, and the final GPDs are given in Sec. V. We also present a comparison with the unpolarized and helicity GPDs. Sec. VI shows our results for the Mellin moments of GPDs and quasi-GPDs, and Sec. VII summarizes our findings.

II. METHODOLOGY
The most computationally expensive aspect of this work is the calculation of the proton matrix elements of non-local operator containing a Wilson line. Without loss of generality, the Wilson line is in the z-direction, W (0, z), which is the same as the direction of the momentum boost for the proton. The operator under study is the tensor with a Dirac structure of the form σ 3j , where j is in the x-or y-direction. Under these constraints, the matrix element reads |N (P i ) and |N (P f ) represent the initial (source) and final (sink) state of the proton labeled by its momentum. We calculate the matrix elements h 1 T and h 2 T separately, because they do not contribute to the same kinematic setup, and can be used as independent equations for disentangling the GPDs (see, e.g., Eqs. (6) - (13)). The matrix elements have dependence on the parity projection, Γ ν , which is implied in the right-hand-side of Eq. (1) for simplicity. We will discuss this below and in Sec. III. Also, in this discussion, we consider h j T as the renormalized matrix element in a given scheme and at a scale µ 0 , entering through the renormalization function, Z T (z, µ 0 ). More details on the renormalization procedure are given in Sec. III B.
GPDs require off-forward kinematics, that is, P f − P i ≡ ∆ = 0. In fact, GPDs depend on the 4-vector momentum transfer squared, t, and not on the individual nucleon momenta. We note that the matrix element h j T depends on the source and sink momenta. In the boosted frame, t is defined as t ≡ −∆ 2 + (E(P i ) − E(P f )) 2 . E(p) is the energy of the proton at momentum p given by the dispersion relation, E(p) = m 2 + p 2 , and m is the mass of the proton.
The standard definition of the light-cone GPDs is in the symmetric (Breit) frame, which requires that P f = P + ∆ 2 and P i = P − ∆ 2 , where P represents the proton momentum boost, P = (0, 0, P 3 ). Besides t, the GPDs have implicit dependence on the momentum transfer in the direction of the boost via the parameter skewness. On the lattice, the relevant quantity is the quasi-skewness defined as The skewness is an important parameter of GPDs, as it separates the x region into two parts, that is the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) region [99][100][101][102], and the Efremov-Radyushkin-Brodsky-Lepage (ERBL) [103,104] region, defined as DGLAP region : x > |ξ| , ERBL region : Each region has a physical interpretation [105]. In the positive-x (negative-x) DGLAP region, the GPDs correspond to the amplitude of removing a quark (antiquark) of momentum p from the hadron, and then inserting it back with momentum p + ∆ (∆: Minkowski momentum transfer). In the ERBL region, the GPD is the amplitude for removing a quark-antiquark pair with momentum −∆. By definition, the ERBL region becomes trivial at ξ = 0.
As mentioned above, the matrix elements depend on the details of the kinematic setup, P f , P i , while the GPDs depend on t and ξ; the remaining dependence on the setup is absorbed into the coefficients of the GPDs that appear in the decomposition. Since there are four transversity GPDs, H T , E T , H T , and E T , one needs four independent matrix elements h j T (Γ ν , z, P f , P i ) to disentangle them; this can be controlled by the choice of the operator (j), parity projector (ν), initial and final momenta. Note that the decomposition is independent of z, and is applied at each value of z separately. The decomposition of the matrix elements is based on continuum parametrizations, which for the transversity case take the following form in Euclidean space where O ≡ Γ νūN (P f , s ) O u N (P i , s) with u N the proton spinors. Also, P = P f +Pi 2 and ∆ = P f − P i . F G plays the role of a form factor, which gives the quasi-GPD of G, G q , once the Fourier transform is taken (G : H T , E T , H T , E T ). The parametrization of Eq. (3), in its general form, is very complicated. Here, we give the relevant expressions for the class of momentum transfer that we use. We apply four different parity projectors, that is, the unpolarized, Γ 0 , and three polarized, Γ k , The first class of momenta we employ is ∆ = (0, q, 0), which correspond to zero skewness. In this case, the initial and final momenta are P i = (0, − q 2 , P 3 ) and P f = (0, q 2 , P 3 ), respectively. For these momenta, we have nonzero contributions from four matrix elements, that is: where C 0 = 2m 2 E(E+m) for zero skewness, and E denotes the energy (E f = E i ≡ E). Note that F H T and F E T are obtained directly from Eq. (6) and Eq. (7), respectively. F E T and F H T are disentangled using Eqs. (8) -(9) together with Eq. (6).
For nonzero skewness, we employ ∆ = (0, q y , q z ) with |q z | = |q y | = q (q > 0), that is, P f = (0, qy 2 , P 3 + qz 2 ) and P i = (0, − qy 2 , P 3 − qz 2 ). For these momenta, we have nonzero contributions from four matrix elements, that is: where , P f z = P 3 + qz 2 , and P iz = P 3 − qz 2 . Unlike the case ∆ = (0, q, 0), here all matrix elements enter in the decomposition of all four GPDs. Since we are using positive and negative values for the momentum transfer, one has to be careful with the signs in the decomposition. In addition to the signs of the kinematic factors, one also has to consider that H T q , E T q , and H T q are even functions of ξ, while E T q is odd [106,107], which also holds for the quasi-GPDs [108] and F G . For example, F E T (z, −ξ, t, P 3 ) = −F E T (z, ξ, t, P 3 ), which is taken into account in the decomposition for the negative value of ξ.
Once F H T , F E T , F H T , and F E T are disentangled from the renormalized matrix elements, we transform them in momentum (x) space to obtain the x-dependence using the Backus-Gilbert (BG) method [109], as described in Sec. III C. This procedure gives the quasi-GPDs, H T q , E T q , H T q , and E T q . Note that we use the subscript q to denote the quasi-GPDs. Finally, the light-cone GPDs H T , E T , H T , and E T are obtained after application of the matching procedure, outlined in Sec. III D.

A. Matrix elements
In this work, we focus on the isovector flavor combination u−d for the transversity GPDs, which requires calculation of only the connected diagram shown in Fig. 1. The matrix elements are constructed from the two-point and threepoint correlation functions, where N α (x) = abc u a α (x) d b T (x)Cγ 5 u c (x) is the interpolating field for the proton, and τ is the current insertion time. Without loss of generality, we take the source to be at (0, 0). The three-point functions are calculated for the up-and down-quark/antiquark fields, (ψ,ψ), which are combined to form the u − d isovector contribution.
N (x, t) N (0, 0) For nonzero momentum transfer, one must form an optimized ratio to cancel the time dependence in the exponentials and the overlaps between the interpolating field and the nucleon states, namely In the limit (t s − τ ) a and τ a, the ratio of Eq. (16) becomes time-independent and the ground state matrix element is extracted from a constant fit in the plateau region, that is In this work, we choose t s = 12a, which was also used in our previous work for the unpolarized and helicity GPDs [68]. An extensive study on the excited-states effect for the forward limit of non-local operators was done in Ref. [58]. In that work, we demonstrated that a source-sink separation above 1 fm is sufficient to obtain the ground state contribution within the reported uncertainties. h j,B T denotes the bare matrix element, while Eq. (1) is the renormalized one. These are related multiplicatively, using the renormalization function, Z T (z, µ 0 ), obtained non-perturbatively, Note that the multiplication is complex. We refer the reader to Sec. III B for more details.
To improve the overlap with the proton ground state, we construct the proton interpolating field using momentumsmeared quark fields [110], on APE-smeared gauge links [111]. The momentum smearing technique is essential to suppress gauge noise in matrix elements with boosted hadrons, and in particular for non-local operators [40]. The momentum smearing approach allowed us to obtain GPDs for protons boosted up to P 3 = 1.67 GeV. Beyond that momentum, it is unfeasible to obtain the matrix element with controlled statistical uncertainties at a reasonable computational cost. This is in agreement with other calculations of non-local operators with boosted hadrons (see, e.g., Table 1 of Ref. [82]). The momentum smearing function S on a quark field, ψ, reads where α G is a parameter of the Gaussian smearing [112,113], U j is a gauge link in a spatial j-direction. P is the momentum of the proton (either at the source, or at the sink) and ξ is a free parameter that can be tuned to achieve maximal overlap with the proton boosted state. For ξ = 0, Eq. (19) reduces to the Gaussian smearing function. The fact that the exponent in Eq. (19) depends on the momentum of the proton state, means that separate quark propagators are needed for every ∆, because the gauge links are modified every time by a different complex phase. In our implementation, we keep ξP parallel to the proton momentum at the source and at the sink. This strategy avoids potential problems due to rotational symmetry breaking. It also has the benefit that every correlator entering the ratio of Eq. (16) is optimized separately. The effectiveness of the momentum smearing has been demonstrated in our previous work for PDFs [40,58], as well as for GPDs [68]. In fact, for the unpolarized GPDs, we found that the statistical noise is suppressed by a factor of 4-5 in the real part, and 2-3 in the imaginary part, depending on the value of z.
The analysis is performed using a gauge ensemble of twisted-mass fermions with a clover improvement, and Iwasakiimproved gluons. The ensemble has two dynamical degenerate light quarks plus a strange and a charm quark (N f = 2 + 1 + 1) [114] in the sea. The quark masses have been tuned so that the pion mass is about 260 MeV. The lattice volume is 32 3 × 64 and the lattice spacing a = 0.093 fm. The three-point correlators are obtained for a source-sink separation of t s = 1.12 fm. In Table I, we summarize the statistics for each value of the nucleon momentum boost P 3 , momentum transfer ∆ and t, as well as skewness ξ. The GPDs have definite symmetry with respect to ξ → −ξ, and therefore, we combine the data at +1/3 and −1/3. Below, we also compare with the transversity PDF, h 1 , obtained on the same ensemble with the statistics shown in Table II.

B. Renormalization
The matrix elements h j T are renormalized non-perturbatively with the renormalization function Z T , which is defined in an RI-type scheme at some scale µ 0 . The vertex functions of the non-local tensor operator are calculated using the momentum source method [115,116] that suppresses statistical noise. We work in the twisted basis to calculate the matrix elements h j T , and, therefore, the operator σ 3j (physical basis) renormalizes with the operator γ 5 σ 3j . We apply the following condition to the vertex functions of γ 5 σ 3j at each value of z separately, We also calculate the quark propagator that is needed for the quark field renormalization, Z q , V(p, z) (S(p)) is the amputated vertex function of the operator (fermion propagator) and S Born (p) is the tree-level of the propagator. In the prescription of Eq. (20), the vertex functions are projected with the so-called minimal projector, which defines Z T (z). The use of this definition is necessary, as the matching formalism is only known for this scheme [117]. Note that we use the symbol Z T for the renormalization function prior to taking the chiral limit of Eq. (24). Similarly for the fermion field renormalization, Z q . In a nutshell, the vertex function of the tensor non-local operator σ 0l (l = j = 3 = l) 1 with a Wilson line in the z-direction contains contributions from three structures, that (see, e.g., Eq. (76) of Ref. [118]). The projector is defined such that it isolates the tree-level contribution of the vertex function, S 1 . For the operator under study, we use the projector In this work, we calculate the vertex functions with the Wilson line in all spatial directions projected with the equivalent P T . Since S 1 is independent of the direction of the Wilson line, we average over the three directions. Numerically, we find that the estimates of Z T (z) using the minimal projector are similar to the estimates obtained by projecting with the tree-level value. This is an indication that the contamination from S 3 in the vertex function is small 2 .
The prescription of Eq. (20) is mass-independent, and therefore Z T should not depend on the quark mass. However, there might be residual cut-off effect of the form am q . To eliminate any systematics related to such an effect, we extract Z T using five degenerate-quark-mass ensembles (N f = 4) with the same lattice spacing as the ensemble we use for h j T . These N f = 4 ensembles correspond to a pion mass in the range 350 -520 MeV. The estimates of Z T from each ensemble are used for a chiral extrapolation. More details on this procedure can be found in Ref. [58].
Z T is scheme-and scale-dependent, and therefore, is defined at some RI scale µ 0 . We use several values of µ 0 , chosen to be isotropic in the spatial directions, which suppresses discretization effects. Furthermore, the vertex momentum is such that the ratio p 4 (p 2 ) 2 is less than 0.35 [119]. In this work, we use different values of µ 0 ((a µ 0 ) 2 ∈ [1, 5]) to check the dependence of the matching formalism on µ 0 . For each value of µ 0 , we apply a chiral extrapolation using the fit to extract the mass-independent Z RI T (z, µ 0 ). For our final results, we use the renormalization functions defined on a single RI renormalization scale, (aµ 0 ) 2 ≈ 2.57. This scale also enters the matching equations, which connect the quasi-GPDs in the RI at a scale of µ 0 to the GPDs in the MS at a scale of 2 GeV. We find negligible dependence in the final GPDs when varying the initial scale µ 0 in the quasi-GPDs.

C. Reconstruction of x-dependence
The quantities F G 3 , where G = H T , E T , H T , E T , are related to the quasi-distributions, G q , via a Fourier transform, as the latter are expressed in momentum space, Therefore, extracting the quasi-GPDs requires integration over a continuum range of z, while the lattice provides only a discrete set of determinations of F G , for integer values of z/a up to roughly half of the lattice extent in the direction of the boost, L/2a. Thus, obtaining the quasi-GPDs, or for that matter any x-dependent distributions, poses a mathematically ill-defined problem, as discussed in detail in Ref. [50]. The inverse problem originates from incomplete information, i.e. attempting to reconstruct a continuous distribution from a finite number of input data points. As such, its solution necessarily requires making additional assumptions that provide the missing information. These assumptions should be mild and preferably model-independent -else, the reconstructed distribution may be biased.
In this work, we use the Backus-Gilbert (BG) method [109], which was also proposed in Ref. [50]. The method relies on a model-independent criterion to choose from among the infinitely many possible solutions to the inverse problem, namely that the variance of the solution with respect to the statistical variation of the input data should be minimal. While the BG method is superior to the naive Fourier transform, there are limitations due the small number of the lattice data, and the BG would be improved if a larger volume and finer lattice spacing ensemble is used. The reconstruction is done separately for each value of x. In practice, we separate the exponential of the Fourier transform into its cosine and sine parts, related to the real and imaginary parts of the matrix elements, respectively. We define a vector a K (x), where K is either the cosine or sine kernel, of dimension d equal to the number of available input matrix elements, i.e. d = z max /a + 1; the matrix elements for z beyond z max are neglected, assuming that they are approximately zero within uncertainties. The BG procedure consists in finding the vectors a K (x) for both kernels according to the variance minimization criterion. The vector a K (x) is an approximate inverse of the cosine/sine kernel function K(x), that is where are elements of a d-dimensional vector of discrete kernel values corresponding to integer values of z/a entering the reconstruction. Therefore, the function ∆(x − x ) is an approximation to the Dirac delta function δ(x − x ). The quality of this approximation depends on the achievable dimension d at given simulation parameters.
The vectors a K (x) are identified from optimization conditions based on the BG criterion. For more details, see Ref. [50]. Below, we summarize the methodology. We define a d × d-dimensional matrix M K (x), with matrix elements where x c is the maximum value of x for which the quasi-distribution is taken to be non-zero (i.e. its reconstruction proceeds for x ∈ [0, x c ]). The parameter ρ regularizes the matrix M K . This regularization was proposed by Tikhonov [120], and is suggested as a possible way to make M K invertible [50,121,122]). The value of ρ determines the resolution of the method and should be taken as rather small, in order to avoid a bias. We use ρ = 10 −3 , which leads to reasonable resolution and is large enough to avoid oscillations in the final distributions related to the presence of small eigenvalues of M K . We have checked that the dependence on ρ is negligible. Additionally, we define a d-dimensional vector u K , with elements Applying the aforementioned optimization conditions leads to Finally, the BG-reconstructed quasi-distributions are given by

D. Matching Procedure
Following the reconstruction of the x-dependence of the quasi-GPDs, we proceed with obtaining the light-cone GPDs. Contact between the physical GPDs and the quasi-GPDs is established through a perturbative matching procedure. The general factorization formula reads where C G is the matching kernel and is known to one-loop level in perturbation theory. The involved renormalization scales are: µ 0 -RI renormalization scale, its z-component (µ 0 ) 3 (with r = µ 2 0 /(µ 0 ) 2 3 ), and µ -final MS scale; here we choose µ = 2 GeV. This formula establishes that quasi-distributions are equal to light-cone distributions up to powersuppressed corrections (nucleon mass (m) corrections and higher-twist corrections). The matching coefficient for the GPDs was first derived for flavor non-singlet unpolarized and helicity quasi-GPDs in Ref. [97] and for transversity quasi-GPDs in Ref. [98], using the transverse momentum cutoff scheme. Recently, a matching formula was also derived for all Dirac structures [117] relating quasi-GPDs renormalized in a variant of the RI/MOM scheme to MS light-cone PDFs (minimal projector of Eq. (23)). In these calculations, it was shown that the matching for GPDs at zero skewness is the same as for PDFs. It was also demonstrated that, to one-loop level, the H-type and E-type GPDs have the same matching formula. The matching kernel for the transversity GPDs and parton momentum p 3 reads The functions G 1 , G 2 , G 3 for the matching of bare quasi-GPDs can be found in Ref. [117], while the one-loop RI counterterm f P T for the RI/MOM variant that we employ (minimal projector, P T ) is given in Ref. [49]. The plus prescription is defined as and it combines the so-called "real" (vertex) and "virtual" (self-energy) corrections. Below we give the expressions for the functions G i for completeness:

IV. NUMERICAL RESULTS
We begin our presentation with the bare matrix elements for the ground state, as extracted from Eq. (17). In Fig. 2, we plot the four matrix elements contributing to ∆ = 2π 32 (0, 2, 0) (t = −0.69 GeV 2 ), that is Eqs. (6) -(9). We compare the signal for the three values of P 3 employed. As expected, the statistical uncertainties increase with the momentum. We find that the matrix elements of h 1 T (Γ 2 ) have the most dominant contributions in both the real and imaginary parts, followed closely by h 2 T (Γ 1 ). We remind the reader that h 1 T (Γ 2 ) is directly related to the leading H T -GPD, see Eq. (6). h 2 T (Γ 0 ) has a smaller signal than the above matrix elements, but it is clearly non-negligible. On the contrary, the matrix element contributing to E T , h 1 (Γ 3 ), has negligible contribution for both the real and imaginary parts, with the exception of the real part for P 3 = 1.67 GeV, which slightly deviates from zero. We note that a convergence with respect to P 3 is not necessarily anticipated in the matrix elements, but rather at the level of the final matched GPDs. As can be seen from Eqs. (6) - (13), there is a dependence on the kinematic setup, which includes P 3 through the energies, and in some cases, directly. We remind the reader that the matching also contains the momentum boost P 3 .   (9). Results correspond to t = −0.69 GeV 2 , ξ = 0 and three nucleon momenta. The color notation is the same as in Fig. 2.
Upon renormalization of the matrix elements of Fig. 2, we disentangle the four F G that will be eventually matched to each transversity GPD. We demonstrate the dependence of F G on P 3 in Fig. 3. For z = 0, F G are independent of P 3 ; this does not hold for z = 0 due to the breaking of Lorentz invariance. In fact, the values at z = 0 correspond to the tensor form factors, which are the lowest moments of the transversity GPDs. Further discussion can be found in Sec. VI. Focusing on the highest momentum, we find signal for F H T , F E T and F H T . As expected from the behavior of h 1 (Γ 3 ), F E T is suppressed compared to the other ones. The imaginary part of F E T and F H T is also zero within uncertainties.
It is interesting to compare the matrix elements contributing to H T for different values of the momentum transfer. In Fig. 4 we show h 1 T (Γ 2 ) at P 3 = 1.25 GeV for −t = 0 , 0.69, 1.02 GeV 2 . For the case of ξ = 0 (−t = 0 , 0.69 GeV 2 ), the matrix element is proportional to F H T , while for ξ = 0 (−t = 1.02 GeV 2 ) it receives contributions from F E T and F E T . The most notable feature of h 1 T (Γ 2 ) is the lowering of its value with the increasing of −t for both the real and imaginary part. The real part becomes compatible with zero at z/a = 9, 8, 6 for −t = 0 , 0.69, 1.02 GeV 2 , respectively. For the imaginary part, we find that compatibility with zero is at z/a = 14, 12, 8 for −t = 0 , 0.69, 1.02 GeV 2 , respectively. The decomposed renormalized F G are shown in Fig. 5 for |ξ| = 1/3 and −t = 1.02 GeV 2 . F E T and F H T are the most dominant contributions in the matrix elements, followed by F H T . The ξ-odd F E T is compatible with zero. We also find that F H T (F E T ) has the highest (lowest) signal-to-noise ratio. Based on these results, we expect that the final E T -GPDs will have a signal compatible with zero, and E T will have enhanced statistical uncertainties as compared to H T .   V. x-DEPENDENCE OF GPDS As mentioned in Secs. II -III, the quasi-distribution approach relates the lattice data at a given value of the momentum boost to the light-cone GPDs. Therefore, the final light-cone GPDs should be momentum-independent. Practically, this argument is not exact, because the matching kernel is known only to one-loop level, and there are systematic effects, such as higher-twist contamination. In this work, we use three value of P 3 to check for convergence in the final GPDs with respect to the momentum boost. Choosing the right value for the cutoff z max in the reconstruction of the x-dependence is also an important aspect of the analysis. The criterion is not unique, and one can use the z-behavior of each F G as a guidance. Based on our results, we choose z max such that the functions F G (z max ) become zero. According to this criterion, we find that appropriate choices for H T at P 3 = 0.83, 1.25, 1.67 GeV and ξ = 0 are z max /a = 13, 9, 7, respectively. This holds for both the real and imaginary part. As expected, the increase of P 3 results in a faster decrease of the matrix elements. For the real part of E T and H T , we choose z max /a = 7. Our results for E T and H T indicate that the imaginary part is compatible with zero within errors, and is, hence, neglected. Some fluctuations at large z max are due to the rapid increase of the renormalization functions. For all the GPDs at |ξ| = 1/3, we use z max /a = 7 for both the real and imaginary parts. We remind the reader that the distribusions at nonzero skewness, G(x, 1/3, t) as already been combined with G(x, −1/3, t), which is symmetric for the three GPDs we show here.
The convergence of H T is shown in the left panel of Fig. 6 for ξ = 0 and t = −0.69 GeV 2 . The bands include only statistical uncertainties. We find that convergence is achieved for the two highest values of P 3 , implying that the reconstructed H T is momentum-independent even when the matrix elements have a momentum boost of 1.25 GeV. This conclusion is based on the current statistical uncertainties and the one-loop truncation of the matching formalism. For H T , a momentum of 0.83 GeV is also compatible with the higher momenta up to around x = 0.4. Beyond that point, the distribution is lower than its value for the higher momenta. In the right panel of Fig. 6, we compare, at the highest momentum P 3 , H T (x, 0, −0.69 GeV 2 ) with its forward limit, h 1 (x). We find that for the small and intermediate x region, h 1 (x) is higher than H T , which is expected. After x = 0.4 the two distributions are compatible. The same equality seems to hold numerically for the whole anti-quark region. The large-x behavior of PDFs and GPDs for the unpolarized case has been studied using a power counting analysis [123]. While similar arguments do not exist for the transversity GPDs, our data indicate that there is no t dependence for x → 1, similar to the unpolarized H-GPD, but unlike the helicity H-GPD [68].  At P 3 = 1.25 GeV, we have results for both zero and nonzero skewness, which are compared in Fig. 7. In the ERBL region, there is a significant decrease of the distribution as −t increases. However, the distribution in the DGLAP region shows less sensitivity in t. We note that the discontinuity at x = ±ξ is not physical, as twist-2 GPDs are continuous functions at the boundaries of the ERBL region [52,108]. The observed effect is due to uncontrolled higher-twist contamination, which cannot be treated by the matching formalism as it contains only the leading twist.
The data for E T and H T at ξ = 0 are shown in Fig. 8 for the two higher momenta. For these GPDs, we do not show results for P 3 = 0.83 GeV, as the matrix elements for F E T and F H T do not decay to zero. This is an indication that a boost of 0.83 GeV is not large enough. As expected from the decomposition of the matrix elements in coordinate space, the uncertainties on these quantities are significantly enhanced compared to H T . Thus, one will need considerably larger statistics to address them in the future. At the present stage, the qualitative conclusion that can be drawn is the approximate symmetry between the quark and antiquark regions, originating from the imaginary part of the respective matrix elements being compatible with zero (see Fig. 3). This also implies a much larger magnitude of the antiquark part for these two GPDs as compared to H T . We also find that H T is negative. Similar qualitative conclusions are observed in the scalar diquark model of Ref. [108]. Comparing the distributions for the two momenta, we find compatibility within the large uncertainties.    9 compares E T (left) and H T (right) at the same boost of 1.25 GeV for zero and nonzero skewness. The behavior is similar to H T , that is, the increase of −t reduces the magnitude of the distributions, and the introduction of skewness leads to nonphysical discontinuities at x = ±ξ due to higher-twist effects. However, due to the large uncertainties in E T and H T , the function left and right of the boundaries x = ±ξ appears to be continuous within uncertainties. Here we do not show E T , as the signal is weak and zero within uncertainties (see, e.g., Fig. 5).
We also explore the combination E T + 2 H T which is related to the transverse spin structure of the proton, and is considered a more fundamental quantity than E T [7]. The E T + 2 H T combination has the physical interpretation of the lateral deformation in the distribution of transversely polarized quarks in an unpolarized proton. Also, according to Ref. [8], the lowest Mellin moment (n = 0 in Eq. (37)) of E T + 2 H T in the forward limit is the transverse spin-flavor dipole moment in an unpolarized target [8], k T . The first non-trivial moment of E T + 2 H T (n = 1 in Eq. (37)) is related to the transverse-spin quark angular momentum in an unpolarized proton. In Fig. 10, we show the combination E T + 2 H T for P 3 = 1.25 GeV at zero and nonzero skewness. Our results for the two values of ξ are compatible within uncertainties, which are rather large. We do find that the distribution for |ξ| = 1/3 tends to be systematically lower than the one at ξ = 0, but further study is needed to control the uncertainties and reach more meaningful conclusions. Since we have results for the unpolarized and helicity GPDs on the same ensemble and kinematic setup [68], it is interesting to compare how the momentum transfer affects these distributions. To this end, we plot the unpolarized (f 1 (x)), helicity (g 1 (x)) and transversity (h 1 (x)) PDFs for P 3 = 1.67 GeV in the top panel of Fig. 11. All distributions are of similar magnitude and shape. f 1 (x) and h 1 (x) decay faster to zero as x increases, while g 1 (x) has a comparatively slower decay. The slope of f 1 (x) and  Fig. 11). We note that this plot corresponds to the maximum available momentum, P 3 = 1.25 GeV. However, we observed convergence with momentum for the leading GPDs and their PDFs, so comparison with the upper and lower left plots is acceptable. As in the previous two plots, the unpolarized H(x, 1/3, −1.02 GeV 2 ) is lower than the other two distributions. In the ERBL region, the transversity is slightly higher than H and H. In the DGLAP region, H and H T are compatible, while H is a bit lower. We want to emphasize again that these observations are only qualitative.

VI. MOMENTS OF GPDS
The Mellin moments of GPDs are defined via These are interesting in their own right, as they are related to form factors and towers of generalized form factors.
Recently, there has been exploration of the Mellin moments of quasi-GPDs, and their relation to the moments of GPDs. Of particular relevance is the work of Ref. [108], which derives relations for the Mellin moments of the transversity quasi-GPDs and GPDs using model-independent arguments. Also, a numerical analysis is presented using the diquark spectator model. The following model-independent relations are given for the n = 0 Mellin moments for both the GPDs and the quasi-GPDs 4 , As can be seen, the lowest moments of GPDs are independent of ξ, and the lowest moments of quasi-GPDs are, in addition, P 3 -independent. The form factors A T 10 , B T 10 , A T 10 are extracted from the matrix element of the local tensor operator as defined in Ref. [7]. Note that the lowest moment of E T is zero due to time reversal symmetry [6].
The corresponding relation for the n = 1 Mellin moments of the transversity GPDs are related to the generalized form factor of the one-derivative tensor operator, that is For the two lowest moments, that is, n = 0, 1 (Eqs. (38) -(45)), a ξ-dependence appears only in the n = 1 moment of E T , as it is the only ξ-odd GPD. H T , E T and H T are even functions of ξ.
Here, we calculate the moments of the transversity GPDs as a consistency check of our results. Our goal is not to provide numerical results for the form factors and generalized form factors, as the calculation of these quantities are still at an exploratory stage, but to perform a number of checks for the Mellin moments using our results: For the H T -GPD, we find the following values using our lattice data For the E T -GPD, we have and for the H T -GPD, we find Similar to the Mellin moments of GPDs, the form factors do not depend on the momentum boost of the proton. Since Eqs. (55) -(57) are extracted directly from the matrix elements, we can also provide estimates for B T 10 and A T 10 at P 3 = 0.83 GeV.
1. For the quasi-GPDs at t = −0.69 GeV 2 and ξ = 0, we have three momenta for H T q (Eq. (46)) and two momenta for E T q and H T q (Eq. (49) and Eq. (52)). The two lowest P 3 for H T are in agreement, and in slight tension with P 3 = 1.67 GeV. The values for E T q between P 3 = 1.25 GeV and P 3 = 1.67 GeV are consistent within the uncertainties. It should be mentioned, however, that the uncertainties are much larger than for H T q . Similar conclusions to the ones for E T q are also valid for H T q .
We observe that the agreement of both the n = 0, 1 moments of the GPDs for different P 3 values is better than for the quasi-GPDs. This is an indication that the matching procedure removes the bulk of the P 3 -dependence.
2. The n = 0 moments of quasi-GPDs for a given value of P 3 are fully compatible with the results of the n = 0 moment of the corresponding GPDs for the same value of P 3 .

3.
For all the n = 0 moments that we present here, we find that the values at t = −1.02 GeV 2 are lower than those at t = −0.69 GeV 2 , as expected. For n = 1, we observe a flatter behavior with increase of −t in H T . This is similar to the t-dependence of past calculations, for example, of A T 20 using one-derivative operators [124,125]. For E T and H T , the signal decays to zero at t = −1.02 GeV 2 .

4.
Another outcome of the numerical analysis is the fact that the n = 1 moment of a given GPD is suppressed compared to n = 0. This is expected, as the higher moments have support at higher values of x, where the GPDs decay.

5.
Finally, we compare the n = 0 moments (Eqs. (38) - (40)) with the value of the matrix elements at z = 0 (Eqs. (55) - (57)). We find that these are in excellent agreement with the values obtained from the integrals, for both t = −0.69 GeV 2 and t = −1.02 GeV 2 . Regarding the case of E T , the results are consistent with zero.
The above conclusions are highly nontrivial 5 , as the extraction of the Mellin moments from the final GPDs includes the reconstruction of the x-dependence and the matching. Therefore, these results serve as very important cross-checks of the validity of our results.

VII. SUMMARY
In this paper, we present the first lattice QCD calculation of transversity GPDs for the proton, employing the quasidistribution approach. GPDs are defined in the Breit frame, which we employ in this work. We use kinematic setups for both zero and nonzero skewness. In particular, we present results for {ξ, t} = {0, −0.69 GeV 2 } using momentum boosts P 3 = 0.83, 1.25, 1.67 GeV. For nonzero skewness, we have {|ξ|, t} = {1/3, −1.02 GeV 2 } for P 3 = 1.25 GeV. The matrix elements are renormalized in position space using a variation of the RI-MOM scheme, the so-called minimal projector. The choice of the projector is such that it isolates the tree-level contributions from the vertex functions of the operator. This is necessary, as the available matching formulas have been developed for that scheme [117]. To compute the x-dependence of the GPDs, we apply the Backus-Gilbert method to obtain the quasi-GPDs. Finally, we apply the perturbative matching equations to extract the light-cone GPDs. In particular, the analytic equations of the matching relate the quasi-GPDs defined in the RI scheme at a scale µ 0 , to the physical GPDs in the MS scheme at 2 GeV.
We use a combination of operators, momentum source and sink, as well as parity projectors, so that we can disentangle the four transversity GPDs, H T , E T , H T and E T . For the latter, we find zero signal within uncertainties, as it is suppressed compared to the other GPDs. The P 3 -dependence, at fixed −t = 0.69 GeV 2 and ξ = 0, is investigated boosting the proton at P 3 = 0.83, 1.25 and 1.67 GeV. Our results in Fig. 6 and Fig. 8 show that momentum convergence in H T is observed at the two highest boosts. A much larger statistics is needed to fully establish such a conclusion for E T and H T , that suffer from large statistical errors. Nevertheless, there is qualitative agreement between our lattice results and the analysis of GPDs in the scalar diquark model of Ref. [108], where, for example, H T is negative, in agreement with our findings. At P 3 = 1.25 GeV, we also extract the GPDs at |ξ| = 1/3 and −t = 1.02 GeV 2 . At non-zero ξ, there is a non-trivial distinction between the ERBL (−ξ < x < +ξ) and DGLAP (−1 < x < −ξ, ξ < x < 1) regions, and we find that the t-dependence of the GPDs is more prominent in the ERBL region.
In addition to the individual GPDs, we extract the combination E T + 2 H T (see Fig. 10), both at zero and non-zero skewness and for the two t-values considered in this work. This quantity provides the transverse spin-flavor dipole moment in an unpolarized target, k T , through its lowest moment and in the forward limit (t = 0). At the present stage, we are in no position to estimate k T , because that would require the knowledge of E T and H T for multiple t-values to extract their values at t = 0 through fits. This is certainly a very interesting direction, that we will pursue in the future.
Our results for the transversity GPDs are combined with the unpolarized and helicity GPDs from Ref. [68], that were calculated on the same ensemble and for the same kinematic setup. We compare the three types of PDFs, and the effect of introducing momentum transfer and nonzero skewness (see Fig. 11). As expected, the GPDs, H, H and H T , are suppressed compared to their PDF counterparts.
Another aspect of our analysis is the calculation of the two lowest Mellin moments for the GPDs. We also extract the lowest moment of quasi-GPDs using the relations of Ref. [108]. This is an important part of this work, leading to a number of conclusions that are consistent with the expected relations. In a nutshell, we find that the n = 0 Mellin moments of quasi-GPDs do not depend on P 3 , even though the quasi-GPDs have an explicit dependence on P 3 . In addition, the moments of GPDs obtained at different momenta are consistent. The expectation that the n = 0 Mellin moments of GPDs and quasi-GPDs are the same, is confirmed by our results numerically. Also, the Mellin moments have the expected t-dependence, that is, they decrease as −t increases. Another conclusion is that going from n = 0 to n = 1 results in decreasing the values for the moment. Last, but not least, the n = 0 moments are fully consistent with their extraction from the matrix elements at z = 0. These conclusions hold for all transversity GPDs, except E T , which is consistent with zero within our precision.
The calculation presented here is the first of a series of studies aiming at the calculation of GPDs on several ensembles, in order to quantify systematic uncertainties such as pion mass dependence and discretization effects. Having results from larger-volume ensembles will allow us to obtain the GPDs for several values of t and fit the t-dependence. As previously mentioned, this is important for obtaining the forward limit for the GPDs that drop out of the matrix element at t = 0. In this way, lattice QCD can provide a robust way of probing the three-dimensional structure of the nucleon and complement the rich experimental programs aiming at unraveling this structure.