Nucleon Momentum Fraction, Helicity and Transversity from 2+1-flavor Lattice QCD

High statistics results for the isovector momentum fraction, $\langle x \rangle_{u-d}$, helicity moment, $\langle x \rangle_{\Delta u-\Delta d}$, and the transversity moment, $\langle x\rangle_{\delta u-\delta d}$, of the nucleon are presented using seven ensembles of gauge configurations generated by the JLab/W&M/LANL/MIT collaborations using $2+1$-flavors of dynamical Wilson-clover quarks. Attention is given to understanding and controlling the contributions of excited states. The final results are obtained using a simultaneous fit in the lattice spacing $a$, pion mass $M_\pi$ and the finite volume parameter $M_\pi L$ keeping leading order corrections. The data show no significant dependence on the lattice spacing and some evidence for finite-volume corrections. The main variation is with $M_\pi$, whose magnitude depends on the mass gap of the first excited state used in the analysis. Our final results, in the $\overline{\rm MS}$ scheme at 2 GeV, are $\langle x \rangle_{u-d} = 0.160(16)(20)$, $\langle x \rangle_{\Delta u-\Delta d} = 0.192(13)(20)$ and $\langle x \rangle_{\delta u-\delta d} = 0.215(17)(20)$, where the first error is the overall analysis uncertainty assuming excited-state contributions have been removed, and the second is an additional systematic uncertainty due to possible residual excited-state contributions. These results are consistent with other recent lattice calculations and phenomenological global fit values.


I. INTRODUCTION
Steady progress in both experiment and theory is providing an increasingly detailed description of the hadron structure in terms of quarks and gluons. The distributions of quarks and gluons within nucleons are being probed in experiments at the Relativistic Heavy Ion Collider (RHIC) at BNL [1,2], Jefferson Lab [3] and the Large Hadron Collider at CERN. Experiments at the planned electron-ion collider [4] will significantly extend the range of Bjorken x and Q 2 and further improve our understanding. From these data, and using higher order calculations of electroweak and strong corrections, the phenomenological analyses of experimental data (global fits) are providing parton distribution functions (PDFs) [5,6], transverse momentum dependent PDFs (TMDs) [7], and generalized parton distributions (GPDs) [8]. These distributions are not measured directly in experiments [9,10], necessitating phenomenological analyses that have involved different theoretical inputs.
Lattice QCD calculations are beginning to provide such input, and a review of the cross-fertilization between the two efforts has been presented in Refs. [11,12]. With increasing computing power and advances in algorithms, the precision of lattice QCD calculations has increased significantly and there now exist many quantities for which there is good agreement with experimental results, and for some, the lattice results are the most precise as reviewed in the recent Flavor Averaging Group (FLAG) 2019 report [13].
In this work, we present high-statistics lattice data for the isovector momentum fraction, and the helicity and transversity moments, whose calculation is now reaching precision comparable to that for nucleon charges, which are the zeroth moments of the distributions and obtained from the matrix elements of local quark bilinear operators [13,14].
Calculations for the three first moments, x u−d , x ∆u−∆d and x δu−δd , have been done on seven ensembles generated using 2+1-flavors of Wilson-clover quarks by the JLab/W&M/LANL/MIT collaborations [15]. The data at three values of lattice spacings a, two values of the pion mass, M π ≈ 170 and ≈ 270 MeV, and on a range of large volumes, characterized by M π L, allow us to carry out a simul-taneous fit in these three variables to address the associated systematic uncertainties. In the analysis of the two-and three-point correlation functions from each ensemble, we carry out a detailed investigation of the dependence of the results on the spectra of possible, a priori unresolved, excited states included in the fits to remove excited-state contamination (ESC). Concretely, the full analysis is carried out using three strategies to estimate the mass gap of the first excited state from the two-and threepoint correlation functions, and we use the spread in the results to assign a second systematic uncertainty to account for possible remaining contributions from excited-states.
Our final results, given in Eq. (20), are x u−d = 0.160 (16) (20), x ∆u−∆d = 0.192 (13) (20) and x δu−δd = 0.215 (17) (20) in the MS scheme at 2 GeV. These estimates are in good agreement with other lattice and phenomenological global fit results as discussed in Sec. VI. The most extensive and precise results from global fits are for the unpolarized moments of the nucleons, the momentum fraction x q , while those for the helicity fraction, the polarized moment x ∆q , have a large spread and our lattice results are consistent with the smaller error global fit values at the lower end. Lattice QCD results for the transversity x δq are a prediction due to lack of sufficient experimental data [11,12].
The paper is organized as follows: In Sec. II, we briefly summarize the lattice parameters and methodology. The definitions of moments and operators investigated are given in Sec. III. The twoand three-point functions calculated, and their connection to the moments, are specified in Sec. IV. The analysis of excited state contributions and the extraction of the ground state matrix elements is presented in Sec. V. Results for the moments after the chiral-continuum-finite-volume (CCFV) extrapolation are given in Sec. VI, and compared with other lattice calculations and global fit values. We end with conclusions in Sec. VII. The data and fits used to remove excited-state contamination are shown in Appendix A and the calculation of the renormalization factors, Z V D,AD,T D , for the three operators is discussed in Appendix B.

