Nucleon Helicity Generalized Parton Distribution at Physical Pion Mass from Lattice QCD

The generalized parton distributions (GPDs) offer a window on three-dimensional imaging of the nucleon, providing understanding of how the fundamental properties of the nucleon, such as its mass and spin, arise from the underlying quark and gluon degrees of freedom. In this work, we present the first lattice calculation of the nucleon isovector helicity GPD at physical pion mass, using an $a \approx 0.09$ fm lattice ensemble with 2+1+1 flavors of highly improved staggered quarks generated by MILC Collaboration. We perform the GPD calculation in Breit frame using averaged nucleon boost momentum $P_z \approx 2.2$ GeV with nonzero momentum transfers in $[0.2,1.0]\text{ GeV}^2$. Nonperturbative renormalization in RI/MOM scheme is used to obtain the quasi-distribution before matching to the lightcone GPDs. The three-dimensional distribution $\tilde{H}(x,Q^2)$ is presented, along with the three-dimensional nucleon tomography and impact-parameter--dependent distribution for selected Bjorken $x$ at $\mu=3$ GeV in $\overline{\text{MS}}$ scheme.


INTRODUCTION
Three-dimensional images of the nucleon are critically important for understanding the origin of mass and orbital angular momentum as well as other properties of the nucleon. Generalized parton distributions (GPDs) use a single set of three-dimensional functions [1,2] to describe the spatial and momentum structures of the nucleon, including in its limits the parton distributions functions (PDFs) and the elastic form factors (FFs). Experimentally, GPDs can be accessed in exclusive processes such as deeply virtual Compton scattering or meson electroproduction [3]. Worldwide experimental collaborations and facilities have been devoted to searching for these last unknowns of the nucleon; these include HERMES at DESY, COMPASS at CERN, J-PARC in Japan, Halls A, B, and C at Jefferson Laboratory, and PHENIX and STAR at RHIC (Brookhaven National Laboratory) in the US. The pursuit of GPDs and their imaging has led the hadronic-physics community worldwide to plan future experiments: the Electron-Ion Collider [4], Electron-Ion Collider in China (EicC) [5,6], and the Large Hadron-Electron Collider (LHeC) in Europe [7,8]. However, obtaining high-quality images from experimental data remains decades away. Meanwhile, large-scale numerical lattice-QCD calculation can fill in some of the gaps where experiments cannot reach, such as zero-skewness limit. In this work, we focus on the longitudinally polarized GPD functionsH(x, ξ, t).

