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

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) Müller et al. (1994); Ji (1997a); Radyushkin (1996). 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 Burkardt (2000) and in revealing its spin structure Ji (1997a). Experimentally, GPDs can be accessed in exclusive processes such as deeply virtual Compton scattering Ji (1997b) or meson production Kriesten et al. (2020). 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) National Academies of Sciences and Medicine (2018) at Brookhaven National Laboratory, an EIC in China (EicC) Chen (2018); Chen et al. (2020), and the Large Hadron-Electron Collider (LHeC) in Europe Abelleira Fernandez et al. (2012); Agostini et al. (2020).

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 three-dimensional 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” Ji (2013, 2014); Ji et al. (2017), 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  Lin (2014a); Lin et al. (2015); Chen et al. (2016); Lin et al. (2018a); Alexandrou et al. (2015, 2017a, 2017b); Chen et al. (2018a); Alexandrou et al. (2018a); Chen et al. (2018b); Zhang et al. (2019a); Alexandrou et al. (2018b); Lin et al. (2018b); Fan et al. (2018); Liu et al. (2018a); Wang et al. (2019); Lin and Zhang (2019); Chen et al. (2019); Xiong et al. (2014); Ji and Zhang (2015); Ji et al. (2015a); Xiong and Zhang (2015); Ji et al. (2015b); Lin (2014b); Monahan (2018); Ji et al. (2018a); Stewart and Zhao (2018); Constantinou and Panagopoulos (2017); Green et al. (2018); Izubuchi et al. (2018); Xiong et al. (2017); Wang et al. (2018); Wang and Zhao (2018); Xu et al. (2018a); Zhang et al. (2017); Ishikawa et al. (2016); Chen et al. (2017); Ji et al. (2018b); Ishikawa et al. (2017, 2019); Li (2016); Monahan and Orginos (2017); Radyushkin (2017a); Rossi and Testa (2017); Carlson and Freid (2017); Ji et al. (2017); Briceño et al. (2018); Hobbs (2018); Jia et al. (2017); Xu et al. (2018b); Jia et al. (2018); Spanoudes and Panagopoulos (2018); Rossi and Testa (2018); Liu et al. (2018b); Ji et al. (2019a); Bhattacharya et al. (2019); Radyushkin (2019); Zhang et al. (2019b); Li et al. (2019); Braun et al. (2019); Detmold et al. (2019); Ebert et al. (2020a); Ji et al. (2019b); Bali et al. (2018a, b); Sufian et al. (2019); Bali et al. (2019a); Orginos et al. (2017); Karpie et al. (2018a, b, 2019); Joó et al. (2019a, b); Radyushkin (2018); Zhang et al. (2018); Lin et al. (2020); Zhang et al. (2020a); Fan et al. (2020a); Sufian et al. (2020); Shugert et al. (2020); Green et al. (2020); Chai et al. (2020); Shanahan et al. (2020); Braun et al. (2020); Bhattacharya et al. (2020); Ji et al. (2020a); Ebert et al. (2020b); Lin (2020a); Joó et al. (2020); Bhat et al. (2020); Fan et al. (2020b). 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. Joó et al. (2019b); Lin and Zhang (2019). Machine-learning algorithms have been applied to the inverse problem Karpie et al. (2019); Zhang et al. (2020b) and to making predictions for larger boost momentum and larger Wilson-displacement Zhang et al. (2020c). On the lattice discretization errors, a $N_{f}=2+1+1$ superfine ($a\approx 0.042$ fm) lattice at 310-MeV pion mass was used to study nucleon PDFs in Ref. Fan et al. (2020a), and results using multiple lattice spacings were reported in Refs. Zhang et al. (2020b); Lin et al. (2020); Sufian et al. (2020). The first attempt to obtain strange and charm distributions of the nucleon was recently reported Zhang et al. (2020a). However, beyond one-dimensional hadron structure, there is little work available. Last spring, the first lattice study of GPDs was made for pions Chen et al. (2019). During 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_{\pi}\approx 260$ MeV Alexandrou et al. (2019a, 2020a). 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,\xi,t)$ and $E(x,\xi,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 $\xi,t\to 0$, $H$ reduces to the usual unpolarized parton distributions while the information encoded in $E$ cannot be accessed, since they are multiplied by the four-momentum transfer $\Delta$. Only in exclusive processes with a nonzero momentum transfer can $E$ be probed. The one-loop matching Ji et al. (2015a); Liu et al. (2019) 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\approx 0.09$ fm, spatial (temporal) extent around 5.8 (8.6) fm, and with the physical pion mass $M_{\pi}\approx 135$ MeV and $N_{f}=2+1+1$ (degenerate up/down, strange and charm) flavors of highly improved staggered dynamical quarks (HISQ) Follana et al. (2007) generated by MILC Collaboration Bazavov et al. (2013). The gauge links are one-step hypercubic (HYP) smeared Hasenfratz and Knechtli (2001) 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 Mondal et al. (2020); Jang et al. (2020); Gupta et al. (2018a); Lin et al. (2018c); Gupta et al. (2018b, 2017); Yoon et al. (2017); Bhattacharya et al. (2016, 2015, 2014); Briceno et al. (2012); Bhattacharya et al. (2012); Lin (2020b). Gaussian momentum smearing Bali et al. (2016) 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 $\left\langle P_{f}\left|\bar{\psi}\left(-\frac{z}{2}\right)\Gamma L\left(-\frac{z}{2},\frac{z}{2}\right)\psi\left(\frac{z}{2}\right)\right|P_{i}\right\rangle$ with projection operators $\frac{1+\gamma_{t}}{2}(1+i\gamma_{5}\gamma_{x,y,z})$. 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}=|\frac{\vec{P_{i}}+\vec{P_{f}}}{2}|=|\frac{2\pi}{L}\{0,0,10\}a^{-1}|$, and vary spatial momentum transfer $\vec{q}=\vec{P_{f}}-\vec{P_{i}}=\frac{2\pi}{L}\{n_{x},n_{y},0\}a^{-1}$ with integer $n_{x,y}$ and $n_{x}^{2}+n_{y}^{2}\in\{0,4,8,16,20\}$ with four-momentum transfer squared $Q^{2}=-q_{\mu}q^{\mu}=\{{0,0.19,0.39,0.77,0.97}\}$ GeV². 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 Yoon et al. (2016); Chen et al. (2018b); Lin et al. (2018b); Liu et al. (2018a).