II. LATTICE METHODOLOGY
This work follows closely the methodology described in Ref. [16], with two major differences. The first is the calculation here uses 2+1-flavors of Wilson-clover fermions in a clover-on-clover unitary formulation of lattice QCD, whereas the clover-on-HISQ formulation was used in Ref. [16]. The clover action includes one iteration of stout smearing with weight ρ = 0.125 for the staples [17]. The tadpole corrected tree-level Sheikholeslami-Wohlert coefficient c SW = 1/u 0 [18], where u 0 is the fourth root of the plaquette expectation value, is very close to the nonperturbative value determined, a posteriori, using the Schrödinger functional method [19], a consequence of the stout smearing. The update of configurations was carried out using the rational hybrid Monte Carlo (RHMC) algorithm [20] as described in Ref. [21].
The parameters of the seven clover ensembles generated by the JLab/W&M/LANL/MIT collaborations [15] and used in the analysis are summarized in Table I. The range of lattice spacings covered is 0.071 ≤ a ≤ 0.127 fm and of lattice size is 3.7 ≤ M π L ≤ 6.2. So far simulations have been carried out at two pion masses, M π ≈ 270 and 170 MeV. These seven data points allow us to perform chiralcontinuum-finite-volume fits to obtain physical results.
The second improvement is higher statistics data that lead to a more robust analysis of three strategies for evaluating ESC. Table I also gives the number of configurations, the source-sink separations τ , high precision (HP) and low precision (LP) measurements made to cost-effectively increase statistics using the bias-corrected truncated-solver method [22,23].
The parameters used to construct the Gaussian smeared sources [14,16,24,25], are given in Table II. To construct the smeared source, the gauge links were first smoothened using twenty hits of the stout algorithm with ρ = 0.08 and including only the spatial staples. The root-mean-square size of the Gaussian smearing, dr r 4 S † S/ dr r 2 S † S with S(r) the value of the smeared source at radial distance r, was adjusted to be between 0.72-0.76 fm to reduce ESC. The quark propagators from these smeared sources were generated by inverting the Dirac operator (same as what was used to generate the lattices) using the multigrid algorithm [26][27][28]. These propagators are then used to construct the two-and three-point correlation functions.

III. MOMENTS AND MATRIX ELEMENTS
The first moments of spin independent (or unpolarized), q = q ↑ + q ↓ , helicity (or polarized), ∆q = q ↑ − q ↓ , and transversity, δq = q + q ⊥ distributions, are defined as  Table I. Lattice parameters of the 2+1-flavor clover ensembles generated by the JLab/W&M/LANL/MIT collaboration and analyzed in this study. We specify the lattice spacing a, pion mass Mπ, lattice size L 3 × T , the values of source-sink separation τ simulated, the number of configurations analyzed, and the total number of high precision (HP) and low precision (LP) measurements made.  )   Table II. The parameters used in the calculation of the clover propagators. The hopping parameter for the light quarks, κ l , in the clover action is given by 2κ l = 1/(m l +4). cSW is the Sheikholeslami-Wohlert improvement coefficient in the clover action. The parameters used to construct Gaussian smeared sources [29], {σ, NKG} in the Chroma convention [30], are given in the fourth column where NKG is the number of applications of the Klein-Gordon operator and σ controls the width of the smearing. The resulting root-meansquare radius of the smearing in lattice units, defined as dr r 4 S † S/ dr r 2 S † S with S(r) the value of the smeared source at radial distance r, is given in the last column.

