Pion generalized parton distribution from lattice QCD

We present the first direct lattice calculation of the valence-quark generalized parton distribution (GPD) of the pion using the large-momentum effective theory (LaMET) approach. We focus on the zero-skewness limit, where the GPD has a probability-density interpretation in the longitudinal Bjorken $x$ and the transverse impact-parameter distributions. Our calculation is done using clover valence fermions on an ensemble of gauge configurations with $2+1+1$ flavors (degenerate up/down, strange and charm) of highly improved staggered quarks (HISQ) with lattice spacing $a = 0.12$~fm, box size $L \approx 3$~fm and pion mass $m_\pi \approx 310$~MeV. The parton distribution function and the form factor are reproduced as special limits of the GPD as expected. Due to the large errors, this exploratory study does not show a clear preference among different model assumptions about the kinematic dependence of the GPD. To discern these assumptions, future studies using higher-statistics data will be crucial.


I. INTRODUCTION
In the past few decades, extensive studies of parton distribution functions (PDFs) have provided us with detailed knowledge of the longitudinal momentum distribution of quarks and gluons, and therefore a onedimensional picture of hadrons. To map out the multidimensional partonic structure of hadrons, which is an important goal for experiments carried out at the existing facilities in DESY, JLab, BNL, CERN or the planned Electron-Ion Collider, we need to study quantities exhibiting the transverse structure of hadrons. One such quantity that has attracted a lot of interest in the past few years are the generalized parton distributions (GPDs) [1,2] (see also [3]).
The GPDs unify seemingly different physical quantities, such as the PDFs and hadron form factors, into the same framework. They offer a description of the correlations between the transverse position and longitudinal momentum of quarks and gluons inside the nucleon, thereby giving access to quark and gluon orbital angular momentum contributions to the nucleon spin [1]. Experimentally, the GPDs can be accessed through hard exclusive processes like deeply virtual Compton scattering or meson production. Useful constraints on the forms of the nucleon GPDs have been obtained from measurements of such processes at DESY [4][5][6] and JLab [7][8][9][10]. However, as the GPDs usually contribute to experimental observables through convolutions and they have more complicated kinematic dependence than the PDFs, extracting the GPDs from these experimental measurements is in general rather difficult. Therefore, inputs from theory are In recent years, a new theory framework has been developed that allows for lattice calculations of the xdependence, instead of the moments, of parton quantities [17,18]. This theory is now known as the largemomentum effective theory (LaMET). In this approach, a parton observable such as the PDFs or the GPDs can be accessed from lattice QCD in the following manner: 1) Constructing an appropriate static-operator matrix element (a quasi-observable) that approaches the parton observable in the large-momentum limit of the external hadron. The quasi-observable constructed this way is usually hadron-momentum-dependent but timeindependent, and, therefore, can be readily computed on the lattice. 2) Calculating the quasi-observable on the lattice. 3) Converting it to the parton observable through a factorization formula accurate up to power corrections that are suppressed by the hadron momentum. The existence of such a factorization is ensured by construction; for a proof, see Refs. [19][20][21].
As for the GPDs, the factorization of the isovector quark quasi-GPDs has been proven to leading-power accuracy using the operator product expansion [21], and the corresponding hard matching function was also computed both in a cutoff scheme [24,25] and in a regularization-independent momentum-subtraction (RI/MOM) scheme [21] (for studies of quasi-GPDs in diquark models see e.g. [62,84]). This allows us to perform exploratory studies on the quark GPDs once we have lattice simulations of the corresponding quasi-GPD matrix elements.
In this paper, we carry out the first direct lattice calculation of the valence-quark GPD of the pion using the LaMET approach. As a first step, we focus on the zeroskewness limit, that is, where the momentum transfer between the initial and final states is purely transverse. In this limit, the quark GPD is related to the impactparameter distribution of quarks that has a probabilitydensity interpretation [85] (see also Ref. [86]). Our calculation is done using clover valence fermions on an ensemble of gauge configurations with 2 + 1 + 1 flavors (degenerate up/down, strange and charm) of highly improved staggered quarks (HISQ) [87] generated by the MILC Collaboration [88] with lattice spacing a = 0.12 fm, box size L ≈ 3 fm and pion mass m π ≈ 310 MeV.

