Isovector parton distribution functions of the proton on a superfine lattice
Zhouyou Fan, Xiang Gao, Ruizi Li, Huey-Wen Lin, Nikhil Karthik, Swagato Mukherjee, Peter Petreczky, Sergey Syritsyn, Yi-Bo Yang, Rui Zhang
MMSUHEP-20-007
Isovector Parton Distribution Functions of Proton on a Superfine Lattice
Zhouyou Fan, Xiang Gao,
2, 3, ∗ Ruizi Li, Huey-Wen Lin,
1, 4
Nikhil Karthik, SwagatoMukherjee, Peter Petreczky, Sergey Syritsyn,
5, 6
Yi-Bo Yang,
7, 1 and Rui Zhang
1, 4 Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA Physics Department, Tsinghua University, Beijing 100084, China Department of Computational Mathematics, Michigan State University, East Lansing, MI 48824, USA RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY, 11973, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,Chinese Academy of Sciences, Beijing 100190, China (Dated: May 26, 2020)We study isovector unpolarized and helicity parton distribution functions (PDF) of proton withinthe framework of Large Momentum Effective Theory. We use a gauge ensemble, generated by theMILC collaboration, with a superfine lattice spacing of 0 .
042 fm and a pion mass of 310 MeV,enabling us to simultaneously reach sub-fermi spatial separations and larger nucleon momenta. Wecompare the spatial dependence of quasi-PDF matrix elements in different renormalization schemeswith the corresponding results of the global fits, obtained using 1-loop perturbative matching. Wepresent determinations of first four moments of the unpolarized and helicity PDFs of proton fromthe Ioffe-time dependence of the isovector matrix elements, obtained employing a ratio-based renor-malization scheme.
I. INTRODUCTION
Decades of deep inelastic scattering (DIS) and semi-inclusive DIS (SIDIS) data over wide kinematic ranges haveprovided us insight into the structure of nucleon. Significant progress also has been made in recent years. For example,the determination of the polarized gluon distribution at small- x [1] based on the inclusive jet and pion production datafrom polarized p - p collisions at the Relativistic Heavy-Ion Collider (RHIC) [2–4] and double spin asymmetries fromopen-charm muon production at COMPASS [5], and the constraints on the polarization of sea quarks and antiquarkswith longitudinal single-spin asymmetries in W ± -boson production [6, 7]. In the future, the kinematic coverage ofnucleon PDFs will be be greatly extended by the data from from the Jefferson Lab 12-GeV program [8] and theElectron-Ion Collider (EIC) [9]. On the energy frontier, nucleon PDF not only was a critical input for the discovery ofthe Higgs boson at the Large Hadron Collider (LHC) [10, 11], but also is expected to play critical roles in determiningthe Standard-Model backgrounds during LHC’s search for physics beyond the Standard Model in future Runs 3–5.Despite great progress in the experimental and phenomenological sides, non-perturbative determinations of thePDFs starting from the microscopic theory of quantum chromodynamics (QCD) remains a challenge. To obtainthe quark PDF one has to calculate the matrix element with the quark fields are separated along the lightconebetween the hadronic states. Due to the lightcone separation, straightforward calculation of PDF is not possibleusing lattice QCD, a technique based on Euclidean time formulation. One can bypass this obstacle by calculating asimilar matrix element with spatially separated quark fields at equal time within highly boosted hadron states, whichdefines the so-called quasi-PDF (qPDF) [12, 13]. For large hadron momenta this matrix element can be related toPDF [12, 13]. The Large Momentum Effective Theory (LaMET) provides a systematic way to relate the qPDF atlarge, but finite, hadron momentum to the PDF order by order in perturbation theory [13]. Related approaches toconnect PDF to matrix elements of boosted hadrons calculable in the Euclidean time lattice computations, such as“the good lattice cross-section” [14, 15] and the pseudo-PDF [16, 17], have also been proposed. Renormalization of theunderlying boosted hadron matrix elements, usually referred as the Ioffe-time distributions (ITD), involves Wilson-line. The multiplicative renormalizability of the ITD to all orders of perturbation theory has been proven [18, 19].A practical ways to implement renormalization on the lattice, such as the use of RI-MOM scheme [20–24] andreduced Ioffe-time distributions [17], have been established. Relation between different theoretical approaches also isnow understood [25]. Based on these theoretical developments, unpolarized and polarized nucleon PDFs have beencalculated on the lattice [24, 26–33]. Furthermore, lattice calculations of the valence pion PDF also appeared [34–38]. ∗ [email protected] a r X i v : . [ h e p - l a t ] M a y The status of this field is well summarized in recent review papers [39–42]. All these calculations for the nucleon, sofar, have been carried out with lattice spacing a > .
08 fm.Having small lattice spacing plays a crucial role in calculation of PDF within the LaMET framework. To suppress thetarget mass and higher twist corrections the hadron momentum, P z , should be large. But to avoid large discretizationeffects one must ensure aP z (cid:28)
1. Furthermore, to obtain lightcone-PDF from qPDF one needs perturbative matching,which, presently, is known only up to 1-loop order. Applicability of 1-loop perturbative matching can be guaranteedonly for spatial separations z Λ QCD (cid:28) a = 0 .
042 fm. The lattice spacing used in this study is at least twice smaller thanthat used in any previous lattice calculations of the nucleon PDF. The unpolarized and helicity PDFs of the nucleonare well constrained through global fits to experimental results. Thus, we study the systematic of our calculationsby comparing P z - and z -dependence of renormalized qPDF matrix elements with the same reconstructed from thewell-known phenomenological PDFs using the LaMET framework.The rest of the paper is organized as follows. In section II we discuss the general features of LaMET and our latticesetup. In section III we discuss the nucleon 2-point functions for large values of P z and the determination of the energylevels of the fast moving nucleon. Section IV is dedicated to the analysis of the nucleon 3-point functions and thecalculations of bare qPDF. Section V describes the non-perturbative RI-MOM renormalization. The comparison ofthe lattice results on qPDF with the results of global analysis of unpolarized and helicity PDF is discussed in sectionsVI and VII, respectively. Different from RI-MOM renormalization, we discuss the analysis of ratios of nucleon matrixelements in Section VIII. Finally, section IX contains our conclusions. II. LATTICE SETUP AND LAMET
In this paper, we report the results of a lattice QCD calculation using clover valence fermions on an ensemble of N f = 2 + 1 + 1 gauge configurations with lattice spacing a = 0 .
042 fm, with space-time dimensions of 64 × M π ≈
310 MeV in the continuum limit. The gauge configurations have been generated using HighlyImproved Staggered Quarks (HISQ) [43] by the MILC Collaboration [44]. The gauge links entering the clover Wilson-Dirac operator have been smeared using hypercubic (HYP) smearing [45]. We used tree-level tadpole improved resultfor the coefficient of the clover term and the bare quark mass has been tuned to recover the lowest pion mass of thestaggered quarks in the sea [46–49]. We use only one step of HYP smearing to improve the discretization effects, sinceit is possible that multiple applications of smearing could alter the ultraviolet results for the PDF. We use multigridalgorithm [50, 51] in Chroma software package [52] to perform the inversion of the clover fermion matrix allowing uscollect relatively high statistics sample. We collected a total of 3258 measurements using 6 sources per configurationand 543 gauge configurations. In the following, we elaborate on the steps of our computation.
A. Nucleon two-point correlators
The two crucial components of the lattice computation are the two-point function and the three-point function in-volving boosted nucleon and the qPDF operator. The two point function for the nucleon boosted to spatial momentum P is the standard operator C ( t s ) = (cid:68) ˆ N s (cid:48) ( P , t s ) ˆ N † s ( P , (cid:69) , ˆ N s ( P , t ) = (cid:88) x (cid:15) abc u ( s ) a (˜ x ) (cid:16) u ( s ) b (˜ x ) T Cγ d ( s ) c (˜ x ) (cid:17) e − i P · x , (1)where ˜ x = ( x , t ) and t s is the source-sink separation along the Euclidean time direction. The index ‘ s ’ refers to thekind of quark smearing that is applied to improve the signal-to-noise of the boosted nucleon states. We either usedpoint quark operators ψ ( x ) or we used the Gaussian momentum smeared [53] for the quark fields, ψ ( s ) ( x ) that entersˆ N s , ψ ( s ) (˜ x ) = S mom ψ (˜ x ) = 11 + 6 α ψ (˜ x ) + α (cid:88) j U j (˜ x ) e i k · ˆ j ψ (˜ x + ˆ j ) , (2)where k is the momentum of the quark field, U j (˜ x ) are the gauge links in the ˆ j direction, and α is a tunable parameteras in traditional Gaussian smearing. The quark momentum should be chosen such that the signal-to-noise ratio isoptimal for the given nucleon momentum. Naively one would expect that | k | should be one third of the nucleonmomenta [53]. For this particular study, we use j = z and k z = 4 π/L , and a large Gaussian-smearing parameter α = 10. Such a momentum source is designed to align the overlap with nucleons of the desired boost momentum, andwe are able to reach higher boost momentum for the nucleon states with reasonable signals. In the nucleon two-pointcorrelators, we can study multiple values of the nucleon momentum, P = { , , P z } with P z = n z πL , n z ∈ [0 , , (3)without a significant increase in computational needs. These values of n z from 1 to 6 correspond to P z =0 . , . , . , . , .
31 and 2 .
77 GeV in physical units respectively. We either used smeared fields for both thesource and sink, which we refer to as SS, or smeared fields only for the source and point fields for the sink which werefer to as SP in the rest of the paper.
B. Nucleon three-point function
The three point function we compute is of the form C ( t s , τ ) = PP (cid:68) ˆ N s ( P , t s ) O Γ ( z ; τ ) ˆ N † s ( P , (cid:69) (4)where O Γ ( z ; τ ) is the u − d isovector qPDF operator O Γ ( z ; τ ) = (cid:88) x u (˜ x + z )Γ W z (˜ x + z, ˜ x ) u (˜ x ) − (cid:88) x d (˜ x + z )Γ W z (˜ x + z, ˜ x ) d (˜ x ) . (5) W z is the straight Wilson line along the spatial z -direction, connecting lattice sites ˜ x and ˜ x + z . The Dirac Γ used willdetermine the quantum numbers of the PDF — Γ = γ t for the unpolarized case and Γ = γ z γ for the longitudinallypolarized case. The projector operator, PP, is given by PP = γ t for the unpolarized case and PP = iγ z γ γ t for thelongitudinally polarized case, respectively. We only use smeared quark sources for the computation of C . In orderto reduce the computational cost, we only computed the C for two large values P z = 1 .
84 and 2 .
31 GeV, and forsource-sink separations t s = 16 a, a, a . C. Extraction of nucleon matrix element and perturbative matching to PDF
Using the three-point and two-point functions whose calculations are described above, we can extract the barematrix element h ( z, P z , Γ) = (cid:104) P z | O Γ ( z ) | P z (cid:105) , (6)formally in the infinite source-sink separation t s limit of their ratio R ( z, P z , Γ; τ, t s ) = C ( τ, t s ) C ( t s ) . (7)To obtain the matrix element h ( z, P z ) from the above ratio, we calculate the nucleon three-point function with insertionof O Γ ( z ) operator at three nucleon three-point source-sink separations, approximately t s = 0 . , . , .
84 fm, anddescribe its t s - and τ -dependence through 2- and 3-state ansatz. In Sec. IV, we describe our extraction of bare matrixelement from various extrapolations in detail.The next step of the computation is the renormalization of the bare matrix element h . One possible choice for O Γ is O γ z . However, for this case of Γ = γ z there is a mixing with the quark bilinear operator containing the unit matrix,Γ = 1 if Wilson fermions are used [20, 21, 54]. This mixing is absent if we use Γ = γ t , and we will use this choice forthe unpolarized PDF in this study. One way to perform the renormalization procedure on the lattice to use RI-MOMscheme [20, 22], where in the renormalized matrix element is defined as h R ( z, P z , µ R , p Rz ) = Z ( z, µ R , p Rz ) h ( z, P z , Γ) , (8)where Z is the non-perturbatively determined RI-MOM renormalization factor at a renormalization scale whosenorm is µ R and the z -component is p Rz . The dependence on p Rz arises because the z -component of the momentumnow plays a special role. We will discuss the details of the RI-MOM renormalization in section V. We will alsoconsider an alternate ratio scheme that has a well defined continuum limit in Sec. VIII. Here, the multiplicativerenormalization factor Z ratio ( z ) can be taken as the hadron matrix element at a different fixed momentum P (cid:48) z i.e., Z ratio ( z ) = ( h ( z, P (cid:48) z , Γ)) − .After the RI-MOM renormalization one obtains the renormalized matrix element h R ( z, P z , µ R , p Rz ), from which wecan define the qPDF as a function of Bjorken- x ˜ q ( x, P z , µ R , p Rz ) ≡ (cid:90) ∞−∞ dz π e ixP z z h R ( z, P z , µ R , p Rz ) . (9)From this formula it is clear that h R ( z, P z , µ R , p Rz ) can be considered as the coordinate space qPDF. For finitemomentum P z , ˜ q ( x, P z , µ R , p Rz ) has support in −∞ < x < ∞ . Unlike the physical PDF, which is frame independent,the qPDF has a nontrivial dependence on the nucleon momentum P z . When the nucleon momentum P z (cid:29) { M, Λ QCD } with M being the nucleon mass, the qPDF in RI-MOM scheme can be matched to the PDF defined in MS-scheme, q ( x, µ ) through the factorization theorem [12, 13, 25],˜ q ( x, P z , p Rz , µ R ) = (cid:90) − dy | y | C (cid:18) xy , r, yP z µ , yP z p Rz (cid:19) q ( y, µ ) + O (cid:18) M P z (cid:19) + O (cid:32) Λ P z (cid:33) , (10)where r = ( µ R /p Rz ) and C is the perturbative matching coefficient, O ( M /P z ) is the target-mass correction due tothe non-zero nucleon mass, and O (Λ /P z ) stands for higher-twist contributions. The flavor indices of q , ˜ q , and C are implied. In what follows we will discuss the non-singlet case, and therefore, mixing with gluon and sea-quarkPDFs is absent in the above formula. We use 1-loop expression of the kernel C . (The 1-loop matching including forthe singlet case also has been worked out in Ref. [55, 56].)The matching kernel C ( x, r, P z /µ, P z /p Rz ) for Γ = γ t was derived in Ref. [24] and depends on details of the RI-MOMscheme. It can be written in the following form C (cid:18) x, r, P z µ , P z p Rz (cid:19) = δ (1 − x ) + (cid:20) f , Γ (cid:18) x, P z µ (cid:19) − (cid:12)(cid:12)(cid:12)(cid:12) P z p Rz (cid:12)(cid:12)(cid:12)(cid:12) f , Γ , P (cid:18) P z p Rz ( x − , r (cid:19)(cid:21) + . (11)The subscript ‘+’ stands for the plus-prescription. Both the functions, f , Γ and f , Γ , P , depend on the choice of the Γin the operator insertion [24]. On the other hand, f , Γ is independent of the projection operator ( P ) used in definingthe RI-MOM renormalization condition, but f , Γ , P is different for different choices of the RI-MOM renormalizationcondition [24]. We also note that it is also possible to convert h R ( z, P z , µ R , p Rz ) to MS-scheme and define thecorresponding qPDF ˜ q ( x, P z , µ ) that then can be directly matched to MS PDF [22].To study the longitudinally polarized quark PDF one can use Γ = γ z γ or Γ = γ t γ . In the case Γ = γ z γ there isno mixing with quark bilinear operators with Γ = 1 [22]. Therefore, we will use this choice to study the longitudinallypolarized quark PDF and qPDF. The bare matrix element of O γ z γ can be renormalized using RI-MOM scheme andthen match to PDF in the same manner as this was done for unpolarized. The RI-MOM renormalization for thelongitudinally polarized case will be discussed in section V, while details of the matching procedure, including theformulas for f and f functions will be give in section VII. III. ANALYSIS OF THE NUCLEON TWO-POINT FUNCTION
For the extraction of the qPDF matrix element of the nucleon at large momenta it is important to understandthe contribution of different energy states to the nucleon two-point correlation function. We calculated nucleon two-point function using smeared source and smeared sink (SS correlator), as well as smeared source and point sink (SPcorrelator), for seven values of the momenta aP z = 2 π/L · n z , n z = 0 , , , , , C i ( t s , P z ), i = SS or SP, we define the effective mass aE eff ( t s , P z ) = ln (cid:32) C i ( t s /a, P z ) C i ( t s /a + 1 , P z ) (cid:33) . (12)Our results for the effective masses are shown in Fig. 1 for the SP and SS correlators.The effective mass should approach a constant corresponding to the ground state energy E ( P z ) at sufficiently large t s . The momentum dependence of the ground state energy is expected to be described by the dispersion relation E ( P z ) = (cid:112) P z + M , with M being the nucleon mass. Therefore, in Fig. 1, we show the expected ground stateenergy at different P z obtained from the dispersion relation as horizontal lines at the right for comparison. Along t s /a E e ff ( P z ) a N state = 2, pseudo N state = 2, true N state = 3 t s /a E e ff ( P z ) a N state = 2, pseudo N state = 2, true N state = 3 FIG. 1. The effective masses obtained from SP (left) and SS (right) correlators for different momenta. The bands come fromthe results of two-state ( N state = 2) and three-state ( N state = 3) fits. For N state = 2, ‘pseudo’ indicates that the effectivepseudo-plateau in range 5 a < t min < a for the first excited state E have been used, and ‘true’ indicates that the true plateauvalue of E in range t min > a have been used (see text for details). with the expected asymptotic values at large t s , we also show the t s -dependence of the the effective mass based onan effective two-state fit to the two-point function, as we will explain shortly. Indeed, we see that the effective massesapproach the corresponding values. The effective masses corresponding to the SP correlator reach a plateau at aslightly larger t s than the SS correlators. On the other hand, at small t s , the effective masses for the SP correlatorsare smaller than those for SS correlators. This implies that the contribution of the excited states is smaller for the SPcorrelator, for which a plausible reason could be that the different excited states contribute with different signs to thecorrelator. Thus, even though the ground and the excited state energies are the same in the SP and SS correlators,the two are affected differently by the higher excited states, which we can take advantage of to obtain the excitedstate spectrum reliably.In order to determine the energy levels, we fit the spectral decomposition of C ( t s ), C ( t s ) = N state − (cid:88) n =0 A n e − E n t s , (13)truncated at N state to the two-point function data over a range of values of t s between [ t min , a ]. Since the latticeextent in the time direction is 192, we did not find any effect of lattice periodicity in this range of t s to be important.We performed this fitting with one-state ( N state = 1), two-state ( N state = 2), and three-state ( N state = 3) Ans¨atze.The ground state energies, E from the fits of SS correlators for n z = 3 and 4 are shown in left panels of Fig. 2 asfunction of t min , where t min indicates that only C ( t s > t min ) have been fitted. Similar results were obtained at theother values of the momenta. The horizontal lines in the figures correspond to the results from the dispersion relationfor E . The single exponential fits give a good description of the SS correlator for t min > a , while two exponentialfits give stable results for the ground state energy already for t min > a .We found the determination of the excited state energies from the SS correlators to be more problematic thanfrom SP correlators. The excited state energy for SS is not well-constrained by simple two exponential fits, andit is also not very stable with respect to the variation of t min . Since the SP and SS correlators receive differentcontributions from excited states, we performed a combined analysis of them to obtain more reliable results for theexcited state energies. Since we were able to obtain the ground state energy E reliably from one or two exponentialfits to both the SS and SP correlators and they agree with the expectation from the dispersion relation well, weused E as a prior to performed more stable two-exponential fits. The results from the two-state exponential fits,with E as prior, for n z = 3 and n z = 4 are shown in middle and right panels Fig. 2 for the SP and SS correlators,respectively. For the SP correlators, the excited state energy E seems to approach a plateau smoothly for t min > a .It is interesting to note that, empirically, we observe the values of the plateaus agree with the dispersion relation E ( P z ) = (cid:112) P z + E ( P z = 0) , which are shown as the horizontal lines. While being an interesting observation,such a stringent identification of this state is not important to our analysis and requires further studies to rigorouslyestablish this. For the SS correlator, E develops a pseudo-plateau for 5 a < t min < a and it relaxes to the trueplateau (i.e., as identified from the SP case) for t min > a . For n z = 4, it is actually difficult to identify the trueplateau. To model the the excited state contributions to R ( τ, t s ) in the range 0 < τ < t s /
2, with t s = 16 a, a, a , t min /a E ( P z ) ( G e V ) SS, P z = 1.38 GeV E , N state = 1 E , N state = 2 t min /a E ( P z ) ( G e V ) SP, P z = 1.38 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z ) ( G e V ) SS, P z = 1.38 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z ) ( G e V ) SS, P z = 1.84 GeV E , N state = 1 E , N state = 2 t min /a E ( P z ) ( G e V ) SP, P z = 1.84 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z ) ( G e V ) SS, P z = 1.84 GeV E , N state = 2 E , N state = 3 E , N state = 3 FIG. 2. Fit results for n z = 3 (up) and n z = 4 (down) the nucleon two-point function. Left panels are for the ground state ( E )from one-state and two-state fits. Middle and right panels are for the first ( E ) and second ( E ) excited states, determined bytwo-state and three-state prior-based fits (see text for details). The horizontal lines are the values calculated from the dispersionrelation. one might consider using the well-determined values of E and E from the SP correlator at large t s . However, aswe will demonstrate now, such choices provide a less accurate description of the SS two-point function in the range5 a < t s < a . A better description of the excited-state contributions to the C (5 a < t s < a ) can be obtained byusing the effective pseudo-plateau value of E in the range of 5 a < t min < a .Since, we observe E to be well-described by a particle-like dispersion relation for sufficiently large t s , we performthree-state fits for both SP and SS correlators by imposing a prior on E as well, using its best estimate from the SPcorrelator. The results are shown in middle and right panels in Fig. 2. We see that with the prior-based three-stateexponential fits, we can obtain stable results for the first excited state energy, E ( P z ), already for relatively small t min which agrees with the dispersion relation value that we input via the prior. The value of the second excited stateis also shown in Fig. 2 and it roughly agrees with the values of E from the two-exponential fit (with prior only on E ) at smaller t min . Since the value of E is quite large, the third exponential probably corresponds to a combinationof several excited states. In Fig. 1, we show the 1- σ bands for the effective mass corresponding to: (1) Two-state fitthat uses values of E and the true value of E ; (2) two-state fit obtained by setting E to be the effective value inthe range from 5 a < t min < a ; (3) three-state fit that we described above. We find that the curves (2) and (3)agree quite well with each other in the range of 5 a < t s < a and they extrapolate in the similar fashion to theasymptotic value E . However, the curve (1) fails in capturing the data in the range 5 a < t s < a . Since for ourthree-point calculations the source-sink separations were chosen to be t s = 16 a, a, a , we must model the effectiveexcited state contributions to the three-point functions in the range 0 < τ < t s /
2. Thus, through this analysis on SPand SS correlators, we numerically demonstrate the usage of an effective value of E in the range of 5 a < t s < a that is higher than the true value of E is justified, and is the best extrapolation one could perform for the extractionof bare matrix elements in the absence of enough data to perform a three-state fit.Let us now, summarize the analysis of the nucleon two-point function. Using boosted Gaussian sources we wereable to extract ground state energy levels up to momenta 2 . P z seem to follow the continuum dispersion relation. Using this fact, we performed prior-basedfits using the energy from the dispersion relation as a prior and extracted the excited state energies as function of P z .For SP correlator the extracted value of E agrees well with the one from the dispersion relation. We show this inthe left panel of Fig. 3. Furthermore, we were able to extract an effective third energy level. These results are shownas blue points ( E ) and orange points ( E ) in the left panel of Fig. 3. Our results for the energy levels obtained fromSS as function of P z are summarized in right panel of Fig. 3. Here, the effective values of E from the pseudo-plateauand the true values are both showed. The main point of the elaborate analysis is that even though a third excitedstate contributes in the relatively shorter range of t s we use in the paper, it possible to describe the the SS correlatorvery well by a two state form with an effective value of E , which is larger than the energy of the physical excitedstate. Further details on the analysis of the two-point functions are provided in Appendix A. P z (GeV) E ( P z )( G e V ) E , SP E , SP, pr( E ) E , SP, pr( E , E ) P z (GeV) E ( P z ) ( G e V ) E , SS E , SS, pr( E , E ) E , SS, pr( E , E ) E , SS, pr( E ), pseudo FIG. 3. The energies of different states as function of P z . On the left panel, the P z dependence of E , E and E for SPcorrelators are shown. The values of E were obtained from two-state fit with prior only on E , and E from a three-statefit with priors on both E and E . On the right panel, the P z dependence of true values (blue points) and effective values of E (purple points) for SS correlators are both shown (see text for details). The lines indicate the corresponding continuumdispersion relations. IV. NUCLEON THREE-POINT CORRELATORS
In order to obtain the nucleon qPDF matrix element we consider the ratio of the 3-point function to 2-point function, R ( z, P z ; t s , τ ), at different source sink separations, t s , and operator insertion, τ . At fixed ( z, P z ) we are interested infitting the ( t s , τ )-dependence as expected from the spectral decomposition of R . If only two states contribute to thecorrelation functions, the dependence of this ratio on τ and t s is given by the following form: R fit3 ( t s , τ ) = B + e − ∆ Et s / (cid:0) B e − ∆ E ( t s / − τ ) + B e ∆ E ( t s / − τ ) (cid:1) + B e − ∆ Et s A A e − ∆ Et s . (14)Here B is the desired matrix element h , and ∆ E = E − E . Generically, B and B are independent fit parameters,except at z = 0, where B = B . If we assume that the terms proportional to A are small, the denominator can beexpanded to leading order to obtain a simpler form R fit2 ( t s , τ ) = B + e − ∆ Et s / (cid:16) B e − ∆ E ( t s / − τ ) + B e ∆ E ( t s / − τ ) (cid:17) + B e − ∆ Et s . (15)Finally, if the term proportional to B is also small compared to other terms, we get an even simpler expression thatdepends only on three parameters, B , B and B : R fit1 ( t s , τ ) = B + e − ∆ Et s / (cid:16) B e − ∆ E ( t s / − τ ) + B e ∆ E ( t s / − τ ) (cid:17) . (16)For each ( z, P z ), we fitted the ( t s , τ )-dependence of R ( z, P z ; t s , τ ) to Eqs. (14, 15, 16) and determined B in each case.In all these fits we used a fixed value of ∆ E ( P z ) = E ( P z ) − E ( P z ), with the pseudo-plateau values of E ( P z ) andthe ground state energies E ( P z ) determined from the 2-point SS correlation function, as shown in Fig. 3(right).In the following, we discuss the ratio of the 3-point function to 2-point function, R ( z, P z ; t s , τ ), and the correspondingfits for Γ = γ t and n z = 4. In Fig. 4, we show the lattice data on this ratio, together with the fit results for tworepresentative values of z , namely z = 0 and z = 8 a . The t s dependence of the lattice results is small compared to thestatistical errors. In particular, the difference between t s = 12 a and t s = 16 a data is quite small. This means that thecontribution of the excited states is not large even though the source-sink separation is below 1 fm. Given that the t s -dependence of the ratio is small, it is natural to set B = 0 since the term is suppressed by e − t s ∆ E , and perform fitsusing R fit3 ( t s , τ ). We performed fits of the lattice results using the three different fit forms above and the value of ∆ E ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 0 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 8 t s =16a t s =18a t s =20a FIG. 4. The real part ratio of the three-point function to two-point function at z = 0 (left) and z = 8 a (right) as a function of τ and for different t s . The red, orange and gray bands correspond to the bare matrix elements B extracted from fits to R fit1 , R fit2 and R fit3 , respectively. z/a R e ( h ( z , P z , t )) P z = 1.84 GeV R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a R fit2 , min =2a, t min =6a R fit3 , min =2a, t min =6a 0 5 10 15 20 25 z/a I m ( h ( z , P z , t )) P z = 1.84 GeV R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a R fit2 , min =2a, t min =6a R fit3 , min =2a, t min =6a FIG. 5. The z dependence of the bare matrix element for n z = 4. The left panel corresponds to the real part and the right panelcorresponds to the imaginary part. The different colors are the matrix elements obtained by various extrapolation methods(denoted by R fit1 , R fit2 and R fit3 ), the number of operator insertion points skipped near source and sink (denoted by τ min ), andthe t min value in two-point function fit from which the excited state E was obtained from. obtained from two-state fits of the 2-point functions with t min = 6 a . We used operator insertion time τ (cid:62) τ min = 2 a in the fits. The matrix elements, B , were obtained using the three fit functions agree within errors. The real andimaginary part of the ratio R should be symmetric and anti-symmetric with respect to z = 0 at any fixed t s and τ .In general, our lattice data is compatible with this expectation, but in certain cases we see some tension between thevalues calculated for positive and negative z . This is probably due to limited statistics. Hence, we reduce this effectby symmetrizing and antisymmetrizing the real and imaginary part of the ratio with respect to z = 0, respectively.The z -dependence of the bare matrix elements is shown in Fig. 5 for all three types of fits. We see again that all threefits give the consistent results within errors. Since B obtained from R fit3 and R fit1 are consistent with that obtainedfrom R fit1 , but with larger errors, in the following we will focus on the results obtained from R fit1 . We also carried outadditional checks for any systematic effects, as discussed below.We performed several checks to understand the systematic effects in R fit1 . First, we studied the dependence of theextracted matrix element on τ min and found no significant dependence on it. Second, we performed the fits using onlya single source-sink separation t s and compared the corresponding results from the three values of t s . Interestingly, thematrix elements calculated for source sink separation t s = 16 a, a and 20 a agree within errors, though the t s = 20 a results have very large errors. We also studied the variation of the extracting matrix element on ∆ E by using E obtained using different values of t min . We found no significant variation. Finally, we used the summation methodto obtain the matrix element. This determination has very large statistical errors but it is still compatible with allother determinations. The above checks of systematic effects are discussed further in the Appendix B. We performed p Rz (GeV) Z m p z/a = 16 Re( Z mp )-Im( Z mp ) p Rz (GeV) Z z z/a = 16 Re( Z z )-Im( Z z ) FIG. 6. The renormalization constant Z mp and Z γ z γ at µ R = 4 GeV and the Wilson link length z = 16 a ≈ .
67 fm as afunction of p Rz . a similar analysis for the three-point functions corresponding to the helicity qPDF, i.e. for Γ = γ z γ . Details of thoseanalysis also are discussed in the Appendix B. V. NONPERTURBATIVE RENORMALIZATION
We calculated the nonperturbative renormalization of the qPDF operator in the RI-MOM scheme using off-shellquark states in the Landau gauge [20, 22]. The matrix element of O Γ ( z ) in an off-shell quark state | p (cid:105) isΛ( p, z, Γ) = (cid:104) S ( p ) (cid:105) − (cid:42)(cid:88) w γ S † ( p, w + zn ) γ Γ W z ( w + zn, w ) S ( p, w ) (cid:43) (cid:104) S ( p ) (cid:105) − , (17)where n µ = (0 , , ,
1) is the unit vector along the z direction, and the summation is over all lattice sites w . Thequark propagators are defined as S ( p, x ) = (cid:88) y e ipy (cid:104) ¯ ψ ( x ) ψ ( y ) (cid:105) , S ( p ) = (cid:88) x e − ipx S ( p, x ) , (18)and γ is inserted on both sides of S † ( p, w + zn ) in Eq. 17 to get the necessary propagator (cid:80) y e − ipy (cid:104) ¯ ψ ( y ) ψ ( w + zn ) (cid:105) .For the unpolarized qPDF, we use the RI-MOM renormalization constant defined via Z mp ( z, p Rz , a − , µ R ) = Tr[ P Λ tree ( p, z, γ t )]Tr[ P Λ( p, z, γ t )] (cid:12)(cid:12)(cid:12)(cid:12) p = µ R , p z = p Rz , (19)where Λ tree ( p, z, γ t ) = γ t e − izp z is the tree level matrix element in the momentum space. Furthermore, P = γ t − ( p t /p x ) γ x is the projection operator corresponding to the so-called minimal projection, where only the term with theDirac structure proportional to γ t is kept [23, 24]. Hence, we use the subscript ‘mp’ for the renormalization constant.The renormalization constant Z mp ( z, p Rz , a − , µ R ) depends on the lattice spacing a , as well as the other two scales p Rz and µ R .We followed a similar procedure for the longitudinally polarized case, where the RI-MOM renormalization constantis defined as Z γ z γ ( z, p Rz , a − , µ R ) = Tr[ P Λ tree ( p, z, γ z γ )]Tr[ P Λ( p, z, γ z γ )] (cid:12)(cid:12)(cid:12)(cid:12) p = µ R , p z = p Rz , (20)with Λ tree ( z, p z , γ z γ ) = γ z γ e − ip z z . The projection operator P in this case was chosen to be P = γ γ z / z = 0 , , z , so it is enough at the present stage. We used the following values of the momenta for the off-shell quarkstate: ap = πL (5 , , , πL (6 , , , /
3) and πL (7 , , , / L = 64 being the spatial size the of the lattice. Thesemomenta correspond to µ R = | p | = 3 .
99 GeV, 3 .
94 GeV and 3 .
97 GeV, i.e. to µ R ∼ . z -direction and, therefore, with the abovechoice of the three momenta we have p Rz = 0 . × { , , , . . . , } GeV.0 z (fm) Z m p h ( z , P z , t ) P z = 1.84 GeV R = 4 GeV, p Rz = 0 ReIm z (fm) Z m p h ( z , P z , t ) P z = 2.31 GeV R = 4 GeV, p Rz = 0 ReIm z (fm) Z z h ( z , P z , z ) P z = 1.84 GeV R = 4 GeV, p Rz = 0 ReIm z (fm) Z z h ( z , P z , z ) P z = 2.31 GeV R = 4 GeV, p Rz = 0 ReIm
FIG. 7. Top panels: The z -dependence of the real (left) and imaginary (right) parts of the RI-MOM renormalized (modulo thewavefunction renormalization, Z q ) unpolarized qPDF matrix element. Bottom panels: Similar results for the real (left) andimaginary (right) parts of the helicity matrix element. The renormalization constant at z = 16 a ≈ .
67 fm is plotted in Fig. 6 as a function of p Rz . Due to the lineardivergence, the renormalization constant can be far from 1 at a large z ≈ .
67 fm, making the nonperturbativerenormalization unavoidable. Fig. 6 also shows that the renormalization constant will be sensitive to the value of p Rz ,while such a dependence should be canceled by the matching in the continuum. We will consider the residual p Rz dependence in the final PDF prediction as a systematic uncertainty.Having determined the renormalization constants Z mp and Z γ z γ we obtained the renormalized matrix elements,i.e., coordinate space qPDF. For the unpolarized case, h R ( z, P z , µ R , p Rz ) = Z q Z mp ( z, p Rz , a − , µ R ) h ( z, P z , γ t ) , (21)and for longitudinally polarized case,∆ h R ( z, P z , µ R , p Rz ) = Z q Z γ z γ ( z, p Rz , a − , µ R ) h ( z, P z , γ z γ ) . (22)In the above equations, Z q is the quark wavefunction renormalization factor.In Fig. 7 we show the renormalized matrix elements, modulo the factor Z q , in the RI-MOM scheme at p Rz = 0, µ R = 4 GeV. We find that the errors are large. We can achieve substantial error reductions at z (cid:54) = 0, by redefiningthe renormalized matrix elements as h R ( z, P z , µ R , p Rz ) ≡ h R ( z, P z , µ R , p Rz ) h R ( z = 0 , P z , µ R , p Rz ) , and ∆ h R ( z, P z , µ R , p Rz ) ≡ ∆ h R ( z, P z , µ R , p Rz )∆ h R ( z = 0 , P z , µ R , p Rz ) . (23)The errors of the matrix elements for z (cid:54) = 0 are reduced due to the strong correlations between z (cid:54) = 0 and z = 0 matrixelements for each gauge configurations. The effectiveness of this procedure in can be seen from Figs. 12 and 15. Asone can see the error reduction due to this division is very significant. In fact, with this method, the errors are reducedenough that the z -dependence of the matrix element is well constrained also for n z = 5. Since for the extraction of the1 x u - d = 3.2 GeV NNPDF31NNPDF31, P z = 1.84 GeVNNPDF31, P z = 2.31 GeVCT18CT18, P z = 1.84 GeVCT18, P z = 2.31 GeV FIG. 8. The NNLO isovector nucleon PDFs CT18 [57] and NNPDF3.1 [58] (solid lines), and the corresponding target-masscorrected ones (dashed lines), at a scale µ = 3 . x u - d R = 4 GeV = 3.2 GeV P z = 1.84 GeV NNPDF31qPDF, p Rz = 0qPDF, p Rz = 0.47 GeVqPDF, p Rz = 0.93 GeV 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 x u - d R = 4 GeV = 3.2 GeV P z = 2.31 GeV NNPDF31qPDF, p Rz = 0qPDF, p Rz = 0.47 GeVqPDF, p Rz = 0.93 GeV FIG. 9. qPDF corresponding to NNPDF3.1 for P z =1.84 GeV (left) and P z =2.3 GeV (right) with α s = 0 .
25 and 3 differentRI-MOM renormalization condition. qPDF we are only interested in the z -dependence of the matrix element, and we know that the unpolarized isovectornucleon matrix element at z = 0 is the isospin of the nucleon, which is unity (in our convention, c.f. Eq. 5 ) afterrenormalization, we can consider the above improved ratio of renormalized matrix elements. However, the effect oftaking this ratio is not trivial in the case of the matrix element of the helicity qPDF— the value of the renormalizedmatrix element at z = 0 should be g A ≈ .
3; this procedure is equivalent to studying a helicity PDF with the firstmoment normalized to unity, i.e., in a normalization where g A = 1. VI. UNPOLARIZED PDF: PERTURBATIVE MATCHING AND COMPARISONS WITH h R ( z , P z ) In this section, we will discuss how the renormalized coordinate-space qPDF, h R ( z, P z , µ R , p Rz ), can be related andcompared with phenomenological unpolarized nucleon PDF, such as the CT18 [57] and NNPDF3.1 [58], extractedfrom the global analysis of experimental data. The unpolarized quark PDF in the valence region is well constrainedthrough global analysis. Therefore, it is natural to start from these phenomenological PDFs as a function of Bjorken- x ,use the perturbative matching to reconstruct the corresponding coordinate-space qPDF as a function of z for different P z values, and compare with our results for h R ( z, P z , µ R , p Rz ). The reason for comparing in the z -space, rather thanconstructing the x -dependent PDF from our h R ( z, P z , µ R , p Rz ) and then comparing with the phenomenological PDFs,is the following: As can seen from Figs. 12 and 15, h R ( z, P z , µ R , p Rz ) is quite noisy for z (cid:62) . x -space is difficult to perform. Similar approach also hadbeen used for pion PDF [35].Even at the leading α s order the qPDF and the PDF differ due to the trace term in the small z -expansion [12, 25].2 zP z R e ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF31, P z = 1.84 GeVNNPDF31, P z = 2.31 GeVCT18, P z = 1.84 GeVCT18, P z = 2.31 GeV P z = 1.84 GeV P z = 2.31 GeV zP z I m ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF31, P z = 1.84 GeVNNPDF31, P z = 2.31 GeVCT18, P z = 1.84 GeVCT18, P z = 2.31 GeV P z = 1.84 GeV P z = 2.31 GeV FIG. 10. Comparisons of the real (left panel) and imaginary (right panel) parts of the ITDs, in the RI-MOM renormalizationat the scales p Rz = 0 and µ R = 4 GeV, with that obtained from the CT18 and NNPDF3.1unpolarized isovector nucleon PDFs. zP z R e ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF31, P z = 1.84 GeVNNPDF31, P z = 2.31 GeVCT18, P z = 1.84 GeVCT18, P z = 2.31 GeV 0 10 20 30 40 50 zP z I m ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF31, P z = 1.84 GeVNNPDF31, P z = 2.31 GeVCT18, P z = 1.84 GeVCT18, P z = 2.31 GeV FIG. 11. Real (left) and imaginary (right) parts of the ITDs corresponding to target-mass corrected CT18 and NNPDF3.1PDFs, for p Rz = 0, in an extended range of the Ioffe-time. This difference was explicitly calculated in Ref. [59]. In the context of DIS, such corrections have been studied longago [60], and are known as target-mass corrections. Following Ref. [59], we introduce the target-mass corrected PDF q (cid:48) ( x, P z ) = 1 √ c (cid:20) f + q (cid:18) xf + (cid:19) − f − q (cid:18) − xf − (cid:19)(cid:21) , (24)where c = M / (4 P z ) , f ± = √ c ± q ( x ) is the usual PDF that corresponds to P z → ∞ . In our analysis weuse two sets of NNLO PDF for the u and d quark and anti-quark distributions, the CT18 [57], and NNPDF3.1 [58],evaluated at scale µ = 3 . µ at which the PDF was evaluated. Since the matchingonly known to 1-loop order, we chose a scale µ = 3 . u quark is calculated as q u ( x ) = u ( x ), x > q u ( x ) = − ¯ u ( − x ), x <
0. The isovector nucleon PDF, q (cid:48) u ( x ) − q (cid:48) d ( x ) isshown in Fig. 8.In Fig. 8, we also show the target-mass corrected isovector nucleon PDF for the two momenta used in our study,namely 1.84 GeV and 2.31 GeV. We see from the figure that target-mass correction is small for the values of P z used in this study. The CT18 and NNPDF3.1 results for isovector nucleon PDF agree well for positive x , while thereis a significant difference between the two for − . < x <
0. Using the target-mass corrected NNPDF3.1 isovectornucleon PDF obtained from Eq. (24) and the 1-loop matching to RI-MOM we obtained the corresponding qPDF for P z = 1 .
84 GeV and P z = 2 .
31 GeV, µ R = 4 GeV, and p Rz = 0 , . , .
93 GeV. The functions f ,γ t and f ,γ t , mp inEq. (11) for the 1-loop matching to RI-MOM scheme with minimal projection were taken from Eq. (28) Eq. (31) ofRef. [24]. Fig. 9 shows comparisons of the NNPDF3.1 with the corresponding qPDFs. In these comparisons α s was3 z (fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0 P z = 1.84 GeV CT18NNPDF31ReIm 0.0 0.2 0.4 0.6 0.8 1.0 z (fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0.93 GeV P z = 1.84 GeV CT18NNPDF31ReIm0.0 0.2 0.4 0.6 0.8 1.0 z (fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0 P z = 2.31 GeV CT18NNPDF31ReIm 0.0 0.2 0.4 0.6 0.8 1.0 z (fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0.93 GeV P z = 2.31 GeV CT18NNPDF31ReIm
FIG. 12. Comparisons of the qPDF with the ones obtained from the global analysis for two values of the RI-MOM renor-malization scale, p Rz = 0 GeV (left row) and p Rz = 0 .
93 GeV (right row), and for two values of the nucleon boost momenta, P z = 1 .
84 GeV (upper column) and P z = 2 .
31 GeV (lower column). evaluated at scale µ = 3 . α s = 0 .
25. We see significant differences between the PDFand qPDF. For large positive x , the qPDF is larger than PDF, while for negative x qPDF can turn negative. TheqPDF strongly depends on the choice of the RI-MOM scales. It is possible to choose the RI-MOM scale such thatthe qPDF is negative for x < − .
2, even though the PDF is positive.By Fourier transforming the CT18 and NNPDF3.1 target-mass corrected qPDFs with respect to x we obtained thecorresponding distributions as a function of the so-called Ioffe-time, zP z , i.e. the corresponding ITDs [61]. Since thematching is only up to 1-loop order, the scale entering α s is not fixed. We considered three choices of the scale for α s , namely µ/ , µ, µ . The corresponding variations in the ITDs can be considered as estimates of the perturbativeuncertainties, and are shown as bands in Fig. 10. In the same figure, also we compare with the lattice results for theITDs in RI-MOM renormalization, at the renormalization scales of µ R = 4 GeV and p Rz = 0 GeV. Albeit large errors,for both values of P z the real parts of the ITDs compare well at least up to zP z (cid:46)
5. However, lattice results for theimaginary parts of ITDs undershoot the phenomenological ITDs even for zP z (cid:38) x region (c.f. Fig. 8), Fig. 10do not show any visible difference in their corresponding ITDs. To understand this better, Fig. 11 we explore theseITDs in an extended range of Ioffe-time. The difference between the PDFs in the negative- x region is only reflectedin <
10% difference in the imaginary part of the ITDs for zP z >
25, with essentially showing no difference in the realpart ITDs even up to zP z = 50.To explore the dependence of the lattice results on the choice of RI-MOM scale p Rz and the range of validityof the 1-loop matching, in Fig. 12 we show comparisons between the qPDFs as a function of z obtained in thelattice calculations and from the global analysis of PDF for two different choices of renormalization scale, namely p Rz = 0 , .
93 GeV. Very little dependence on the p Rz was observed. While the real part of the the qPDF obtainedfrom the global analysis agrees with the lattice results up to z ∼ z (cid:46) . P z = 2 .
31 GeV the agreements seem to extend to larger values of z , partly becauseof larger errors. However, it is encouraging that the central value seems to shift towards the global analysis results4 x u - d = 3.0 GeV NNPDF11polNNPDF11pol, P z = 1.84 GeVNNPDF11pol, P z = 2.31 GeVJAM17JAM17, P z = 1.84 GeVJAM17, P z = 2.31 GeV FIG. 13. NNPDF1.1pol and JAM17 isovector helicity PDF at a scale µ = 3 GeV. Also, shown are the corresponding target-masscorrected isovector helicity PDFs (dashed lines) for P z = 1 .
84 GeV and P z = 2 .
31 GeV. as P z is increased from 1 .
84 GeV to 2 .
31 GeV. In any case, at large z , we see clear tension between the imaginarypart of the lattice qPDF lattice and the results of global analysis. This suggests that the range of applicability of1-loop matching is perhaps limited to z (cid:46) . z . This observation has an important implication for our ability to described the x -dependence of PDF within the LaMET framework. For example, if the 1-loop perturbative matching works onlyfor z (cid:39) . x (cid:39) . P z (cid:38)
10 GeV.
VII. HELICITY PDF: PERTURBATIVE MATCHING AND COMPARISONS WITH ∆h R ( z , P z ) Our analysis of helicity qPDF closely follows the analysis performed in the unpolarized case, namely we startfrom the helicity PDF obtained in global analyses, reconstruct the corresponding target-mass corrected qPDF, andthen compare with the lattice results. The helicity PDF have been extracted from the global analysis by NNPDFcollaboration using DIS, inclusive W ± and jet production data from RHIC, as well as the open charm data fromCOMPAS resulting in NNPDFpol1.1 [62]. The JAM collaboration used the DIS and SIDIS data in their global analysis,combined with e + e − data to constrain the fragmentation functions at NLO [63]. The resulting PDF parameterizationis called JAM17. In Fig. 13, we show the isovector helicity PDF ∆ q u − ∆ q d . The positive- x region corresponds toquark contribution, while the negative- x region corresponds to anti-quark region. The target-mass corrected helicityPDF, ∆ q (cid:48) ( x, P z ), was obtained from helicity PDF, ∆ q ( x ), following Ref. [59]:∆ q (cid:48) ( x, P z ) = 11 + 4 c (cid:20) f + q (cid:18) xf + (cid:19) + f − q (cid:18) − xf − (cid:19)(cid:21) − (cid:90) x ±∞ c (1 + 4 c ) / (cid:20) ∆ q (cid:18) yf + (cid:19) + ∆ q (cid:18) − yf − (cid:19)(cid:21) , (25)where c = M / P z , f ± = √ c ±
1, and for the integration limits + ∞ (- ∞ ) correspond to x > x < γ γ z Λ( p, z, γ z γ )] = Tr[ γ z Λ( p, z, γ z )].Thus, the 1-loop matching of the helicity qPDF in the RI-MOM scheme with minimal projection is same as thatfor the unpolarized qPDF with Γ = γ z (instead of the Γ = γ t used before), and with the RI-MOM renormalizationcondition corresponding to the projection operator P = γ z (instead of the minimal projection). The 1-loop matchingfor the Γ = γ z operator is known for two different RI-MOM projections, the minimal projection and the /p projection,corresponding to P = γ z − ( p z /p x ) γ x and P = /p/ (4 p z ), respectively [24]. The function that depends on the RI-MOMprojection operator, i.e. f ,γ z ,γ z , entering the matching coefficient in Eq. 11 was simply deduced from these knownresults. The Lorentz structure of Λ( p, z, γ α ) for a general γ α , α = x, y, z, t is given byΛ (1) ( p, x, γ α ) = γ α (cid:104) ˜ f t ( x, ρ ) (cid:105) + + γ z p α p z (cid:104) ˜ f z ( x, ρ ) (cid:105) + + /pp α p (cid:104) ˜ f p ( x, ρ ) (cid:105) + (26)and f ,γ z , mp = ˜ f t + ˜ f z and f ,γ z ,/p = ˜ f t + ˜ f z + ˜ f p [24]. Here, the subscript ‘+’ refers to the standard plus-prescriptionand ρ = − p /p z . The functions ˜ f t , ˜ f z and ˜ f p have been calculated in Ref. [24], and we use the same notations here.5 z P z R e ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF11pol, P z = 1.84 GeVNNPDF11pol, P z = 2.31 GeVJAM17, P z = 1.84 GeVJAM17, P z = 2.31 GeV P z = 1.84 GeV P z = 2.31 GeV z P z I m ( h R ( z , P z , P R )) R = 4 GeV, p Rz = 0 NNPDF11pol, P z = 1.84 GeVNNPDF11pol, P z = 2.31 GeVJAM17, P z = 1.84 GeVJAM17, P z = 2.31 GeV P z = 1.84 GeV P z = 2.31 GeV FIG. 14. Comparisons of the real (left panel) and imaginary (right panel) parts of the isovector helicity ITDs, in the RI-MOMrenormalization at the scales p Rz = 0 and µ R = 4 GeV, with that obtained from the NNPDF1.1pol and JAM17. Therefore, for the case of P = γ z the RI-MOM projection-dependent function is given by f ,γ z ,γ z = ˜ f t + ˜ f z + ( p z /p ) ˜ f p = f ,γ z , mp + (cid:16) f ,γ z ,/p − f ,mp (cid:17) /r. (27)Thus, for the helicity qPDF the 1-loop matching RI-MOM function in the minimal projection scheme is the same asin Eq. 11, but with f ,γ z , mp given by Eq. 27, and f ,γ z , f ,γ z , mp and f ,γ z ,/p are given by Eqs. (A6-A8) of Ref. [24].Using the matching discussed above, we can obtain the isovector helicity qPDF from the target-mass correctedNNPDF1.1pol and JAM17. As before, the 1-loop matching we used α s evaluated at scale µ = 3 . µ/ µ to estimate the the scale uncertainty. We found noticeable difference betweenthe isovector helicity PDFs and the corresponding qPDFs: for positive x the qPDFs are larger than the PDFs, whilefor x < p Rz and µ R considered by us. By Fourier transforming the qPDFswe obtained the isovector helicity ITDs and compared it with our lattice results in Fig. 14. Since we normalized ourlattice results by the value of matrix element at z = 0, we normalized the phenomenological ITDs by dividing with g A = 1 .
25. Within the large statistical errors, we did not find a significant P z dependence of the lattice results. Whilethe real parts of the lattice results agree with that obtained from the phenomenological PDFs up to zP z (cid:46)
3, theimaginary parts do not agree quantitatively, but found to be qualitatively similar. We also explored the dependenceof our result on the choice of RI-MOM scales. In Fig. 15, we compare the qPDFs for µ R = 4 GeV, and p Rz = 0 GeVand p Rz = 0 .
93 GeV. From the figure, we see that the comparison between the results of lattice calculation, as well asthe global analyses are not sensitive to the choice of the renormalization scales. For both values of P z , the agreementbetween the lattice and the global analyses extends to values of | z | of about 0 . VIII. MOMENTS OF PDF FROM RATIO OF IOFFE-TIME DISTRIBUTIONS
In the previous sections, we analyzed the boosted nucleon matrix matrix elements renormalized in the RI-MOMscheme and matched it to the PDFs in the MS scheme. Due to the multiplicative renormalizability of h ( z, P z , γ t )and h ( z, P z , γ z γ ), we can form well-defined renormalized quantities by taking the ratios of matrix elements at twodifferent momenta P z and P (cid:48) z as M ( z, P z , P (cid:48) z , Γ) = h ( z, P z , Γ) h ( z, P (cid:48) z , Γ) h (0 , P (cid:48) z , Γ) h (0 , P z , Γ) . (28)The second factor on the right hand side of the above definition normalizes the z = 0 matrix element to unity, as wedid in the case of the RI-MOM scheme. The choice P (cid:48) z = 0 in the ratio is usually referred to as the reduced Ioffe-timedistributions [16], and one should think of P (cid:48) z (cid:54) = 0 as a generalization of this choice. Here, we take P z = 2 .
31 GeVand P (cid:48) z = 1 .
84 GeV, respectively. Since both P z , P (cid:48) z > Λ QCD and the nucleon mass, we expect this ratio to be simplydescribed by the leading twist expression [25], M ( z, P z , P (cid:48) z , Γ) = (cid:80) n =0 c n ( µz ) c ( µz ) ( − izP z ) n n ! (cid:104) x n (cid:105) P z ( µ ) (cid:80) n =0 c n ( µz ) c ( µz ) ( − izP (cid:48) z ) n n ! (cid:104) x n (cid:105) P (cid:48) z ( µ ) . (29)6 z(fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0 P z = 1.84 GeV NNPDF11polJAM17ReIm 0.0 0.2 0.4 0.6 0.8 1.0 z(fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0.93 GeV P z = 1.84 GeV NNPDF11polJAM17ReIm0.0 0.2 0.4 0.6 0.8 1.0 z(fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0 P z = 2.31 GeV NNPDF11polJAM17ReIm 0.0 0.2 0.4 0.6 0.8 1.0 z(fm) h R ( z , P z , P R ) R = 4 GeV, p Rz = 0.93 GeV P z = 2.31 GeV NNPDF11polJAM17ReIm
FIG. 15. Comparisons of the isovector helicity qPDF with the ones obtained from the global analysis for two values of theRI-MOM renormalization scale, p Rz = 0 GeV (left row) and p Rz = 0 .
93 GeV (right row), and for two values of the nucleon boostmomenta, P z = 1 .
84 GeV (upper column) and P z = 2 .
31 GeV (lower column).
Following Ref. [59], the target-mass corrected unpolarized PDF moments (cid:104) x n (cid:105) can be obtained by relation: (cid:104) x n (cid:105) P z (cid:104) x n (cid:105) = (cid:98) ( n +1) / (cid:99) (cid:88) i =0 C in − i +1 c i (30)and for the helicity case, (cid:104) x n (cid:105) P z (cid:104) x n (cid:105) = (cid:98) n/ (cid:99) (cid:88) i =0 (cid:18) n − i + 1 n + 1 (cid:19) C in − i c i (31)where C in is the binomial function, c = M / P z . In Eq. 29, c n ( µz ) is the 1-loop order Wilson coefficients in theMS scheme. The Wilson coefficients describes the z dependence of the twist-2 local operator associated with the n th moment of the PDF, (cid:104) x n (cid:105) ( µ ), in the MS scheme and at a factorization scale µ . As in our RI-MOM analysis, we willuse µ = 3 . µ = 3 GeV for the helicity case in the following analysis.Now, we can perform an independent analysis that avoids the usage of RI-MOM procedure completely and comparethe outcome to the prediction for M ( z, P z , P (cid:48) z , Γ) from the knowledge of NNPDF and CTEQ PDF moments. Weperform such a comparison in Fig. 16. For this, we used the values of (cid:104) x n (cid:105) ( µ ) up to an order n = n max for NNPDF31 inEq. 29, , and the complete result for CT18, to obtain the phenomenological expectation for the ratio M ( z, P z , P (cid:48) z , γ t ).The results obtained by using the truncation order n max = 2 , , ,
20 using the NNPDF31 values for (cid:104) x n (cid:105) are shown asdifferent colored bands in Fig. 16. It is clear that inclusion of up to n max = 20 moments is sufficient for convergence tothe correct PDF within z (cid:54) . z < . N = 4 is sufficient to describe the lattice results. This gives us an idea of which moments are being probedby our lattice data at different z . We observe some discernible differences between the phenomenological expectationsand our lattice M ( z, P z , P (cid:48) z , γ t ) for z > . z (fm) ( z , P z , P z , t ) NNPDF31, O2NNPDF31, O3NNPDF31, O4NNPDF31, O20CT18, O20Re( ( z , P z , P z , t )) 0.0 0.1 0.2 0.3 0.4 0.5 z (fm) ( z , P z , P z , t ) NNPDF31, O2NNPDF31, O3NNPDF31, O4NNPDF31, O20CT18, O20Im( ( z , P z , P z , t )) FIG. 16. The real (left) and imaginary(right) parts of M ( z, P z , P (cid:48) z , γ t ) is shown for P z = 2 .
31 GeV and P (cid:48) z = 1 .
84 GeV. Thedata points are from our lattice calculations, whereas the various colored bands are the corresponding results from the isovectorunpolarized PDFs from NNPDF3.1 and CT18. The band in these phenomenological expectations arise due to variations of α s ( µ ) within the scale scale µ/ µ . For NNPDF3.1, we also show results by truncating the expansion in Eq. 29 at variousorders, n = n max , in the PDF moments; these results are denoted by ‘O n max ’. understand this, we estimate the values of the moments (cid:104) x n (cid:105) that best describe our lattice data. To avoid overfittingthe data, we truncate the expansion in Eq. 29 at most by n = 4. In order to avoid lattice artifacts that might bepresent for z of the order of lattice spacing, we fit the data only from z = 2 a to a value z max . The variation of the bestfit values of (cid:104) x n (cid:105) with z max is a source of systematic error. In Fig. 17, we show the z max dependence of our estimatesfor (cid:10) x (cid:11) , (cid:10) x (cid:11) , (cid:10) x (cid:11) and (cid:10) x (cid:11) . From Fig. 16, we note the noisy determination of the imaginary part of M . As aconsequence, we find our estimates of (cid:10) x (cid:11) and (cid:10) x (cid:11) to be noisy as well. On the contrary, we were able to determine (cid:10) x (cid:11) and (cid:10) x (cid:11) reasonably well. In addition to z max dependence, we also studied whether our determination of themoments is affected by the order of truncation used in Eq. 29. We observe no significant variations with truncation.For comparison, the NNPDF and CT18 values of these moments are shown by the horizontal lines. Further, when wefix the values of (cid:10) x (cid:11) and (cid:10) x (cid:11) from NNPDF to reduce the number of fit parameters, we find the estimates for (cid:10) x (cid:11) tobe slightly elevated in value and in the direction away from NNPDF,CT18 value. To a small extent, this is seen in (cid:10) x (cid:11) as well. Thus, the observed difference between our lattice result and the NNPDF, CT18 results could be attributedto this tendency for our lattice values of (cid:10) x (cid:11) , (cid:10) x (cid:11) to be slightly higher than the corresponding phenomenologicalvalues.We repeated similar analysis for the helicity matrix element, Γ = γ γ z . In this case, the Wilson coefficients c n ( µz )are the same as in the case of unpolarized case with Γ = γ z . Since we are setting the value of the matrix elements at z = 0 to be 1 through the ratio, we only obtain the values of (cid:104) x n (cid:105) / (cid:10) x (cid:11) in the expansion Eq. 29, with (cid:10) x (cid:11) = g A . InFig. 18, we compare the results corresponding to the NNPDF11pol and JAM17 with the lattice result for the ratio. Asin the case of the unpolarized matrix element, we also test the dependence of this comparison on the truncation order n max . The sensitivity to higher moments is a bit more than that for the unpolarized case, and we find convergenceat only n max = 6 at z < . (cid:10) x (cid:11) /g A , (cid:10) x (cid:11) /g A , (cid:10) x (cid:11) /g A and (cid:10) x (cid:11) /g A that describes our lattice data via Eq. 29 truncated at most by 4 th order.In Fig. 19, we show the results as a function of the largest z used in the fits, z max . Like the unpolarized PDF case, (cid:10) x (cid:11) /g A is noisy, but seems agree with the global fit results. The more precisely determined value of (cid:10) x (cid:11) /g A is quiterobust to various ways of fitting the data and agrees nicely with the global fit values. IX. SUMMARY AND CONCLUSIONS
In this paper we studied isovector unpolarized and helicity PDFs of proton using the LaMET approach. The latticecalculations have been performed for an unphysically large pion mass of 310 MeV. On the other hand, our latticestudy was carried out using lattice spacing a = 0 .
042 fm, which is the smallest lattice spacing used in such studies.We argued that such small lattice spacing is essential for the validity of 1-loop perturbative matching between PDFand qPDF, which is a key ingredient of LaMET.Extracting the nucleon matrix elements for such large momenta and small lattice spacing is challenging because8 z max (fm) x O2, z min / a = 2O3, z min / a = 2O4, z min / a = 2 0.05 0.10 0.15 0.20 0.25 0.30 z max (fm) x O2, z min / a = 2O3, z min / a = 2O4, z min / a = 2 O2, z min / a = 2, Fix( x )O4, z min / a = 2, Fix( x , x )0.10 0.15 0.20 0.25 0.30 z max (fm) x O3, z min / a = 2O4, z min / a = 2 0.10 0.15 0.20 0.25 0.30 z max (fm) x O4, z min / a = 2O4, z min / a = 2, Fix( x , x ) FIG. 17. The moments of isovector unpolarized PDF, (cid:10) x (cid:11) (top-left), (cid:10) x (cid:11) (top-right), (cid:10) x (cid:11) (bottom-left), and (cid:10) x (cid:11) (bottom-right), that best describes the ratio M ( z, P z , P (cid:48) z , γ t ) with P z = 2 .
31 GeV and P (cid:48) z = 1 .
84 GeV. In each of the panels, the moment (cid:104) x n (cid:105) is shown as a function of z max of the fit using the functional form in Eq. 29 over a range [2 a, z max ] of the data. The resultsfrom fits using only moments up to n = 2 as free parameters in Eq. 29 are labeled ‘O2’, and those up to n = 4 are labeled ‘O4’.The results from fits that fix the moment (cid:10) x (cid:11) , or (cid:10) x (cid:11) and (cid:10) x (cid:11) , to their global fit values are also shown. For comparisons,results from CT18 and NNPDF3.1 are shown as the horizontal lines. of poor signal to noise ratio. To deal with this problem we performed a detailed study of the nucleon two-pointfunction with momentum smeared source and sink, as well as with momentum smeared source and point sink tobetter control the excited state contributions. We demonstrated that the ground state can be reliably isolated up tothe highest momenta used in this study. Furthermore, for the Euclidean time separations used that are relevant forour lattice analysis the two-point function is very well described by the ground state and and an ’effective’ excitedstate contribution, with the energy that is larger than the true excited state energy. Therefore, we argued that thetwo-state Ans¨atze is sufficient to describe the dependence of the 3-point function on the source-sink separation andon the operator insertion time. We showed that the qPDF matrix elements can be extracted in this way, and theresults do no depend on the choices of the fit interval used in our study, demonstrating the robustness of our analysisprocedure.After non-perturbative RI-MOM renormalizations we compared the lattice calculations of the spatial, z , dependenceof qPDFs with that from the phenomenological PDFs, obtained from the global pQCD-based analyses of pertinentexperimental data performed by different collaborations. Working in z -space allowed us to test the LaMET approach.The comparisons showed that there is a rough agreement between the lattice results and the results of global analysis,but only at quite small distances. Even for the very small lattice spacing used in this study, there was not enough datapoints to constrain the x -dependence of the PDFs. Instead, to translate our z -space comparisons to x -dependence, weintroduced a new ratio-based renormalization scheme for the Ioffe-time distributions. Using our lattice calculationsfor Ioffe-time distributions, renormalized via this new ratio-based scheme, we determined the first moments of theisovector unpolarized and helicity PDFs of proton, and compared these moments with that from the correspondingphenomenological PDFs.9 z (fm) ( z , P z , P z , z ) NNPDF11pol, O2NNPDF11pol, O3NNPDF11pol, O4NNPDF11pol, O20JAM17, O20Re( ( z , P z , P z , z )) 0.0 0.1 0.2 0.3 0.4 0.5 z (fm) ( z , P z , P z , z ) NNPDF11pol, O2NNPDF11pol, O3NNPDF11pol, O4NNPDF11pol, O20JAM17, O20Im( ( z , P z , P z , z )) FIG. 18. The real (left) and imaginary(right) parts of M ( z, P z , P (cid:48) z , γ γ z ) is shown for P z = 2 .
31 GeV and P (cid:48) z = 1 .
84 GeV.The data points are from our lattice calculations, whereas the various colored bands are the corresponding results from theisovector helicity PDFs from NNPDF1.1pol and JAM17. For NNPDF1.1pol, we also show results by truncating the expansionin Eq. 29 at various orders, n = n max , in the PDF moments; these results are denoted by ‘O n max ’. z max (fm) x / g A O2, z min / a = 2O3, z min / a = 2O4, z min / a = 2 0.05 0.10 0.15 0.20 0.25 0.30 z max (fm) x / g A O2, z min / a = 2O3, z min / a = 2O4, z min / a = 2 O2, z min / a = 2, Fix( x )O4, z min / a = 2, Fix( x , x )0.10 0.15 0.20 0.25 0.30 z max (fm) x / g A O3, z min / a = 2O4, z min / a = 2 0.10 0.15 0.20 0.25 0.30 z max (fm) x / g A O4, z min / a = 2O4, z min / a = 2, Fix( x , x ) FIG. 19. The moments of helicity PDF, (cid:10) x (cid:11) (top-left), (cid:10) x (cid:11) (top-right), (cid:10) x (cid:11) (bottom-left), and (cid:10) x (cid:11) (bottom-right), thatbest describes the ratio M ( z, P z , P (cid:48) z , γ γ z ) with P z = 2 .
31 GeV and P (cid:48) z = 1 .
84 GeV. In each of the panels, the moment (cid:104) x n (cid:105) isshown as a function of z max of the fit using the functional form in Eq. 29 over a range [2 a, z max ] of the data. The results fromfits using only moments up to n = 2 as free parameters in Eq. 29 are labeled ‘O2’, and those up to n = 4 are labeled ‘O4’. Theresults from fits that fix the moment (cid:10) x (cid:11) , or (cid:10) x (cid:11) and (cid:10) x (cid:11) , to their global fit values are also shown. For comparisons, resultsfrom JAM17 and NNPDF1.1pol are shown as the horizontal lines. ACKNOWLEDGMENTS
This work was supported by: (i) The U.S. Department of Energy, Office of Science, Office of Nuclear Physicsthrough the Contract No. de-sc0012704; (ii) The U.S. Department of Energy, Office of Science, Office of NuclearPhysics and Office of Advanced Scientific Computing Research within the framework of Scientific Discovery throughAdvance Computing (SciDAC) award Computing the Properties of Matter with Leadership Computing Resources;(iii) The Brookhaven National Laboratorys Laboratory Directed Research and Development (LDRD) project No.16-37. (iv) The work of ZF, RL, HL, YY and RZ are supported by the US National Science Foundation under grantPHY 1653405 CAREER: Constraining Parton Distribution Functions for New-Physics Searches. (v) The work of XGis partially supported by the NSFC Grant Number 11890712. (vi) SS also acknowledges support by the RHIC PhysicsFellow Program of the RIKEN BNL Research Center and by the National Science Foundation under CAREER AwardPHY-1847893.This research used awards of computer time provided by the INCITE program at Oak Ridge Leadership ComputingFacility, a DOE Office of Science User Facility operated under Contract No. DE-AC05-00OR22725.We are grateful to Yong Zhao for advice on perturbative matching of the helicity qPDF. [1] D. de Florian, R. Sassot, M. Stratmann, and W. Vogelsang, Phys. Rev. Lett. , 012001 (2014), arXiv:1404.4293 [hep-ph].[2] L. Adamczyk et al. (STAR), Phys. Rev. Lett. , 092002 (2015), arXiv:1405.5134 [hep-ex].[3] A. Adare et al. (PHENIX), Phys. Rev.
D90 , 012007 (2014), arXiv:1402.6296 [hep-ex].[4] A. Adare et al. (PHENIX), Phys. Rev.
D93 , 011501 (2016), arXiv:1510.02317 [hep-ex].[5] C. Adolph et al. (COMPASS), Phys. Rev.
D87 , 052018 (2013), arXiv:1211.6849 [hep-ex].[6] L. Adamczyk et al. (STAR), Phys. Rev. Lett. , 072301 (2014), arXiv:1404.6880 [nucl-ex].[7] A. Adare et al. (PHENIX), Phys. Rev.
D93 , 051103 (2016), arXiv:1504.07451 [hep-ex].[8] J. Dudek et al. , Eur. Phys. J.
A48 , 187 (2012), arXiv:1208.1244 [hep-ex].[9] A. Accardi et al. , Eur. Phys. J.
A52 , 268 (2016), arXiv:1212.1701 [nucl-ex].[10] S. Chatrchyan et al. (CMS), Science , 1569 (2012).[11] G. Aad et al. (ATLAS), Science , 1576 (2012).[12] X. Ji, Phys. Rev. Lett. , 262002 (2013), arXiv:1305.1539 [hep-ph].[13] X. Ji, Sci. China Phys. Mech. Astron. , 1407 (2014), arXiv:1404.6680 [hep-ph].[14] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. D98 , 074021 (2018), arXiv:1404.6860 [hep-ph].[15] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. , 022003 (2018), arXiv:1709.03018 [hep-ph].[16] A. V. Radyushkin, Phys. Rev.
D96 , 034025 (2017), arXiv:1705.01488 [hep-ph].[17] K. Orginos, A. Radyushkin, J. Karpie, and S. Zafeiropoulos, Phys. Rev.
D96 , 094503 (2017), arXiv:1706.05373 [hep-ph].[18] X. Ji, J.-H. Zhang, and Y. Zhao, Phys. Rev. Lett. , 112001 (2018), arXiv:1706.08962 [hep-ph].[19] T. Ishikawa, Y.-Q. Ma, J.-W. Qiu, and S. Yoshida, Phys. Rev.
D96 , 094019 (2017), arXiv:1707.03107 [hep-ph].[20] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Phys. Rev.
D97 , 014505 (2018),arXiv:1706.01295 [hep-lat].[21] M. Constantinou and H. Panagopoulos, Phys. Rev.
D96 , 054506 (2017), arXiv:1705.11193 [hep-lat].[22] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and F. Steffens, Nucl. Phys.
B923 , 394 (2017), arXiv:1706.00265 [hep-lat].[23] I. W. Stewart and Y. Zhao, Phys. Rev.
D97 , 054512 (2018), arXiv:1709.04933 [hep-ph].[24] Y.-S. Liu et al. (Lattice Parton), Phys. Rev. D , 034020 (2020), arXiv:1807.06566 [hep-lat].[25] T. Izubuchi, X. Ji, L. Jin, I. W. Stewart, and Y. Zhao, Phys. Rev.
D98 , 056004 (2018), arXiv:1801.03917 [hep-ph].[26] Y.-S. Liu, J.-W. Chen, L. Jin, R. Li, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, (2018), arXiv:1810.05043 [hep-lat].[27] H.-W. Lin, J.-W. Chen, X. Ji, L. Jin, R. Li, Y.-S. Liu, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Phys. Rev. Lett. ,242003 (2018), arXiv:1807.07431 [hep-lat].[28] J.-W. Chen, L. Jin, H.-W. Lin, Y.-S. Liu, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, (2018), arXiv:1803.04393 [hep-lat].[29] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev.
D99 , 114504 (2019), arXiv:1902.00587 [hep-lat].[30] C. Alexandrou, K. Cichy, M. Constantinou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev.
D98 , 091503 (2018),arXiv:1807.00232 [hep-lat].[31] C. Alexandrou, K. Cichy, M. Constantinou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev. Lett. , 112001(2018), arXiv:1803.02685 [hep-lat].[32] B. Jo´o, J. Karpie, K. Orginos, A. Radyushkin, D. Richards, and S. Zafeiropoulos, JHEP , 081 (2019), arXiv:1908.09771[hep-lat].[33] B. Jo, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, and S. Zafeiropoulos, (2020), arXiv:2004.01687 [hep-lat].[34] J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Schr, and Y. Zhao, Phys. Rev. D100 , 034505 (2019), arXiv:1804.01483[hep-lat]. [35] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, and S. Syritsyn, Phys. Rev. D100 ,034516 (2019), arXiv:1905.06349 [hep-lat].[36] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu, and D. G. Richards, Phys. Rev.
D99 , 074507 (2019),arXiv:1901.03921 [hep-lat].[37] B. Jo´o, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Phys. Rev.
D100 ,114512 (2019), arXiv:1909.08517 [hep-lat].[38] R. S. Sufian, C. Egerer, J. Karpie, R. G. Edwards, Jont, Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards, (2020),arXiv:2001.04960 [hep-lat].[39] Y. Zhao, Int. J. Mod. Phys.
A33 , 1830033 (2019), arXiv:1812.07192 [hep-ph].[40] K. Cichy and M. Constantinou, Adv. High Energy Phys. , 3036904 (2019), arXiv:1811.07248 [hep-lat].[41] C. Monahan,
Proceedings, 36th International Symposium on Lattice Field Theory (Lattice 2018): East Lansing, MI, UnitedStates, July 22-28, 2018 , PoS
LATTICE2018 , 018 (2018), arXiv:1811.00678 [hep-lat].[42] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, (2020), arXiv:2004.03543 [hep-ph].[43] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD,UKQCD), Phys. Rev.
D75 , 054502 (2007), arXiv:hep-lat/0610092 [hep-lat].[44] A. Bazavov et al. (MILC), Phys. Rev.
D87 , 054505 (2013), arXiv:1212.4768 [hep-lat].[45] A. Hasenfratz and F. Knechtli, Phys. Rev.
D64 , 034504 (2001), arXiv:hep-lat/0103029 [hep-lat].[46] R. Gupta, Y.-C. Jang, H.-W. Lin, B. Yoon, and T. Bhattacharya, Phys. Rev.
D96 , 114503 (2017), arXiv:1705.06834[hep-lat].[47] T. Bhattacharya, V. Cirigliano, S. Cohen, R. Gupta, A. Joseph, H.-W. Lin, and B. Yoon (PNDME), Phys. Rev.
D92 ,094511 (2015), arXiv:1506.06411 [hep-lat].[48] T. Bhattacharya, V. Cirigliano, R. Gupta, H.-W. Lin, and B. Yoon, Phys. Rev. Lett. , 212002 (2015), arXiv:1506.04196[hep-lat].[49] T. Bhattacharya, S. D. Cohen, R. Gupta, A. Joseph, H.-W. Lin, and B. Yoon, Phys. Rev.
D89 , 094502 (2014),arXiv:1306.5435 [hep-lat].[50] R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, and C. Rebbi,Phys. Rev. Lett. , 201602 (2010), arXiv:1005.3043 [hep-lat].[51] J. C. Osborn, R. Babich, J. Brannick, R. C. Brower, M. A. Clark, S. D. Cohen, and C. Rebbi,
Proceedings, 28thInternational Symposium on Lattice field theory (Lattice 2010): Villasimius, Italy, June 14-19, 2010 , PoS
LATTICE2010 ,037 (2010), arXiv:1011.2775 [hep-lat].[52] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),
Lattice field theory. Proceedings, 22nd International Sympo-sium, Lattice 2004, Batavia, USA, June 21-26, 2004 , Nucl. Phys. Proc. Suppl. , 832 (2005), [,832(2004)], arXiv:hep-lat/0409003 [hep-lat].[53] G. S. Bali, B. Lang, B. U. Musch, and A. Schaefer, Phys. Rev.
D93 , 094515 (2016), arXiv:1602.05525 [hep-lat].[54] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, (2017), arXiv:1710.01089 [hep-lat].[55] W. Wang, J.-H. Zhang, S. Zhao, and R. Zhu, Phys. Rev.
D100 , 074509 (2019), arXiv:1904.00978 [hep-ph].[56] J.-H. Zhang, X. Ji, A. Schr, W. Wang, and S. Zhao, Phys. Rev. Lett. , 142001 (2019), arXiv:1808.10824 [hep-ph].[57] T.-J. Hou et al. , (2019), arXiv:1912.10053 [hep-ph].[58] R. D. Ball et al. (NNPDF), Eur. Phys. J.
C77 , 663 (2017), arXiv:1706.00428 [hep-ph].[59] J.-W. Chen, S. D. Cohen, X. Ji, H.-W. Lin, and J.-H. Zhang, Nucl. Phys.
B911 , 246 (2016), arXiv:1603.06664 [hep-ph].[60] O. Nachtmann, Nucl. Phys.
B63 , 237 (1973).[61] B. L. Ioffe, Phys. Lett. , 123 (1969).[62] E. R. Nocera, R. D. Ball, S. Forte, G. Ridolfi, and J. Rojo (NNPDF), Nucl. Phys.
B887 , 276 (2014), arXiv:1406.5539[hep-ph].[63] J. J. Ethier, N. Sato, and W. Melnitchouk, Phys. Rev. Lett. , 132001 (2017), arXiv:1705.05889 [hep-ph]. t min /a E ( P z )( G e V ) SP, P z = 0 E , N state = 1 E , N state = 2 t min /a E ( P z )( G e V ) SP, P z = 0.46 GeV E , N state = 1 E , N state = 2 t min /a E ( P z )( G e V ) SP, P z = 0.92 GeV E , N state = 1 E , N state = 2 t min /a E ( P z )( G e V ) SP, P z = 1.38 GeV E , N state = 1 E , N state = 2 t min /a E ( P z )( G e V ) SP, P z = 1.84 GeV E , N state = 1 E , N state = 2 t min /a E ( P z )( G e V ) SP, P z = 2.31 GeV E , N state = 1 E , N state = 2 FIG. 20. Ground state energy from unconstrained one state fit and two state fit of the SP correlators. t min /a E ( P z )( G e V ) SP, P z = 0.46 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z )( G e V ) SP, P z = 0.92 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z )( G e V ) SP, P z = 2.31 GeV E , N state = 2 E , N state = 3 E , N state = 3 FIG. 21. The energies of the first ( E ) and second ( E ) excited states from constrained two-state and three-state fits of SPcorrelator for n z = 1 (left), n z = 2 (middle) and n z = 5 (right). The horizontal line is the values calculated from the dispersionrelation. Appendix A: Analysis of the nucleon two point function
In this appendix we discuss some details of the analysis of the SP and SS two point correlators. In Fig. 20 we showthe ground state energy from one and two exponential fits of SP correlators as function of t min . Contrary to the fitsof the SS correlators stable result for the ground state energy, E is only obtained for t min ≥ p z . In Fig.21 we show the results on E ( p z ) for n z = 1 , E approaches the value expected from the dispersion relation for t min >
11 if two exponentialfit is used. For constrained three exponential fits the same value is approached for t min = 2. In Fig. 22 we show theamplitudes, A i , i = 1 , , ... , of different states normalized by the value of the two-point correlator at t = 0, which bydefinitions is equal to (cid:80) i A i . We see that A is slightly higher than A , while A is significantly larger than either A or A .Similar analysis was performed for SS correlators and the results for the excited state energies and amplitudes areshown in Fig. 23 and Fig. 24, respectively. From these figures we see that a pseudo-plateau develops for the firstexcited states for 5 < t min <
10 of 2-state fit. We see that A and A are similar in this case, and A decreases as t min increasing. Appendix B: Analysis of the three point function
In this appendix we discuss further details of the extraction of the bare matrix element of the qPDF operator. First,we show our results for the ratio of the 3-point function to two point function for different source sink separation anddifferent values of z as function of the operator insertion time τ in Fig. 25 for n z = 4. In this figure we also show the3 t min /a A n SP, P z = 0 GeV A A A t min /a A n SP, P z = 0.46 GeV A A A t min /a A n SP, P z = 0.92 GeV A A A t min /a A n SP, P z = 1.38 GeV A A A t min /a A n SP, P z = 1.84 GeV A A A t min /a A n SP, P z = 2.31 GeV A A A FIG. 22. The amplitudes of different states obtained from constrained 3-state fit of SP correlator and normalized by C SP pt ( t = 0)as function of temperature. t min /a E ( P z )( G e V ) SS, P z = 0.46 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z )( G e V ) SS, P z = 0.92 GeV E , N state = 2 E , N state = 3 E , N state = 3 t min /a E ( P z )( G e V ) SS, P z = 2.31 GeV E , N state = 2 E , N state = 3 E , N state = 3 FIG. 23. The energies of the first ( E ) and second ( E ) excited states from constrained two-state and three-state fits of SScorrelator for n z = 1 (left), n z = 2 (middle) and n z = 5 (right). The horizontal line is the values calculated from the dispersionrelation. t min /a A n SS, P z = 0 GeV A A A t min /a A n SS, P z = 0.46 GeV A A A t min /a A n SS, P z = 0.92 GeV A A A t min /a A n SS, P z = 1.38 GeV A A A t min /a A n SS, P z = 1.84 GeV A A A t min /a A n SS, P z = 2.31 GeV A A A FIG. 24. The amplitudes of different states obtained from constrained 3-state fit of SS correlator and normalized by C SS pt ( t = 0)as function of temperature. ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 12 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 1.84 GeV, z/a = 12 t s =16a t s =18a t s =20a FIG. 25. The ratio of the 3-point function to the 2-point function for z=4,8,12 and n z = 4. The upper panels show the realpart, while the imaginary part is shown in the lower panels. The results of R fit are shown as lines. ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 12 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , t ) P z = 2.31 GeV, z/a = 12 t s =16a t s =18a t s =20a FIG. 26. The ratio of the 3-point function to the 2-point function for z=4,8,12 and n z = 5. The upper panels show the realpart, while the imaginary part is shown in the lower panels. The results of R fit are shown as lines. results for R fit . As one can see from the figure R fit can describe the data well for all values of t . In Fig. 26 we showthe same analysis but for n z = 5.As discussed in the main text we performed R fit using single value of source sink separation. The results are shownin Fig. 27 for the real part of the matrix element. As one can see from the figure the results obtained from this fitfor t = 16, 18 and 20 agree within errors. We performed fits using the form f it with τ > τ min and taking the valueof E from the 2-point function fit with t > t min . The results are shown in Fig. 28. We see no significant dependenceon τ min and t min .Another way to obtain the matrix element is to use the summation method. The summation method is illustratedin Fig. 29 for n z = 4. The results obtained from the summation method agree with those from R fit but have muchlarger errors. The statistical errors of the n z = 5 data are too large to use the summation method. Furthermore, wecould also reduce the error in the summation method by dividing by the matrix element at z = 0 as can be seen inFig. 30Similar analysis of the ratio if the three point function to two point function was carries out for longitudinallypolarized qPDF operator. The results are summarized in Figs. 31, 32 33. To take the advatage of correlationbetween different z and cancel the field renormalization factor, we divided the bare matrix elements by the matrix5 z/a R e ( h ( z , P z , t )) P z = 1.84 GeV t s =16a t s =18a t s =20a z/a R e ( h ( z , P z , t )) P z = 2.31 GeV t s =16a t s =18a t s =20a FIG. 27. The z -dependence of the qPDF matrix element obtained using R fit with single value of the source sink separationfor n z = 4 , z/a R e ( h ( z , P z , t )) P z = 1.84 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a 0 5 10 15 20 25 z/a I m ( h ( z , P z , t )) P z = 1.84 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a0 5 10 15 20 25 z/a R e ( h ( z , P z , t )) P z = 2.31 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a 0 5 10 15 20 25 z/a I m ( h ( z , P z , t )) P z = 2.31 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a R fit1 , min =3a, t min =6a FIG. 28. Real (left) and imaginary (right) parts of the bare matrix as function of z . The top panel show the result for n z = 4,the bottom panel show the results for n z = 5. The results for different choices of τ min and t min in the 2-point function fits areshown. element at z = 0. The errors are much smaller after this division as discussed in the main text.6 t s /a S U M z=0, P z = 1.84 GeV SUM(1) z/a R e ( h ( z , P z , t )) P z = 1.84 GeV SUM(1)
FIG. 29. The t dependence of the sum of the ratio of the three point function to two point function (left) and the z dependenceof the matrix element extracted from the summation method (right). SUM(n) means summation fit with n skipped timeinsertion. z/a R e ( h ( z , P z , t )) P z = 1.84 GeV SUM(1)
FIG. 30. The z -dependence of the real part of the bare qPDF matrix element obtained by summation method after divisionby the matrix element for z = 0 at n z = 4. SUM(n) means summation fit with n skipped time insertion. ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 12 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 1.84 GeV, z/a = 12 t s =16a t s =18a t s =20a FIG. 31. The ratio of the 3-point function to the 2-point function corresponding to helicity qPDF for z=4,8,12 and n z = 4.The upper panels show the real part, while the imaginary part is shown in the lower panels. The results of R fit are shown aslines. ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 12 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 4 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 8 t s =16a t s =18a t s =20a ( t s /2) a R ( t , ; z , P z , z ) P z = 2.31 GeV, z/a = 12 t s =16a t s =18a t s =20a FIG. 32. The ratio of the 3-point function to the 2-point function corresponding to helicity qPDF for z = 4 , ,
12 and n z = 5.The upper panels show the real part, while the imaginary part is shown in the lower panels. The results of R fit are shown aslines. z/a R e ( h ( z , P z , z )) P z = 1.84 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a 0 5 10 15 20 25 z/a I m ( h ( z , P z , z )) P z = 1.84 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a0 5 10 15 20 25 z/a R e ( h ( z , P z , z )) P z = 2.31 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a 0 5 10 15 20 25 z/a I m ( h ( z , P z , z )) P z = 2.31 GeV R fit1 , min =2a, t min =4a R fit1 , min =2a, t min =5a R fit1 , min =2a, t min =6a FIG. 33. The real (left) and the imaginary (right) parts of the bare matrix corresponding to helicity qPDF. The upper panelscorrespond to n z = 4, while the lower panels correspond to n zz