ID
where q ↑(↓) corresponds to quarks with helicity aligned (anti-aligned) with that of a longitudinally polarized target, and q (⊥) corresponds to quarks with spin aligned (anti-aligned) with that of a transversely polarized target.
These moments, at leading twist, are extracted from the forward matrix elements of one-derivative vector, axial-vector and tensor operators within ground state nucleons at rest. The complete set of the relevant twist two operators are where q = {u, d} is the isodoublet of light quarks and σ µν = (γ µ γ ν − γ ν γ µ )/2. The derivative consists of four terms defined in Ref. [16]. Lorentz indices within { } in Eq. (4) are symmetrized and within [ ] are antisymmetrized. It is also implicit that, where relevant, the traceless part of the above operators is taken. Their renormalization is carried out nonperturbatively in the regularization independent RI -MOM scheme as discussed in Appendix B. A more detailed discussion of these twist-2 operators and their renormalization can be found in Refs. [31] and [32].
In our setup to calculate the isovector moments, we work with τ a = τ 3 and fix the spin of the nucleon state to be in the "3" direction. With these choices, the explicit operators calculated are The forward matrix elements (M E) of these operators within the ground state of the nucleon with mass M N are related to the moments as follows: The moments are, by construction, dimensionless.

IV. CORRELATION FUNCTIONS AND MOMENTS
To construct the two-and three-point correlation functions needed to calculate the matrix elements, the interpolating operator N used to create/annihilate the nucleon state is where {a, b, c} are color indices, q 1 , q 2 ∈ {u, d} and C = γ 4 γ 2 is the charge conjugation matrix in our convention. The nonrelativistic projection (1±γ 4 )/2 is inserted to improve the signal, with the plus and minus signs applied to the forward and backward propagation in Euclidean time, respectively [31]. At zero momentum, this operator couples only to the spin-1 2 states. The zero momentum two-point and three-point nucleon correlation functions are defined as where α, β are spin indices. The source is placed at time slice 0, the sink is at τ and the one-derivative operators, defined in Sec. III, are inserted at time slice t. Data have been accumulated for the values of τ specified in Table I, and for each τ for all intermediate times 0 < t < τ . To isolate the various contributions, projected 2and 3-point functions are constructed as The projector P 2pt = 1 2 (1 + γ 4 ) in the nucleon correlator gives the positive parity contribution for the nucleon propagating in the forward direction. For the connected 3-point contributions P 3pt = 1 2 (1 + γ 4 )(1 + iγ 5 γ 3 ) is used. With these spin projections, the three moments are obtained using Eqs. 8, 9 and 10. To display the data, we construct the ratios that give the ground state matrix element in the limits t → ∞ and (τ − t) → ∞. These ratios are shown in Figs. 5-10 in Appendix A. We re-emphasize that the ground state matrix element 0|O|0 used in the analysis is obtained from fits to C 3pt O (τ ; t) with input of spectral quantities from C 2pt (τ ). These fits are carried out within a single-elimination jackknife process, which is used to get both the central values and the errors.