LATTICE MATRIX ELEMENTS
In this work, we focus on the nucleon isovector polarized GPDs and their quasi-GPD counterparts defined in terms of spacelike correlations calculated in Breit frame. We use an ensemble at physical pion mass with N f = 2+1+1 (degenerate up/down, strange and charm) flavors of highly improved staggered dynamical quarks (HISQ) [86] generated by MILC Collaboration [87]. The ensemble has lattice spacing a ≈ 0.09 fm, and fourdimensional volume of 64 3 × 96. We use one step of hypercubic (HYP) smearing [88] on the gauge links to suppress discretization effects. We use clover valence fermion action with the clover parameters tuned to recover the lowest sea pion mass of the HISQ quarks. The same "mixed-action" parameter choices have been used in other lattice calculations of nucleon charges, moments and form factors [30,[89][90][91][92][93][94][95][96][97][98][99][100], and there has been promising agreement between the calculated quantities and the experimental data when applicable. We use Gaussian momentum smearing [101] on the quark field to improve the overlap with the boosted-momentum ground-state nucleon.
We calculate matrix elements of the form with projection operators Γ = 1+γt 2 (1 + iγ 5 γ x,y,z ). An illustration of our Breit-frame setup can be found in Fig. 1. To calculate the GPD matrix elements at nonzero momentum transfer, we first calculate the matrix element ψ is the LaMET Wilson-line displacement operator with ψ being either an up or down quark field, and P {i,f } are the initial and final nucleon momenta. We integrate out the spatial dependence and project the baryonic spin, using projection operators P ρ = 1+γt 2 (1 + iγ 5 γ ρ ) with ρ ∈ {x, y, z}, leaving a time-dependent three-point correlator, C 3pt .
where f n,n (P f , P i , E n , E n , t, t i , t f ) contains kinematic factors involving the energies E n and overlap factors A n obtained in the two-point variational method, n and n are the indices of different energy states and Z O is the operator renormalization constant (which is determined nonperturbatively). We use high-statistics measurements, 501,760 total over 1960 configurations, to help with statistical noise at high boost momenta, P z = | , 0, 10}| with L = 5.63 fm. We vary the spatial momentum transfer q = P f − P i = 2π L {n x , n y , 0} with integer n x,y and n 2 = n 2 x + n 2 y ∈ {0, 4, 8, 16, 20} with four-momentum transfer squared Q 2 = −q µ q µ = {0, 0.19, 0.39, 0.77, 0.97} GeV 2 using periodic boundary condition with all |q| at fixed Q 2 , rotationally averaged. We set the quark momentum smearing parameter to {0, 0, 6} 2π L for n 2 = {0, 4, 8}, and {±1, ±1, 6} 2π L for n 2 = {16, 20}. We simultaneous fit the 501,760-measurement three-point correlators with source-sink separation of t sep = [8,12] lattice units to extract the ground-state matrix elements. The details of how we extract nucleon ground-state matrix elements can be found in our previous works [24,42,92].
The ground-state matrix elements are proportional to Because hẼ always couples with q z = (P f − P i ) z , for fixed operator Γ = γ z γ 5 and zero-skewness limit, the second term of the above equation will always be zero. We do not need to solve a linear system of equations to getH in coordinate space, hH (P z , Q 2 , z). We then normalize hH (P z , 0, 0) by Fig. 2. The real matrix elements decrease quickly to zero due to the large boost momentum used in this calculation. Using large boost momentum helps to reduce the contributions from higher-twist effects at O(Λ 2 QCD /P 2 z ).