II. FROM QUASI-GPD TO GPD IN THE PION
The unpolarized quark GPD in the pion is defined on a lightcone as whereq, q denote the quark fields, with P µ = (P 0 , 0, 0, P z ). µ is the renormalization scale in the MS scheme and is the gauge link that ensures gauge invariance of the quark bilinear operator. Eq. (1) is an off-forward matrix element where the momenta for the initial and final states are different. In the forward (∆ µ → 0) limit, it reduces to the PDF. An appropriate quark quasi-GPD that can be computed on the lattice is given bỹ where we have chosen the Dirac matrix as γ t , since it has the advantage of avoiding mixing with the scalar quark operator [30,44] when a non-chiral lattice fermion is used. This choice will be used throughout this paper. The skewness parameter ξ in Eq. (4) is defined as ξ = −∆ z /(2P z ), which differs from the lightcone definition in Eq. (1) by power-suppressed contributions of O(m 2 π /P 2 z ). We have ignored this difference and denoted it with the same label as the skewness in the GPD.μ denotes the renormalization scale in an appropriate renormalization scheme for the quasi-GPD. In the present paper, we will focus on the u − d combination at zero skewness ξ = 0, where the former avoids contributions from disconnected diagrams as well as the mixing with gluon GPDs, while the latter simplifies the kinematic dependence of the quark GPD and is also related to the impactparameter distribution of quarks that has a probabilitydensity interpretation [85,89].
The bare pion matrix element on the right-hand side of Eq. (4) can be calculated on the lattice. In Refs. [31,40,41], it has been shown that the quark bilinear operator definingh is multiplicatively renormalized, and the renormalization factor can be calculated nonperturbatively on the lattice. In our previous study of the pion PDF [73], we chose to calculate the renormalization factor in the RI/MOM scheme, where the counterterm is determined by requiring that it cancels all the loop contributions for the matrix element in an off-shell external quark state at a specific momentum [29,42]. In other words, the renormalization factor iñ is fixed by where Λ γ t is the amputated Green function of the forward quark bilinear operator in Eq. (5) in an off-shell quark state with momentum p. P is a projection operator that defines the RI/MOM renormalization factor, µ R , p R z are renormalization scales introduced in the RI/MOM scheme. After renormalization, all singular dependence on a has been removed, andh R has a well defined continuum limit. We have suppressed the residual a dependence inh R . The Z factor defined in Eq. (7) coincides with that of the quark quasi-PDF. Since the UV divergence of the above hadron matrix element depends only on the operator defining it and not on the external state, the same renormalization factor can be used to renormalize the quark quasi-GPD matrix element. After renormalization,h R (z, P z , ξ, t, p R z , µ R ) can be converted toH π q through a Fourier transform e ixP z zh R (z, P z , ξ, t, p R z , µ R ), (8) which can then be factorized into the normal GPD in the MS scheme convoluted with a perturbative hard matching kernel, up to power corrections that are suppressed by the pion momentum [21] H π u−d,R (x, ξ, t, P z , p R z , µ R ) where µ is the renormalization scale of the GPD. The matching kernel C has been worked out at one-loop in Ref. [21]. At zero skewness ξ = 0, C is the same as the matching kernel for the PDF that is documented in Refs. [60,72]. Ideally, the continuum limit ofH u−d,R should be taken before applying the matching so that lattice artifacts can be removed and the rotational symmetry recovered. However, only a single lattice spacing is used in the present work. The continuum limit can be explored in the future once we have more data at different lattice spacings. For the power corrections, the meson-mass correction associated with the choice of Dirac matrix γ t is identical to that of the helicity distribution worked out in Ref. [36] with the replacement m 2 π →m 2 π = m 2 π − t/4. The O(Λ 2 QCD /P 2 z ) correction is parametrically about the same size as the O(m 2 π /P 2 z ) correction (except for very small/large x), which is negligible compared with other sources of errors.
ForH π + u−d , the matrix elementh is purely real in the isospin symmetric limit adopted in this work, since its imaginary part is related to the inverse Fourier transform , which is 0 after equating u(d) withd(ū) using isospin symmetry. Analogously, the real part of the matrix element is related to the valence quark quasi-GPD asH π The above analysis applies not only to the pion quasi-GPD but also to the pion GPD. In the following, we will present our skewless isovector combination of valence quark GPD for the charged pion π + as where the dependence on the renormalization scale µ is suppressed.