V. CONTROLLING EXCITED STATE CONTAMINATION
A major challenge to precision results is removing the contribution of excited states in the threepoint functions. These occur because the lattice nucleon interpolating operator N , defined in Eq. (11), couples to the nucleon, all its excitations and multiparticle states with the same quantum numbers. Previous lattice calculations have shown that these ESC can be large [16,[33][34][35]. The strategy to remove these artifacts in this work is the same as described in Ref. [16]: reduce ESC by using smeared sources in the generation of quark propagators and then fit the data at multiple source-sink separations τ using the spectral decomposition of the correlation functions (Eqs. (17) and (18)) keeping as many excited states as possible without overparameterizing the fits. In this work, we examine three strategies, {4, 3 * }, {4 N π , 3 * } and {4, 2 free }, that use different estimates of the excited state masses in the fits as described below.
The spectral decomposition of the zeromomentum two-point function, C 2pt , truncated at four states, is given by We fit the data over the largest time range, {τ min − τ max }, allowed by statistics, i.e., by the stability of the covariance matrix, to extract M i and A i , the masses and the amplitudes for the creation/annihilation of the four states by the interpolating operator N . We perform two types of 4state fits. In the fit denoted {4}, we use the empirical Bayesian technique described in the Ref. [21] to stabilize the three excited-state parameters. In the second fit, denoted {4 N π }, we use a normally distributed prior for M 1 , with value given by the lower of the non-interacting energy of N (−1)π(1) or the N (0)π(0)π(0) state 1 . The masses of these two obtained from the two 4-state fits to the two-point functions. The next six columns give the values of the mass gap, a∆M1 ≡ a(M1 − M0), of the first excited state obtained from different fits studied in this work. The notation used is {2} ({4}) is a two-state (4-state) fit to the two-point functions, {4 N π } is a 4-state fit to the two-point functions with a prior with a narrow width for a∆M1 corresponding to the non-interacting N π state. In the three {2 free } cases, the a∆M1 are determined from fits to the three-point functions used to extract the three moments as explained in the text.
states are roughly equal for the seven ensembles and lower than the M 1 obtained from the {4} fit. The lower energy N (−1)π(1) state has been shown to contribute in the axial channel [36], whereas for the vector channel the N (0)π(0)π(0) state is expected to be the relevant one. Since the two states have roughly the same mass, which is all that matters in the fits, we do not distinguish between them and use the common label {4 N π }. We also emphasize that even though we use a Bayesian procedure for stabilizing the fits, the errors are calculated using the jackknife method and are thus the usual frequentist standard errors.
In the fits to the two-point functions, the {4} and {4 N π } strategies cannot be distinguished on the basis of the χ 2 /dof. In fact, the full range of M 1 values between the two estimates, from {4} and {4 N π }, are viable on the basis of χ 2 /dof alone. The same is true of the values for M 2 , indicating a large flat region in parameter space. Because of this large region of possible values for the excited-state masses, M i , we carry out the full analysis with three strategies that use different estimates of M i and investigate the sensitivity of the results on them. The ground-state nucleon mass obtained from the two fits is denoted by the common symbol M N ≡ M 0 and the successive mass gaps by ∆M i ≡ M i − M i−1 . These are given in Table III. standard correlated χ 2 plus the square of the deviation of the parameter from the prior mean normalized by the prior width. This quantity is minimized in the fits. In the following we quote this augmented χ 2 divided by the degrees of freedom calculated without reference to the prior, and call it χ 2 /dof for brevity.
The analysis of the three-point functions, C 3pt O , with insertion of the operators with zero momentum defined in Eqs. (5), (6) and (7), is performed retaining up to three states |i in the spectral decomposition: (18) To get the forward matrix element, we also fix the momentum at the sink to zero. To remove the ESC and extract the desired ground-state matrix element, 0|O|0 , we make a simultaneous fit in t and τ . The full set of values of τ investigated are given in Table  I. In choosing the set of points, {t, τ }, to include in the final fit, we attempt to balance statistical and systematic errors. First, we neglect t skip points next to the source and sink in the fits as these have the largest ESC. Next, noting that the data at smaller τ have exponentially smaller errors but larger ESC, we pick the largest three values of τ that have statistically precise data. Since errors in the data grow with τ , we partially compensate for the larger weight given to smaller τ data by choosing t skip to be the same for all τ , i.e., by including increasingly more t points with larger τ , the weight of the larger τ data points is increased. Most of our analysis uses a 3 *fit, which is a three-state fit with the term involving 2|O|2 set to zero, as it is undetermined and its inclusion results in an overparameterization based on the Akaike information criteria [37].
The key challenge to 3-state fits using Eq. (18) is determining the relevant M i to use because fits to the two-point function show a large flat region in the space of the M i with roughly the same χ 2 /dof.  Table IV. Estimates of the three unrenormalized moments from the three fit strategies, {4 N π , 3 * }, {4, 3 * } and {4, 2 free }, used to analyze the two-and three-point functions and remove ESC. For each fit strategy we give the τ values and the number of time slices t skip omitted next to the source and the sink in the final fits to the three-point functions.
Theoretically, there are many candidate intermediate states, and their contribution to the three-point functions with the insertion of operator O is not known a priori. To investigate the sensitivity of 0|O|0 to possible values of M 1 and M 2 , we carry out the full analysis with the following three strategies: state fit to the two-point function using Eq. (17) and then a {3 * } fit is made to the three-point function using Eq. (18). Both fits are made within a single jackknife loop. This is the standard strategy, which assumes that the same set of states are dominant in the two-and three-point functions.
• {4 N π , 3 * }: The excited state spectrum is taken from a 4-state fit to the two-point function but with a narrow prior for the first excited state mass taken to be the energy of a noninteracting N (p = 1)π(p = −1) state (or N (0)π(0)π(0) that has roughly the same energy). This spectrum is then used in a {3 * } fit to the three-point function. This variant of the {4, 3 * } strategy assumes that the lowest of the theoretically allowed tower of N π (or N ππ) states contribute.
• {4, 2 free }: The only parameters taken from the {4} state fit are the ground state amplitude A 0 and mass M 0 . In the two-state fit to the three-point function, the mass of the first excited state, M 1 , is left as a free parameter, ie, the most important determinant of ESC, M 1 , is obtained from the fit to the three-point function.
The mnemonic {m, n} denotes an m-state fit to the two-point function and an n-state fit to the threepoint function.
The data for the ratios R O (τ ; t) are plotted in Figs. 5-10 in Appendix A for the three operators, the three strategies, and all seven ensembles. We note the following features: • The fractional statistical errors are less than 2% for all three operators and on all seven ensembles. The only exceptions are the τ = 16 (τ = 21) data on the a094m270 and a091m170L (a071m170) ensemble.
• The errors grow, on average, by a factor between 1.3-1.5 for every two units increase in τ /a. This is smaller than the asymptotic factor, e (M N −3Mπ/2)τ , expected for nucleon correlation functions. There is also a small increase in this factor between • The data for all three operators is symmetric about t = τ /2 as predicted by the spectral decomposition. Only on three ensembles, a094m270, a091m170L and a071m170, the symmetry about t = τ /2 is not manifest in the largest τ data. As stated above, these data have the largest errors, and the deviations are within errors.
• In all cases (operators and ensembles), the convergence of the data towards the τ → ∞ value is monotonic and from above. Thus, ESC causes all three moments to be overestimated.
With the data satisfing the expected conditions, we are able to make 3 * -state (three-state fit neglecting the term with 2|O|2 ) fits in almost all cases to data with the largest three values of τ . The one exceptions is the a071m170 ensemble where we neglect the largest τ = 21 data as the statistics are still inadequate. Including it in the fits does not change the results. We have also checked that the results from fits keeping the largest four values of τ overlap with these within 1σ. The data and the result of the fit evaluated for various values of τ are shown in each panel in Figs. 5-10 along with the τ = ∞ value as the blue band. The figure labels also give the ensemble ID, the value of the moment obtained using Eqs. (8), or (9), or (10), the χ 2 /dof of the fit, and the values of τ at which data have been collected and displayed. The three panels in each row of Figs. 5-10 have the same data but show fits with the three strategies that are being compared. The scale for the yaxis is chosen to be the same for all the plots to facilitate this comparison. The values of ∆M 1 entering/determined by the various fits are given in Table III.
As mentioned above, the key parameter needed to control ESC is the mass gap ∆M 1 of the first excited state that provides the dominant contribution. Theoretically, the lightest possible state with positive parity contributing to the forward matrix elements is either N (p = 1)π(p = −1) or N (0)π(0)π(0) depending on the value of M π and the lowest momenta, which is larger than 200 MeV on all seven ensembles. For our ensembles, the non-interacting energies for these two states are roughly equal. Since the fits do not rely on knowing the identity of the state but only on ∆M 1 , we regard these two possible states as operationally the same and label them N π. Thus in the strategy {4 N π , 3 * }, ∆M 1 is approximately the lowest possible value, and accounts for the possibility that one (or both) of these states gives the dominant ESC.
The strategy {4, 3 * } assumes that the relevant states are the same in two-and three-point functions. The strategy {4, 2 free } only takes A 0 and M 0 from the two-point fit, whose determination is robust-the variation in their values between {4, 3 * } and {4 N π , 3 * } is less than a percent as shown in Table III. The value of ∆M 1 is an output in this case. The relative limitation of the {4, 2 free } strategy is that, with the current data, we can only make two-state fits to the three-point functions, ie, include only one excited state.
The data for ∆M 1 summarized in Table III displays the following qualitative features: • The values a∆M N . This suggests that the lowest excited state in the {4} fit to the two-point function is close to the N (1440).
• a∆M is significantly smaller than a∆M {4} 1 as mentioned above.
• On five ensembles, a∆M in Table III to show how much a∆M 1 can vary between a two-and four-state fit.) To check whether this rough agreement is a possibility for the remaining two ensembles, a094m270 and a071m170, we made fits with a range of priors but did not find a flat direction with respect to a∆M from these two ensembles are unexplained, however, as noted previously, the statistical errors in these two ensembles are the largest.
• The a∆M {2 free } 1 for helicity and transversity moments are roughly the same and much larger than even a∆M The unrenormalized results for the three moments obtained using Eqs. (8), or (9), or (10) are given in Table IV along with the values of {t, τ } used. The parameters and the χ 2 /dof of the fits for the various strategies are given in Tables V, VI, and VII. In these tables, we include results with {4, 2} and {4 N π , 2} in addition to the {4, 3 * }, {4 N π , 3 * } and {4, 2 free } strategies to show that the variation on including the second excited state is small, ie, ∆M 1 is the key parameter in controlling ESC.
We draw the following conclusions from the results presented in Tables III-VII and the fits shown in Figs. 5-10: • The statistics on the a091m170L and a071m170 ensembles need to be increased to make the largest τ data useful.
• The χ 2 /dof of most fits are reasonable.
• The {4, 2 free } fits have reasonable χ 2 /dof but do not indicate a preference for the small ∆M N π 1 given in Table III. Their ∆M 1 lie closer to or higher than ∆M (2) 1 .
• The ∆M 1 from a two-state fit is expected to be larger since it is an effective combination of the mass gaps of the full tower of excited states. This is illustrated by the difference between ∆M Based on the above arguments, we will choose the {4, 3 * } results obtained after performing the CCFV fits for the final central value. We will also take half the spread in results between the {4 N π , 3 * } and {4, 2 free } strategies, which is ≈ 0.02 in most cases, as a second uncertainty to account for possible unresolved bias from ESC.
The renormalization of the matrix elements is carried out using estimates of Z V D , Z AD , and Z T D calculated on the lattice in the RI −MOM scheme and then converted to the MS scheme at 2 GeV. Two methods to control discretization errors are described in the Appendix B. The final values of Z V D , Z AD , and Z T D used in the analysis are given in Table XII. The values of the three renormalized moments from the seven ensembles and with  the three strategies are summarized in Tables VIII  for renomalization method A and in Table IX for method B. These data are used to perform the CCFV fits discussed next.