RESULTS ON LATTICE HELICITY GPD
Following recent work [21,24,26], we nonperturbatively renormalize the matrix elements in RI/MOM scheme with µ R = 3.8 GeV. We then Fourier transform the renormalized matrix elements into quasi-GPDs through two approaches: 1) We take the matrix elements z ∈ [−12, 12] and apply the simple but effective "derivative" method,Q = i +zmax −zmax dz e ixPzzh R /x, to obtain the quasi-GPDs. 2) We adopt the extrapolation formulation suggested by Ref. [102] by fitting |z| ∈ {10, 15} using the formula c 1 (−izP z ) −d1 +c 2 e izPz (izP z ) −d2 , inspired by the Regge behavior, to extrapolate the matrix elements into the region beyond the lattice calculation and suppress Fourier-transformation artifacts. Then, both quasi-GPDs are matched to the lightcone GPDs by applying the matching condition [21,24,103]. Figure 3 shows the quasi-and lightcone distribution of theH GPD at momentum transfer Q 2 ≈ 0.4 GeV 2 using P z ≈ 2.2 GeV. We find that quasi-GPD using both derivative and Regge-inspired extrapolation agree in the mid-to large-x regions, but their difference grows as x approaches zero in both the quark and antiquark distribution. This is expected, since they differ mainly in the treatment of the large-z matrix elements in the quasi-GPD Fourier transformation outside the region with available lattice data (z > 15), which contributes more significantly to the small-x distribution. We also noticed some difference in the antiquark region (x < 0): this is also expected based on past work [15,21,24,26] that even higher boost momenta are needed to improve the antiquark region. We will take the difference between the two different ways of determining the quasi-GPD as a systematic reflected in the final uncertainties. We also see that the matching lowers the positive mid-x to large- x quasi-distribution, as expected. As one approaches the lightcone limit, the probability of a parton carrying a larger fraction of its parent nucleon's momentum should become smaller. Note that the matching from quasi-GPD to GPD has residual systematics at O into the systematic errors to estimate this effect, but future calculations including larger momenta will characterize this systematic better.
For convenience, we will focus on showing the GPD results from the derivative method, and use the Reggeinspired extrapolation to estimate small-x in reconstructing GPD moments from our x-dependent GPD functions. For the rest of the work, we will mainly focus on the x > 0.05 region. We repeat a similar analysis for each available Q 2 in this calculation, then use z-expansion up to 3 parameters to interpret the Q 2 dependence of the lightcone GPD functions; selected x results are shown in the top figure of Fig. 4. Through this fit, we can construct the full three-dimensional shape ofH as functions of x and Q 2 , as shown in the bottom of Fig. 4.
The zero-skewness limit of the GPD is related to the Mellin moments by taking the x-moments [104,105]: where the generalized form factors (GFFs)Ã ni (Q 2 ), and C ni (Q 2 ) in the ξ-expansion on the right-hand side are real functions. When n = 1 (n = 2), we get the axial form factors G A (Q 2 ) =Ã 10 (Q 2 ) (GFFsÃ 20 (Q 2 )). There have been a number of calculations of the Mellin moments of the GPDs on the lattice using local matrix elements through the operator product expansion (OPE), so we can make a comparison between LaMET results and moment methods. We show the resultsÃ 10 (Q 2 ) andÃ 20 (Q 2 ) obtained from this work (labeled as "MSULat21 2+1+1") as a function of momentum transfer Q 2 in Fig. 5, along with other lattice calculations near physical pion mass using traditional local-operator methods. We use z-expansion to interpret the Q 2 dependence of our calculation, and the results are shown as the green band in Fig. 5. The inner-band error indicates the statistical errors, while the outer one includes the systematic errors coming from sev- Nucleon isovectorH quasi-GPDs and lightcone GPDs at momentum transfer Q 2 = 0.39 GeV 2 . The orange and green bands are the quasi-GPD and lightcone GPDs from derivative method [15], while the pink band corresponds to the matched GPD using quasi-GPD from the extrapolation formulation suggested by Ref. [102]. We find both methods give reasonable agreement in the x-dependent behavior, except in the small-x and negative-x region, which is dominated by the large-z matrix elements that rely on the extrapolation. eral sources considered. First, the difference between the quasi-GPD methods, discussed previously. We also create pseudo-lattice data using the CT18NNLO PDF [106] with the same lattice z and P z parameters used in this calculation and take the upper limit of the reconstructed and original CT18 moments as an estimate of the systematics introduced by the analysis procedure (e.g. by Fourier truncation). These steps should account for the lack of sensitivity to the negative-and small-x regions in the current GPD extraction. We vary the maximum Wilson-line length z by 2 lattice units and take half the difference as an estimate of the systematic due to finite z used in obtaining the quasi-GPD. We also include an estimate of O(1/P z ) systematics due to potential highertwist effects by comparing our Q 2 = 0 PDFs to those in the previous works with 3 boost momenta [21,24]. The finite-volume effects studied in our past work using the HISQ ensembles [28] and an independent analytic study using ChPT [107] were found to be smaller than typical statistical errors; therefore, these are not considered in this work, given our M π L is close to 4. The final errors are summed in quadrature to create the final error bands shown in Fig. 5.
We found our axial form factors, obtained from theH GPD using LaMET method has slightly higher central values compared with other single-lattice-spacing lattice calculations, ranging from 2-flavor to 2+1+1-flavor ones using the axial-current approach with various choices of fermion and gauge action. However, overall, it is consistent with past lattice calculations, whose results are shown with statistical errors only. In the bottom plot of Fig. 5, we compare our moment results forÃ 20 (Q 2 ) with those obtained from simulations at the physical point by ETMC [108] and RQCD [109] using the OPE approach. We note that even with the same OPE approach oñ A 20 (Q 2 ), lattice results from a single lattice spacing with different actions do not agree; there have been indications that the systematic uncertainties are more complicated for these GFFs than those obtained for local currents, such as axial form factors. The OPE operators are only expected to give the same results after taking the continuum a → 0 limit, and current results only show a single lattice spacing. Nevertheless, we find our results to be consistent with ETMC's 2+1+1-flavor results; this is perhaps a coincidence, since both our results and ETMC's are done using a single ensemble. This first lattice calculation of the full three-dimensional x and Q 2 dependence of theH GPD functions using LaMET approach has a nice agreement with the previous moment approaches to the generalized form factors using using the OPE. The impact-parameter-dependent polarized quark distributions ∆q(x, b) [118] can be obtained fromH using where b is the transverse distance from the center of momentum. ∆q(x, b) has only been obtained from QCD models so far; this is the first lattice-QCD determination, shown in Fig. 6. We show the three-dimensional distribution as a function of x and b, and two-dimensional distributions at x = 0.3, 0.45 and 0.6 in Fig. 6. It describes the probability density for a polarized parton with momentum fraction x at distance b in the transverse plane. Compared with the q(x, b) obtained from the unpolarized GPD H function from earlier work done on the same ensemble, the center of the two-dimensional distribution q(x, b = 0) ≤ ∆q(x, b = 0) over most of the range of x. The probability densities decrease quickly as x and b increase.