The nonperturbatively renormalized matrix elements are then Fourier transformed into quasi-GPDs through two different approaches. Following the recent work Chen et al. (2018b); Lin et al. (2018b); Liu et al. (2018a), we take the matrix elements $z\in[-12,12]$ and apply the simple but effective “derivative” method, $\tilde{Q}=i\int_{-z_{\text{max}}}^{+z_{\text{max}}}\!\!dz\,e^{ixP_{z}z}\tilde{h}^{\prime}_{R}/x$, to obtain the quasi-GPDs. Alternatively, we adopt the extrapolation formulation suggested by Ref. Ji et al. (2020b) by fitting $|z|\in\{10,15\}$ using the formula $c_{1}(-izP_{z})^{-d_{1}}+c_{2}e^{izP_{z}}(izP_{z})^{-d_{2}}$, 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 Stewart and Zhao (2018); Liu et al. (2018b); Chen et al. (2018b); Lin et al. (2018b). Examples of the GPDs at momentum transfer $Q^{2}\approx 0.4\text{ GeV}^{2}$ are shown in Fig. 1. Figure 1 compares the $H$ and $E$ GPDs at $Q^{2}\approx 0.4\text{ GeV}^{2}$ with the quasi-distribution and matched distribution using $P_{z}\approx 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 Lin et al. (2018a); Chen et al. (2018b); Lin et al. (2018b); Liu et al. (2018a) 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 Regge-inspired 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 three-dimensional $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 $\xi\to 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$ Ji (1998); Hagler (2010):

where the generalized form factors (GFFs) $A_{ni}(Q^{2})$, $B_{ni}(Q^{2})$ and $C_{ni}(Q^{2})$ in the $\xi$-expansion on the right-hand 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})+q^{2}F_{2}(Q^{2})/(2M_{N})^{2}$ and $G_{M}(Q^{2})=F_{1}(Q^{2})+F_{2}(Q^{2})$ in 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 Alexandrou et al. (2020b) and the near-physical calculation of RQCD Bali et al. (2019b). 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_{\pi}L$, we expect that the blue points have suppressed systematic uncertainties. 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 Hou et al. (2019) 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 $Q^{2}=0$ PDFs and to those in the previous works with 3 boost momenta Chen et al. (2018b); Lin et al. (2018b). 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,\xi=0,Q^{2})$ gives the impact-parameter–dependent distribution $\mathsf{q}(x,b)$ Burkardt (2003)

where $b$ is the transverse distance from the center of momentum. Figure 4 shows the first results of impact-parameter–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,\xi=0,t=-\mathbf{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]\text{ GeV}^{2}$. We are able to map out the first three-dimensional GPD structures using lattice QCD in the special limit $\xi=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_{\pi}^{\text{val}}L\in\{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. Zhang et al. (2020b); Lin et al. (2020) 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. Alexandrou et al. (2020c) 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 Zhang et al. (2020c)) so that we can push toward reliable determination of the smaller-$x$ and antiquark regions.