III. LATTICE CALCULATION SETUP
In this work, we use a single ensemble of gauge configurations with 2 + 1 + 1 flavors (degenerate up/down, strange and charm) of highly improved staggered quarks (HISQ) [87] generated by the MILC Collaboration [88] with lattice spacing a = 0.12 fm, pion mass m π ≈ 310 MeV, and box size L ≈ 3 fm (M π L ≈ 4.5). Our calculation is done using clover valence fermions on top of one-step hypercubic(HYP)-smeared [90] gauge links, with the clover parameters tuned to recover the lowest pion mass of the staggered quarks in the sea [91][92][93][94]. Then we calculate the time-independent, nonlocal (in space, chosen to be in the z direction) correlators of a pion with a finite-P z boost h lat (z, P z , t, a) where U z is a discrete gauge link in the z direction, P = {0, 0, P z } is the momentum of the pion, Γ = γ t and ∆ is the momentum transfer between initial and final pion. In this work, we only deal with the zeroskewness limit ξ = 0, where the matching coincides with that for the PDF. We use 3 boosted pion momenta, P z = {0, 0, n 2π L } with n ∈ {2, 3, 4}, which correspond to 0.86, 1.32 and 1.74 GeV, respectively. The initial and final pion momenta are obtained from {0, 0, n 2π L } ∓ ∆/2, where ∆ = {n x , n y , 0} 2π L with n 2 x + n 2 y ≤ 5. We carefully tune the Gaussian smearing parameter to best control the excited state and use four source-sink separations, 0.72, 0.84, 0.96 and 1.08 fm to help us remove the excitedstate contamination from our three-point correlators fits to extract the pion matrix element. We then momentumaverage the spatial symmetry around the x-y plane for each momentum transfer values of t.
To make sure that we have full control of excitedstate contamination, we analyze our data using different source-sink separations and using different levels of excited-state treatment. First, we use the "two-sim" analysis described in Ref. [94] to obtain the groundstate pion matrix elements using all 4 source-sink separations. The "two-sim" analysis only takes account of the leading excited-state contamination coming from the excited-and ground-state mixing. This is the same level as other commonly used methods, such as the "summation" method [95]. A second extraction uses only the largest two separations; if there is a significant excitedstate contamination at the smaller source-sink separation, 0.72 and 0.84 fm, we should see inconsistency in the ground-state matrix element between this analysis and earlier ones. Finally, we use the "two-twoRR" analysis (see Ref. [94] for details), which includes an additional matrix element related to excited states. Given the same input source-sink separations for the three-point correlators, the extracted ground-state matrix elements should be noisier, since more fit parameters are used. 1 All the above analyses generate consistent groundstate nucleon matrix elements. In this work, we will use matrix elements from the "two-twoRR" analysis in our PDF analysis. Fig. 1 shows the real part of the bare and RI/MOM renormalized matrix elements for P z = 8π/L ≈ 1.74 GeV, and t = {0, −2, −5} × (2π/L) 2 . The renormalization scales in the RI/MOM parameters have been chosen as µ R = 3.7 GeV, p R z = 6 × 2π/L. The error bars correspond to statistical errors. For t = 0, it is much narrower at small z than those for t = 0, due to charge conservation. The matrix elements at z = 0, which are not changed by the renormalization, are the values of the pion isovector formfactor. It is decreasing in |t| as expected.