VI. CHIRAL, CONTINUUM AND INFINITE VOLUME EXTRAPOLATION
To obtain the final, physical results at M π = 135 MeV, M π L → ∞ and a = 0, we make a simultaneous CCFV fit keeping only the leading correction term in each variable: Note that, since the operators are not O(a) improved in our clover-on-clover formulation, we take the discretization errors to start with a term linear in a.
The fits to the data renormalized using method A for the three strategies are shown in Figs. 1, 2 and 3 and the results are summarized in Table X. The dependence on a is found to be small. The significant variation is with M 2 π , and this is the main discriminant between the three strategies. The smaller the ∆M 1 , the larger is the extrapolation in the ESC fits (difference between the data at the largest τ and τ = ∞ extrapolated value) and a larger slope versus M 2 π in the CCFV fits. The overall consequence for all three moments is that estimates increase by about 0.02 between {4 N π , 3 * } → {4, 3 * } → {4, 2 free }. Based on the observation that the {4, 2 free } fits do not prefer the small ∆M 4 N π 1 , but are closer to ∆M 4 1 (momentum fraction) or larger (helicity and transversity), we take the {4, 3 * } results for our best. However, to account for possible bias due to not having resolved which excited state makes the dominant contribution, we add a second, systematic, error of 0.02 to the final results based on the observed differences in estimates between the three strategies.
The results from the two renormalization methods summarized in Table X overlap-the differences are a fraction of the errors from the rest of the analysis. Also, note that these differences are much smaller than the differences between the Z's from the two methods. This pattern is expected provided the differences in the Z's are largely due to discretization errors that are removed on taking the continuum limit. (52) 0.1166 (54) (17)  Lastly, a comparison between the chiralcontinuum (CC) and CCFV fit results summarized in Table X indicate up to 10% decrease due to the finite volume correction term, however, this is comparable to the size of the final errors. Also, this effect is clear only between the a094m270 and a094m270L data as shown in Tables VIII and IX. Consequently, most of the variation in the CCFV fit occurs for M π L < 4. The result of a CC fit to five larger volume ensembles (excluding a094m270 and a091m170) lie in between the CC and CCFV data shown in Tables VIII and IX. With these caveats, for present, we choose to present final results from the full data set using CCFV fits.
An update of the comparison of lattice QCD calculations on ensembles with dynamical fermions pre-sented in Ref. [16] is shown in the top half of Table XI and in Fig. 4. Our new results, Eq. (20), are consistent with the PNDME 20 values published in Ref. [16]. This is a valuable check since PNDME 20 calculation used the nonunitary cloveron-HISQ lattice formulation. Also, the current calculation provides weak evidence for a finite-volume effect, whereas the PNDME 20 results were obtained using just CC fits. On the other hand, the range of lattice spacings and pion masses simulated is somewhat smaller than in the PNDME 20 calculation.
Our result for the momentum fraction is in very good agreement with estimates from phenomenological global fits reviewed in Ref. [11], summarized in the bottom half of Table XI,

