Nucleon Tomography and Generalized Parton Distribution at Physical Pion Mass from Lattice QCD

We present the first lattice calculation of the nucleon isovector unpolarized generalized parton distribution (GPD) at the physical pion mass using a lattice ensemble with 2+1+1 flavors of highly improved staggered quarks (HISQ) generated by MILC Collaboration, with lattice spacing $a\approx 0.09$~fm and volume $64^3\times 96$. We use momentum-smeared sources to improve the signal at nucleon boost momentum $P_z \approx 2.2$ GeV, and report results at 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 distributions $H(x,Q^2)$ and $E(x,Q^2)$ at $\xi=0$ are 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.

We present the first lattice calculation of the nucleon isovector unpolarized generalized parton distribution (GPD) at the physical pion mass using a lattice ensemble with 2+1+1 flavors of highly improved staggered quarks (HISQ) generated by MILC Collaboration, with lattice spacing a ≈ 0.09 fm and volume 64 3 × 96. We use momentum-smeared sources to improve the signal at nucleon boost momentum Pz ≈ 2.2 GeV, and report results at nonzero momentum transfers in [0.2, 1.0] GeV 2 . Nonperturbative renormalization in RI/MOM scheme is used to obtain the quasi-distribution before matching to the lightcone GPDs. The three-dimensional distributions H(x, Q 2 ) and E(x, Q 2 ) at ξ = 0 are presented, along with the three-dimensional nucleon tomography and impact-parameterdependent distribution for selected Bjorken x at µ = 3 GeV in MS scheme. Nucleons (that is, protons and neutrons) are the building blocks of all ordinary matter, and the study of nucleon structure is a central goal of many worldwide experimental efforts. Gluons and quarks are the underlying degrees of freedom that explain the properties of nucleons, and fully understanding how they contribute to the properties of nucleons (such as their mass or spin structure) helps to decode the Standard Model. In quantum chromodynamics (QCD), gluons strongly interact with themselves and with quarks, binding both nucleons and nuclei. However, due to their confinement within these bound states, we cannot single out individual constituents to study them. More than half a century since the discovery of nucleon structure, our understanding has improved greatly; however, there is still a long way to go in unveiling the nucleon's detailed structure, which is characterized by functions such as the generalized parton distributions (GPDs) [1][2][3]. GPDs can be viewed as a hybrid of parton distributions (PDFs), form factors and distribution amplitudes. They play an important role in providing a three-dimensional spatial picture of the nucleon [4] and in revealing its spin structure [2]. Experimentally, GPDs can be accessed in exclusive processes such as deeply virtual Compton scattering [5] or meson production [6]. Experimental collaborations and facilities worldwide have been devoted to searching for these last unknowns of the nucleon, including HERMES at DESY, COMPASS at CERN, GSI in Europe, BELLE and JPAC in Japan, Halls A, B and C at Jefferson Laboratory, and PHENIX and STAR at RHIC at Brookhaven National Laboratory in the US. There are also plans for future facilities: a US electron-ion collider (EIC) [7] at Brookhaven National Laboratory, an EIC in China (EicC) [8,9], and the Large Hadron-Electron Collider (LHeC) in Europe [10,11].
Although interest in GPDs has grown enormously, we still need fresh theoretical and phenomenological ideas, including reliable model-independent techniques.
Most QCD models have issues associated with threedimensional structure that are not yet fully understood, so a reliable framework for extracting three-dimensional parton distributions and fragmentation functions from experimental observables does not yet exist. Theoretically, there are factorization issues in hadron production from hadronic reactions, and theoretical efforts are striving to answer key questions that lie along the path to a precise mapping of three-dimensional nucleon structure from experiment. It has become common understanding that we need to develop a program in both theory and experiment that will allow an accurate flavor decomposition of the nucleon GPDs, including flavor differences in the quark structure, the gluon structure and the nucleon sea-quark GPDs. Most current theoretical issues are associated with nonperturbative features of QCD, that is, where the strong coupling is too large for analytic perturbative methods to be valid. Using a nonperturbative theoretical method that starts from the quark and gluon degrees of freedom, lattice QCD (LQCD), allows us to compute these properties on supercomputers.
Probing hadron structure with lattice QCD was for many years limited to the first few moments, due to complications arising from the breaking of rotational symmetry by the discretized Euclidean spacetime. The breakthrough for LQCD came in 2013, when a technique was proposed to connect quantities calculable on the lattice to those on the lightcone. Large-momentum effective theory (LaMET), also known as the "quasi-PDF method" [12][13][14], allows us to calculate the full Bjorken-x dependence of distributions for the first time. Much progress has been made since the first LaMET paper . Most work has been done using only one lattice ensemble, but there has been some progress in determining the size of lattice systematic uncertainties. For example, finite-volume systematics were studied in Refs. [31,86]. Machinelearning algorithms have been applied to the inverse arXiv:2008.12474v2 [hep-ph] 23 Jul 2021 problem [84,105] and to making predictions for larger boost momentum and larger Wilson-displacement [106]. On the lattice discretization errors, a N f = 2 + 1 + 1 superfine (a ≈ 0.042 fm) lattice at 310-MeV pion mass was used to study nucleon PDFs in Ref. [91], and results using multiple lattice spacings were reported in Refs. [89,92,105]. The first attempt to obtain strange and charm distributions of the nucleon was recently reported [90]. However, beyond one-dimensional hadron structure, there is little work available. Last spring, the first lattice study of GPDs was made for pions [32]. Dur-ing the completion of this work, ETM Collaboration reported their findings on both unpolarized and polarized nucleon GPDs with largest boost momentum 1.67 GeV at pion mass M π ≈ 260 MeV [107,108]. In this work, we present the first lattice-QCD calculation of the nucleon GPD at the physical pion mass using the LaMET method and study the three-dimensional structure of the unpolarized nucleon GPDs.
The unpolarized GPDs H(x, ξ, t) and E(x, ξ, t) are defined in terms of the matrix elements where L(−z/2, z/2) is the gauge link along the lightcone and In the limit ξ, t → 0, H reduces to the usual unpolarized parton distributions while the information encoded in E cannot be accessed, since they are multiplied by the fourmomentum transfer ∆. Only in exclusive processes with a nonzero momentum transfer can E be probed. The one-loop matching [35,144] for the GPD H and E turns out to be similar to that for the parton distribution.
In this work, we focus on the nucleon isovector unpolarized GPDs and their quasi-GPD counterparts defined in terms of spacelike correlations calculated in Breit frame. We use clover valence fermions on an ensemble with lattice spacing a ≈ 0.09 fm, spatial (temporal) extent around 5.8 (8.6) fm, and with the physical pion mass M π ≈ 135 MeV and N f = 2+1+1 (degenerate up/down, strange and charm) flavors of highly improved staggered dynamical quarks (HISQ) [109] generated by MILC Collaboration [110]. The gauge links are one-step hypercubic (HYP) smeared [111] to suppress discretization effects. The clover parameters are tuned to recover the lowest sea pion mass of the HISQ quarks. The "mixed-action" approach is commonly used, and there has been promising agreement between the calculated lattice nucleon charges, moments and form factors and the experimental data when applicable [112][113][114][115][116][117][118][119][120][121][122][123][124]. Gaussian momentum smearing [125] is used on the quark field to improve the overlap with ground-state nucleons of the desired boost momentum, allowing us to reach higher boost momentum for the nucleon states. We calculate the matrix elements of the form . We also use high-statistics measurements, 501,760 total over 1960 configurations, to drive down the increased statistical noise at high boost momenta, P z = | We solve a set of linear equations to obtain H and E (similar to form-factor extraction) with all |q| at fixed Q 2 . Technical details (such as renormalization) and more information on how the matrix elements are extracted can be found in the supplement and our previous work [24,27,29,126].
The nonperturbatively renormalized matrix elements are then Fourier transformed into quasi-GPDs through two different approaches. Following the recent work [24,27,29], we take the matrix elements z ∈ [−12, 12] and apply the simple but effective "derivative" method, R /x, to obtain the quasi-GPDs. Alternatively, we adopt the extrapolation formulation suggested by Ref. [127] 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 physical GPDs by applying the matching condition [24,27,41,67]. Examples of the GPDs at momentum transfer Q 2 ≈ 0.4 GeV 2 are shown in Fig. 1. Figure 1 compares the H and E GPDs at Q 2 ≈ 0.4 GeV 2 with the quasi-distribution and matched distribution using P z ≈ 2.2 GeV. The matching lowers the positive mid-x to large-x distribution, as expected; as one approaches lightcone limit, the probability of a parton carrying a larger fraction of its parent nucleon's momentum should become smaller. We also find that the 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, which contributes more significantly to the small-x distribution. By repeating a similar analysis for each available Q 2 in this calculation, we can construct the full three-dimensional shape of H and E as functions of x and Q 2 , as shown in Fig. 2. Due to the limited reliable zP z reach of the lattice calculation, we find that the small-x region is unreliable, due to lack of precision lattice data to constrain it. Thus, due to charge conservation, the antiquark (negative-x) distribution can also be sensitive to P z . It has been found in past work [18,24,27,29] that higher boost momenta are needed to improve the antiquark region. Therefore, for the rest of the work, we will mainly focus on the x > 0.05 region. For convenience, we will focus on showing the GPD results from derivative method, and use the Reggeinspired extrapolation to estimate the uncertainties in the small-x region by reconstructing GPD moments from our x-dependent GPD functions.
Since this is the first lattice calculation with full threedimensional x and Q 2 dependence of the H and E GPD functions, we would like to check how the new results using LaMET approach compare with the previous moment approaches to the generalized form factors. In the ξ → 0 limit, the H and E GPDs decrease monotonically as x or Q 2 increases. We take Mellin moments of the GPDs to compare with previous lattice calculations done using local matrix elements through the operator product expansion (OPE). Taking the x-moments of H and E [128,129]: where the generalized form factors (GFFs) A ni (Q 2 ), B ni (Q 2 ) and C ni (Q 2 ) in the ξ-expansion on the righthand side are real functions. When n = 1, we get the Dirac and Pauli electromagnetic form factors F 1 (Q 2 ) = A 10 (Q 2 ) and F 2 (Q 2 ) = B 10 (Q 2 ). To compare with other lattice results, we plot the Sachs electric and magnetic form factors using F 1,2 as G E (Q 2 ) = F 1 (Q 2 ) +  [18], while the pink band corresponds to the matched GPD using quasi-GPD from the extrapolation formulation suggested by Ref. [127]. We find both methods give reasonable agreement in the x-dependent behavior, except in the small-x region, which is dominated by the large-z matrix elements that rely on the extrapolation. Fig. 3 with selected results from near-physical pion masses. PACS has the largest volume among these calculations and is able to probe the smallest Q 2 . Overall, our results are not only consistent within errors with the earlier PNDME study using the same ensemble (but which used local operators) but are also in good agreement with other lattice collaborations. When n = 2, we obtained GFFs of A 20 (Q 2 ) and B 20 (Q 2 ) so that we can compare our moment results with past lattice calculations using the OPE, as shown in Fig. 3. We compare our moment results with those obtained from simulations at the physical point by ETMC using three ensembles [130] and the near-physical calculation of RQCD [131]. We note that even with the same OPE approach by the same collaboration, the two data sets for A 20 in the ETMC calculation exhibit some tension. This is an indication that the systematic uncertainties are more complicated for these GFFs. Given that the blue points correspond to finer lattice spacing, larger volume and larger m π L, we expect that the blue points have suppressed systematic uncer- tainties. Our moment result A 20 (Q 2 ) is in better agreement with those obtained using the OPE approach at small momentum transfer Q 2 , while B 20 (Q 2 ) is in better agreement with OPE approaches at large Q 2 . The comparison between the N f = 2 ETMC data and N f = 2 RQCD data reveals agreement for A 20 and B 20 . However, the RQCD data have a different slope than the ETMC data, which is attributed to the different analysis methods and systematic uncertainties. Both our results and ETMC's are done using a single ensemble; future studies to include other lattice artifacts, such as lattice-spacing dependence, are important to account for the difference in the results.
Note that the error bands in Fig. 3 include the systematics from the following: 1) The systematics due to the negative-and small-x regions of the current GPD extraction. We create pseudo-lattice data using input of CT18NNLO PDF [138] with the same lattice parameters (such as z and P z ) used in this calculation and apply the same analysis procedure described above. We 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); 2) We vary the maximum Wilson-line length z within 2 lattice units and take half the difference as an estimate of the systematic due to finite z; 3) We estimate 1/P z systematics due to higher-twist effects by comparing our LHPC17 2+1f PACS18 2+1f ETMC18 2+1+1f PNDME19 2+1+1f 0.09fm PNDME19 2+1+1f 0.06fm (top) The nucleon isovector electric and magnetic form factor results obtained from this work (labeled as "MSULat20 2+1+1") as functions of transferred momentum Q 2 , and comparison with other lattice works calculated near physical pion mass: N f = 2 ETMC18 [132] N f = 2 + 1 LHPC14 [133], LHPC17 [134], PACS18 [135]; N f = 2 + 1 + 1 ETMC18 [136], PNDME19 [137] with 2 lattice spacings of 0.06 and 0.09 fm. (bottom) The unpolarized nucleon isovector GFFs obtained from this work, compared with other lattice results calculated near physical pion mass as functions of transfer momentum Q 2 : ETMC19 [130], and RQCD19 [131]. Q 2 = 0 PDFs and to those in the previous works with 3 boost momenta [24,27]. The final errors are summed in quadrature to create the final error bands shown in Fig. 3.
The Fourier transform of the non-spin-flip GPD H(x, ξ = 0, Q 2 ) gives the impact-parameter-dependent distribution q(x, b) [139] q where b is the transverse distance from the center of momentum. Figure 4 shows the first results of impactparameter-dependent distribution from lattice QCD: a three-dimensional distribution as function of x and b, and two-dimensional distributions at x = 0.3, 0.5 and 0.7. The impact-parameter-dependent distribution describes the probability density for a parton with momentum fraction x at distance b in the transverse plane. Figure 4 shows that the probability decreases quickly as x increases. Using Eq. 4 and the H(x, ξ = 0, t = −q 2 ) obtained from the lattice calculation at the physical pion mass, we can take a snapshot of the nucleon in the transverse plane to perform x-dependent nucleon tomography using lattice QCD for the first time.
In this work, we compute the isovector nucleon unpolarized GPDs at physical pion mass using boost momentum 2.2 GeV with nonzero momentum transfers in [0.2, 1.0] GeV 2 . We are able to map out the first threedimensional GPD structures using lattice QCD in the special limit ξ = 0. There are residual lattice systematics are not yet included in the current calculation: In our past studies, we found the finite-volume effects to be negligible for isovector nucleon quasi-distributions calculated within the range M val π L ∈ {3.3, 5.5}. We anticipate such systematics should be small compared to the statistical errors. The lattice discretization has been studied by MSULat collaboration in Refs. [89,105] with multiple lattice spacings in the LaMET study of pion and kaon distribution amplitudes and PDFs; similarly, a comparison of nucleon isovector PDFs with 0.045 and 0.12 fm lattice spacing is shown in supplementary materials. There was mild lattice-spacing dependence for a majority of the Wilson-link displacements studied with similar largest boost momenta with same valence/sea lattice setup. EMTC also report LaMET isovector nucleon PDFs in Ref. [140] using twisted-mass fermion actions and reports different findings. Future work will investigate ensembles with smaller lattice spacing to reach even higher boost momentum (either directly or with the aid of machine learning [106]) so that we can push toward reliable determination of the smaller-x and antiquark regions.
We thank the MILC Collaboration for sharing the lattices used to perform this study. The LQCD calculations were performed using the Chroma software suite [141]. In the LaMET (or "quasi-PDF") approach, we calculate time-independent spatially displaced matrix elements that can be connected to the parton distributions. The operator choice is not unique at finite nucleon momentum; a convenient choice for leading-twist PDFs is to take the average of the hadron momentum (0, 0, P z ) and the quark-antiquark separation to be along the z direction where |N (P i,f ) denotes a nucleon state with momenta P µ i = (E i , −q x /2, −q y /2, P z ) and P µ f = (E f , q x /2, q y /2, P z ), ψ (ψ) are the (anti-)quark fields, n µ is a unit vector and W (λẑ, 0) is the Wilson link along the z direction. See Fig. 5 for an illustration. There are multiple choices of operator in this framework that will recover the same lightcone PDFs when the largemomentum limit is taken. For example, Γ can be γ z or γ t [33,57,81,142]; both will give the unpolarized PDFs in the infinite-momentum frame. In this work, we use Γ = γ t for the unpolarized GPDs calculations. Since the systematics of the LaMET methods, such as the higher- , decrease as the momentum increases, we use P z ≈ 2.2 GeV in this calculation. (A study of the P z convergence of the LaMET approach on the same ensembles can be found in Refs. [24,27,29].) This also allows us to reach smaller-x regions of the GPDs functions.
We use Gaussian momentum smearing [125] ψ(x) + α j U j (x)e ikêj ψ(x +ê j ), where k is the input momentum-smearing parameter, U j (x) are the gauge links in the j direction, and α is a tunable parameter as in traditional Gaussian smearing. We vary the α values in this study so that we can have multiple ways of removing the excited-state contributions from matrix elements later on. To better extract the boostedmomentum ground-state nucleon, we apply the variational method [143] to extract the principal correlators corresponding to pure energy eigenstates from a matrix of correlators. We use a 3×3 smeared-smeared two-point correlation matrix, which can be decomposed as with eigenvalues λ n (t, t r ) = e −(t−tr)En by solving the generalized eigensystem problem C ij where V is the matrix of eigenvectors and t r is a reference time slice. The resulting 3 eigenvalues (principal correlators) λ n (t, t r ) are then further analyzed to extract the energy levels E n . Since they have been projected onto pure eigenstates of the Hamiltonian, the principal correlator should be fit well by a single exponential and double checked for the consistency of the obtained energies. The leading contamination due to higher-lying states is another exponential having higher energy; we use a two-exponential fit to help remove this contamination. The overlap factors (A n ) between the interpolating operators and the n th state are derived from the eigenvectors obtained in the variational method.
To calculate the GPD matrix elements at nonzero momentum transfer, first calculate the matrix element χ N ( P f )|O µ |χ N ( P i ) , where χ N is the nucleon spin-1/2 interpolating field, abc [q a (x)Cγ 5 q b (x)]q c (x). O µ = ψγ µ W (z)ψ 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 also use high-statistics measure-ments, 501,760 total over 1960 configurations, to drive down the increased statistical noise at high boost momenta, P z =  Figure 6 shows an example of the ground-state matrix elements (shown as the gray band) from the variational analysis along with the two-state fitted results with P f = {1, 0, 10} 2π L and P i = {−1, 0, 10} 2π L with projection operator 1+γt 2 (1 + iγ 5 γ z ). The ground-state matrix elements can be extracted from v T 0 C ij 3pt (t)v 0 using the eigenvectors from the two-point variational-method analysis (shown in the left panel of Fig. 6), simultaneous two-state fitted results using source-sink separation of t sep = [8,12] lattice units, two-state fitted results as functions of t min sep . The consistency of these methods demonstrates that our extracted ground-state matrix elements are stably determined. For more details on the lattice study of nucleon matrix elements comparing the variational and two-state fit methods, we refer readers to Ref. [126], which contains a very detailed discussion.
The overdetermined system of linear equations (using multiple q and projection operators) allows for solution of h H (P z , Q 2 , z) and h E (P z , Q 2 , z), similar to vector form factor calculations with our chosen projection operator and momentum configuration: The ground-state matrix elements are proportional to with / p = iE N ( p )γ 4 + p · γ. We obtain a linear system of equations by using different projection operators P ρ , momenta P f and P i to solve for H and E in coordinate space, h {H,E} (P z , Q 2 , z). Selected Q 2 values of h {H,E} (P z , Q 2 , z) normalized by h {H,E} (P z , Q 2 = 0, z = 0) are shown in Fig. 7. The real matrix elements decrease quickly to zero due to the large boost momentum used in this calculation. This helps us to use smaller-displacement data to avoid large contributions from higher-twist effects in the larger-z region.
To obtain the quasi-GPDs, we first apply nonperturbative renormalization (NPR) in RI/MOM scheme to the bare matrix elements, using the NPR done in previous work using the same lattice ensembles [24,27,29], itself following the same strategy described in Refs. [22,41]. The RI/MOM renormalization constant Z is calculated nonperturbatively on the lattice by imposing the following momentum-subtraction condition on the matrix ele-ment: where On the lattice, ps|O t (z)|ps is calculated from the amputated Green function of O t with Euclidean external momentum. In this work, we use the same renormalization scales as used in the previous work [24,27,29]: µ R = 3.8 GeV, p R z = 2.2 GeV, µ MS = 3 GeV. The renormalization-scale dependence was studied in Ref. [67]. We also vary the values of P R z , and the results are shown in Fig. 8, focusing on the ξ = 0 GPDs, where the matching formula is the same as that in the PDFs, as discussed in Ref. [144]. We normalize all matrix elements by h H (P z , Q 2 = 0, z = 0), as in our previous PDF work [24,27,29]. Using matrixelement ratios reduces the lattice systematic error, since in the continuum limit h H (P z , Q 2 = 0, z = 0), the vector charge, goes to 1.
FIG. 6: An example of the ground-state matrix element determination using variational method using rotated two-and three-point correlators (left), two-state fit to multiple source-sink separation (middle), along with the ratio plots data and reconstructed fit to the ratio points. The rightmost plot shows how much the ground-state matrix elements change as we reduce the inputs to the fits; we see some fluctuations and increase of errors due to the reduction of the available data but overall steady determination of the ground-state matrix element. The example used here is from P f = {1, 0, 10} 2π L and Pi = {−1, 0, 10} 2π L with projection operator 1+γ t 2 (1 + iγ5γz). We also perform checks of the z max input in the Fourier transformation and lattice-spacing dependence. The ef-fects are documented in Fig. 8.