CONCLUSION AND OUTLOOK
In this work, we compute the zero-skewness isovector nucleon polarizedH GPDs at physical pion mass using boost momentum 2.2 GeV with nonzero momentum transfers in [0.2, 1.0] GeV 2 . We use two different methods to construct the quasi-GPD before matching to its lightcone counterpart. We noted some differences using MSULat'21 2+1+1f  [94,[110][111][112][113][114][115][116], together with the form factor results obtained from this work (labeled as "MSULat21 2+1+1") by taking n = 1 in Eq. 4. All the existing lattice calculations are done at a single lattice spacing, except for PNDME19 [117] with 2 lattice spacings of 0.06 and 0.09 fm, which we distinguish in the legend. (bottom) The linearly polarized nucleon isovector GFFÃ20(Q 2 ) obtained from this work (labeled as "MSULat21 2+1+1f") by taking n = 2 in Eq. 4, compared with other lattice results calculated near physical pion mass as functions of transfer momentum Q 2 : 2+1+1f ETMC19 [108], 2f ETMC19 [108] (only the largervolume results are shown here), 2f RQCD19 [109]. these two methods, and they have been taken account into systematics in the final results. We also included the systematics from the choices of max Wilson-line displacement used in the LaMET operators and estimate the single boosted momentum dependence from previous helicity PDF studies. When taking the integral of ourH in Eq. 4, we are able to compare with previous calculations of the axial form factor G A andÃ 20 GFF with single ensembles near physical pion mass, since this is the first latticeH calculation at physical pion mass. We found consistent results with most of the past work within two sigma from various actions and other lattice parameters. This seems to verify that LaMET GPD cal- culation gives reasonable results for pursuing precision calculations in the future. We are able to map out the three-dimensional GPDH structures using lattice QCD data for the first time. Future work will investigate ensembles with smaller lattice spacing to reach even higher boost momentum so that we can push toward reliable determination of the smaller-x and antiquark regions, adopt the hybrid renormalization scheme to improve the longdistance nonperturbative behavior [102], and expand the current study into the ξ = 0 GPDs.

ACKNOWLEDGMENTS
We thank the MILC Collaboration for sharing the lattices used to perform this study. The LQCD calculations were performed using the Chroma software suite [119][120][121]. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 through ERCAP; facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, and supported in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research (iCER). The work of HL is partially supported by the US National Science Foundation under grant PHY 1653405 "CAREER: Constraining Parton Distribution Functions for New-Physics Searches" and by the Research Corporation for Science Advancement through the Cottrell Scholar Award "Unveiling the Three-Dimensional Structure of Nucleons".