IV. NUMERICAL RESULTS AND DISCUSSIONS
Now we present our numerical results for the valence quark GPD in the pion. As mentioned previously, the bare quasi-GPD matrix element calculated on the lattice can be renormalized by the RI/MOM renormalization factors for the quasi-PDF matrix element, which have been computed in Ref. [60]. The momentum distribution is then given bỹ where only the z and x dependence is shown for simplicity. Following our earlier work [45] on the nucleon PDF, we also apply the "derivative" method [45] to improve the truncation error in the Fourier transform in Eq. (12).
We then apply the one-loop matching and meson-mass corrections, where the latter turn out to be numerically negligible.
In Fig. 2, we show the results of the valence-quark distribution H π v (x, t, µ) for different values of t with the renormalization scale µ = 4 GeV. We have inverted the factorization formula in Eq. (9) by perturbatively expanding the matching kernel C to O(α s ). Also, the meson-mass power corrections have been removed. For the RI/MOM renormalization of the quark quasi-GPD, we have chosen µ R = 3.7 GeV, p R z = 6 × 2π/L. The error bands in Fig. 2 include statistical as well as systematic errors of p R z dependence by varying it between 4 × 2π/L and 8×2π/L. The curve with t = 0 is consistent with our previous pion PDF result [73] within errors. The impact parameter-space distribution can in principle be obtained by Fourier transforming the t dependence to the impact parameter b ⊥ dependence. However, we have only the 1 The detailed procedure for treating excited-state systematics can be found in Ref. [94] for the nucleon-charge case.
FIG. 1. The real part of the bare and renormalized matrix elements,hB(z, P z , t) andhR(z, P z , t), for the zero-skewness isovector valence quark GPD of the pion. The averaged pion momentum in the z-direction is P z = 4 × 2π/L and the momentum transfer squared is t = {0, −2, −5} × (2π/L) 2 from top to bottom, respectively.
results for very few values of t in this work. For the kinematic dependence of H π v (x, t, µ), a naive functional form is where q π v is the pion valence-quark PDF satisfying dx q π v (x, µ) = 1, and F π u−d (t) is the isovector form factor of π + with the normalization useful parametrization is which has been used to parametrize the unpolarized zeroskewness GPD of the nucleon and in discussing the fit to experimental data for the nucleon form factors [96].
With current uncertainties, we are unable to conclude which of the above parametrizations is favored. Although the first naive form might not be favored by the study of the GPD asymptotic behavior at x → 1 [85,96], it cannot be falsified in this exploratory study. Even if we remove all the uncertainties from the renormalization and matching and test the factorization at the bare matrix element level, the factorized form for the bare matrix element h lat (z, P z , t, a) = F π u−d (t)h lat (z, P z , t = 0, a) cannot be falsified. Future studies using higher-statistics data will be crucial to justify or falsify model assumptions such as Eqs. (13) and (14) about the kinematic dependence of the GPD.

V. SUMMARY
We have presented the first direct lattice calculation of the valence-quark generalized parton distribution of the pion using the LaMET approach. We have focused on the zero-skewness limit, where the GPD has a probability-density interpretation in the longitudinal Bjorken x and the transverse impact-parameter distributions. Our calculation is done using clover valence fermions on an ensemble of gauge configurations with 2+1+1 flavors (degenerate up/down, strange and charm) of highly improved staggered quarks (HISQ) with lattice spacing a = 0.12 fm, box size L ≈ 3 fm and pion mass m π ≈ 310 MeV. The parton distribution function and form factor are reproduced as special limits of the GPD as expected. Future studies using higher-statistics data will be crucial to justify or falsify different model assumptions about the kinematic dependence of the GPD.