Collaboration
Ref.   [12] and also the comparison in [47] for x u−d . The JAM17 † estimate for x ∆u−∆d is obtained from [11], where, as part of the review, an analysis was carried out using the data in [48].  Using high statistics data we confirmed the behavior of the correlation functions predicted by the spectral decomposition as shown in Figs. 5-10. These higher precision data allowed us to investigate the systematic uncertainty associated with excited-state contamination when extracting the three moments. We carried out the full analysis with three different estimates of the mass gap of the first excited state that cover a large range of possible values. We use the spread in results to assign a systematic uncertainty to account for possible residual excited-state bias.
To obtain the final result in the continuum limit, we fit the seven points using the ansatz in Eq. (19) that includes the leading order terms in M π , the lattice spacing a and the finite volume parameter M π L. Having two pairs of points, {a094m270, a094m270L} and {a091m170, a091m170L}, that differ only in the lattice volume, allowed us to quantify finite volume corrections in all three moments as shown in Figs. 1, 2 and 3. A comparison of the results with and without the finite-volume correction (CCFV versus CC) are shown in Table X. Based on this analysis, we present final results from the CCFV fits that are about 5% smaller than the CC-fit values for the momentum fraction and the helicity moment. These results are consistent with the PNDME 20 [16] calculation, which used the clover-on-HISQ lattice formulation and were obtained using just CC fits.
In appendix B, we describe two methods for removing the p 2 dependent artifacts in the renormalization constants. The results for the moments from these two methods are given in Table X. The data show that after the continuum extrapolation (CCFV or CC fits), the two estimates overlap even though the renormalization constants themselves differ by   Table XI. The left panel compares results for the momentum fraction, the middle for the helicity moment, and the right for the transversity moment. Our NME 20 result is also shown as the blue band to facilitate comparison.
≈ 5% as shown in Table XII. The better agreement after the continuum extrapolation suggests that the main difference between the two methods are indeed discretization artifacts.
The data at three values of the lattice spacing shown in Figs. 1, 2 and 3 do not exhibit any significant dependence on the lattice spacing a. The main variation is with M 2 π , and its magnitude depends on the mass gap of the first excited state used in the analysis of the ESC. Since the mass gaps obtained from fits to the three-point functions ({4, 2 free } strategy) do not prefer values corresponding to the lowest possible excitations (N π states used in the {4 N π , 3 * } strategy) but are closer to two-state fits to the twopoint function (see Table III), we quote final results from the {4, 3 * } strategy. We add a systematic error of 0.02, based on the observed spread (see Table X), to account for possible unresolved excited-state effects.
Our final results, taken from Table X, are given in Eq. (20). These are compared with other lattice calculations and phenomenological global fit estimates in Table XI and Fig. 4. They are in good agreement with other recent lattice results from the PNDME [16], ETMC [38,39], Mainz [32] and χQCD [40] collaborations. Our estimate for the momentum fraction is in good agreement with most global fit estimates but has much larger error. The three estimates for the helicity moment from global fits have a large spread, and our estimate is consis-tent with the smaller error estimates. Lattice estimates for the transversity moment are a prediction.
Having established the efficacy of the lattice QCD approach to reliably calculate these isovector moments, we expect to make steady progress in reducing the errors by simulating at additional values of {a, M π } and by increasing the statistics. In this appendix, we show plots of the unrenormalized isovector momentum fraction, x u−d , the helicity moment, x ∆u−∆d , and the transversity moment, x δu−δd , for the seven ensembles in Figs. 5-10. The data shown is the ratio C 3pt O (τ ; t)/C 2pt (τ ) multiplied by the appropriate factor given in Eqs. (8)- (10) to get the three x . The three panels in each row show fits with the three strategies: {4 N π , 3 * } (left), {4, 3 * } (middle) and {4, 2 free } (right). The fits to C 3pt O (τ ; t) using Eq. (18) are made keeping data at the largest three values of τ , except for a071m170 as discussed in Sec. V. The results of these fit are shown for various τ by lines with the same color as the data. In all cases, to extract the ground state matrix element (blue band), the fits to C 2pt (τ ) and C 3pt O (τ ; t) are done within a single jackknife loop.
The data show a monotonic convergence in τ towards the τ → ∞ estimate. Also, the data are symmetric about t−τ /2 for all values of τ , except for the largest τ on a071m170, a091m170L and a094m270 ensembles, which are statistics limited. Lastly, the largest extrapolation, ie, the difference between the data at t = τ /2 with the largest τ and the τ = ∞ value, is for the {4 N π , 3 * } strategy since it has the smallest mass gap as shown in Table III. This is most evident on the m π ≈ 170 MeV ensembles. The smallest is for the {4, 2 free } strategy in which the mass gap is the largest.

Appendix B: Renormalization
In this appendix, we describe two methods of calculating the renormalization factors, Z V D,AD,T D , for the three one-derivative operators specified in Eqs. (5), (6) and (7). On the lattice, these Z's are first determined nonperturbatively in the RI −MOM scheme [58,59] as a function of the lattice scale p 2 = p µ p µ , and then converted to MS scheme using 3-loop perturbative factors calculated in the continuum in Ref. [60]. For data at each p, we perform horizontal matching by choosing the MS scale µ = |p|. These numbers are then run in the continuum MS scheme from scale µ to 2 GeV using threeloop anomalous dimensions [60]. The two methods differ in how the dependence of Z MS (2GeV) on p 2 a 2 , a lattice artifact, is removed. For details of the three operators and their decomposition into irreducible representations, we refer the reader to Refs. [31,32].
The data for the renormalization factors Z V D,AD,T D in the MS scheme at µ = 2 GeV is shown in Fig. 11 for the seven ensembles as a function of p 2 -the scale of the RI −MOM scheme on the lattice. For all three operators, the data do not show a window in p 2 where the results are independent of p 2 . We analyze the variation with p 2 as being mainly due to a combination of the breaking of full rotational invariance on the lattice and other p 2 dependent artifacts. Many methods have been proposed to control it, see for example Refs. [32,38,61]. We use the following two: • In method A, we take an average over the data points in an interval of 2 GeV 2 about p 2 = Λ/a, where the scale Λ = 3 GeV is chosen to be large enough to avoid nonperturbative effects and at which perturbation theory is expected to be reasonably well behaved. Also, this choice satisfies bothpa → 0 and Λ/p → 0 in the continuum limit as desired. The window over which the data are averaged (given in column three of Table XII) and the error (half the height of the band) are shown by shaded bands in Figs. 11. This method was used in our previous work presented in Ref. [16].
• In method B, we make a fit to the data using the ansatz Z(p) = Z 0 + a µ p µ p µ + b( µ p µ p µ ) 2 to remove the p 2 dependent artifacts. The starting value of p 2 is taken to be the lower limit used in method A, which is given in column three of Table XII, and by which a roughly linear in p 2 behavior is manifest. The results are shown next to the y-axis in Fig. 11 using the star symbol.
These estimates of Z V D , Z AD and Z T D are summarized in Table XII. The discretization errors are expected to be different in the two methods, so we do not average the values of the renormalization constants but perform the full analysis, including the CCFV fits, for the two methods and compare the values of the moments after the continuum extrapolation. These final results are summarized in Table X and found to be consistent.  Table IV along with the values of τ and t skip used in the fit. The y-interval is selected to be the same for all the panels to facilitate comparison.   Table VI. The rest is the same as in Fig. 5.  Figure 8. Continuation of the data and fits to remove excited-state contamination in the extraction of the helicity moment x ∆u−∆d for a091m170 (top row), a091m170L (middle row), a073m270 (third row), and a071m170 (bottom row). The data for the ratio C 3pt O (τ ; t)/C 2pt (τ ) is scaled using Eq. (9) to give x ∆u−∆d , and the fit parameters are listed in Table VI. The rest is the same as in Fig. 5. . Data and fits to remove excited-state contamination in the extraction of the transversity moment x δu−δd for a127m285 (top row), a094m270 (second row), and a094m270L (bottom row) ensembles. The data for the ratio C 3pt O (τ ; t)/C 2pt (τ ) is scaled using Eq. (10) to give x δu−δd , and the fit parameters are listed in Table VII. The rest is the same as in Fig. 5.  Figure 10. Continuation of the data and fits to remove excited-state contamination in the extraction of the transversity moment x δu−δd for a091m170 (top row), a091m170L (middle row), a073m270 (third row), and a071m170 (bottom row). The data for the ratio C 3pt O (τ ; t)/C 2pt (τ ) is scaled using Eq. (10) to give x δu−δd , and the fit parameters are listed in Table VII. The rest is the same as in Fig. 5.  Table XII. Results for the renormalization factors, ZV D,AD,T D , in the MS scheme at 2 GeV. These are calculated in the RI'-MOM scheme as a function of scale p = √ pµpµ on the lattice, matched to the MS scheme at the same scale µ = p, and then run in the continuum MS scheme from µ to 2 GeV. Results are given for two methods used to remove the p 2 dependent artifacts as described in the text. In method A (columns 4-6), the Z's are obtained by averaging the data shown in Fig. 11 over the range of p 2 specified in the in the third column. Results using method B (columns 7-9) are obtained using fits to the data starting with the lower value of p 2 given in column 3 with the ansatz Z(p) = Z0 +ap 2 +bp 4 .  Figure 11. Nonperturbative renormalization factors for x u−d , (ZV D ), x ∆u−∆d , (ZAD), and x δu−δd , (ZT D ) in the MS scheme at µ = 2 GeV for the seven ensembles. The shaded bands mark the region in p 2 that is averaged and the error in the estimate. The points next to the y-axis with the star symbol give a second estimate obtained from the fit Z(p) = Z0 + ap 2 + bp 4 .