Moments of nucleon isovector structure functions in 2+1+1 -flavor QCD
Santanu Mondal, Rajan Gupta, Sungwoo Park, Boram Yoon, Tanmoy Bhattacharya, Huey-Wen Lin
LLA-UR-20-23801MSUHEP-20-009
Moments of nucleon isovector structure functions in -flavor QCD
Santanu Mondal, ∗ Rajan Gupta, † Sungwoo Park, ‡ BoramYoon, § Tanmoy Bhattacharya, ¶ and Huey-Wen Lin
3, 4, ∗∗ (Precision Neutron Decay Matrix Elements (PNDME) Collaboration) Los Alamos National Laboratory, Theoretical Division T-2, Los Alamos, NM 87545 Los Alamos National Laboratory, Computer Computational and Statistical Sciences, CCS-7, Los Alamos, NM 87545 Department of Physics and Astronomy, Michigan State University, MI, 48824, USA Department of Computational Mathematics, Science and Engineering,Michigan State University, East Lansing, MI 48824 (Dated: May 29, 2020)We present results on the isovector momentum fraction, (cid:104) x (cid:105) u − d , helicity moment, (cid:104) x (cid:105) ∆ u − ∆ d ,and the transversity moment, (cid:104) x (cid:105) δu − δd , of the nucleon obtained using nine ensembles of gaugeconfigurations generated by the MILC collaboration using 2 + 1 + 1-flavors of dynamical highlyimproved staggered quarks (HISQ). The correlation functions are calculated using the Wilson-Cloveraction and the renormalization of the three operators is carried out nonperturbatively on the latticein the RI (cid:48) − MOM scheme. The data have been collected at lattice spacings a ≈ . , . , . , and 0.06 fm and M π ≈ ,
220 and 135 MeV, which are used to obtain the physical values usinga simultaneous chiral-continuum-finite-volume fit. The final results, in the MS scheme at 2 GeV,are (cid:104) x (cid:105) u − d = 0 . (cid:104) x (cid:105) ∆ u − ∆ d = 0 . (cid:104) x (cid:105) δu − δd = 0 . PACS numbers: 11.15.Ha, 12.38.GcKeywords: nucleon momentum distribution, helicity, transversity, lattice QCD
I. INTRODUCTION
The elucidation of the hadron structure in termsof quarks and gluons is evolving from determiningthe charges and form factors of nucleons to includingmore complex quantities such as parton distributionfunctions (PDFs) [1], transverse momentum depen-dent PDFs (TMDs) [2], and generalized parton dis-tributions (GPDs) [3] as experiments become moreprecise [4, 5]. These distributions are not measureddirectly in experiments, and phenomenological anal-yses including different theoretical inputs are neededto extract them from experimental data. Input fromlattice QCD is beginning to play an increasinglylarger role in such analyses [6]. In cases where bothlattice results and phenomenological analyses of ex-perimental data (global fits) exist, one can comparethem to validate the control over systematics in thelattice calculations, and on the other hand providea check on the phenomenological process used to ex-tract these observables from experimental data. Inother cases, lattice results are predictions. The list of ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ∗∗ [email protected] quantities for which good agreement between latticecalculations and experimental results, and their pre-cision, has grown very significantly as discussed inthe recent Flavor Averaging Group (FLAG) 2019 re-port [7]. While steady progress has been made in de-veloping the framework for calculating distributionfunctions using lattice QCD [8, 9], even calculationsof their moments have had large statistical and/orsystematic uncertainties prior to 2018. This was thecase even for the best studied quantity, the isovec-tor momentum fraction (cid:104) x (cid:105) u − d [6]. In this work, weshow that the lattice data for the momentum frac-tion, helicity and transversity moments are now ofquality comparable to that for nucleon charges (ze-roth moments). Together with much more precisedata from the planned electron-ion collider [4] andthe Large Hadron Collider, which will significantlyimprove the phenomenological global fits, we antic-ipate steady progress towards a detailed descriptionof the hadron structure.In this paper we present results on the threemoments from high statistics calculations done onnine ensembles generated using 2+1+1-flavors ofHighly Improved Staggered Quarks (HISQ) [10] bythe MILC collaboration [11]. The data at four val-ues of lattice spacings a , three values of the pionmass, M π , including two ensembles at the phys-ical pion mass, and on a range of large physicalvolumes, characterized by M π L , allow us to carry a r X i v : . [ h e p - l a t ] M a y ensemble a M valπ L × T M val π L τ /a aM N N conf N HP N LP ID (fm) (MeV) a m
310 0 . . .
3) 16 ×
48 3 . { , , , , } . ,
668 122 , a m
310 0 . . .
8) 24 ×
64 4 . { , , , } . ,
104 64 , a m
220 0 . . .
9) 32 ×
64 4 . { , , , } . ,
624 73 , a m L . . .
7) 40 ×
64 5 . { , , , } . ,
000 128 , a m
310 0 . . .
8) 32 ×
96 4 . { , , , } . ,
052 144 , a m
220 0 . . .
8) 48 ×
96 4 . { , , , } . ,
668 122 , a m
130 0 . . .
0) 64 ×
96 3 . { , , , } . ,
328 99 , a m W . . .
2) 48 ×
144 4 . { , , , } . − , a m
135 0 . . .
4) 96 ×
192 3 . { , , , } . ,
008 48 , Table I. Lattice parameters, nucleon mass M N , number of configurations analyzed, and the total number of highprecision (HP) and low precision (LP) measurements made. For the a m W ensemble, HP data were not collected,however, we note that the bias correction factor on all other eight ensemble was negligible. out a simultaneous fit in these three variables to ad-dress the associated systematics uncertainties. Wealso investigate the dependence of the results on thespectra of possible excited states included in thefits to remove excited-state contamination (ESC),and assign a second error to account for the asso-ciated systematic uncertainty. Our final results are (cid:104) x (cid:105) u − d = 0 . (cid:104) x (cid:105) ∆ u − ∆ d = 0 . (cid:104) x (cid:105) δu − δd = 0 . Z V D,AD,T D ,for the three operators in Appendix B.
II. LATTICE METHODOLOGY
The parameters of the nine HISQ ensembles aresummarized in Table I. They cover a range of latticespacings (0 . ≤ a ≤ .
15 fm), pion masses (135 ≤ M π ≤ . ≤ M π L ≤ . c sw = 1 /u , where u is the fourth root of the pla-quette expectation value calculated on the hypercu-bic (HYP) smeared [16] HISQ lattices.The masses of light clover quarks were tuned sothat the clover-on-HISQ pion masses, M val π , matchthe HISQ-on-HISQ Goldstone ones, M sea π . M val π val-ues are given in Table I. M sea π values are available inRef. [14]. All fits in M π to study the chiral behav-ior are made using the clover-on-HISQ M val π sincethe correlation functions, and thus the chiral behav-ior of the moments, have a greater sensitivity to it.Henceforth, for brevity, we drop the superscript anddenote the clover-on-HISQ pion mass as M π . Thenumber of high precision (HP) and low precision(LP) measurements made on each configuration inthe truncated solver bias corrected method [17, 18]for cost-effective increase in statistics are specifiedin Table I. III. MOMENTS AND MATRIX ELEMENTS
In this work, we calculate the first moments of spinindependent (or unpolarized), q = q ↑ + q ↓ , helicity(or polarized), ∆ q = q ↑ − q ↓ , transversity, δq = q (cid:62) + q ⊥ distributions, defined as (cid:104) x (cid:105) q = (cid:90) x [ q ( x ) + q ( x )] dx , (1) (cid:104) x (cid:105) ∆ q = (cid:90) x [∆ q ( x ) + ∆ q ( x )] dx , (2) (cid:104) x (cid:105) δq = (cid:90) x [ δq ( x ) + δq ( x )] dx , (3)where q ↑ ( ↓ ) corresponds to quarks with helicityaligned (anti-aligned) with that of a longitudinallypolarized target, and q (cid:62) ( ⊥ ) corresponds to quarkswith spin aligned (anti-aligned) with that of a trans-versely polarized target.These moments, at leading twist, can be extractedfrom the hadron matrix elements of one-derivativevector, axial-vector and tensor operators at zero mo-mentum transfer. The unpolarized and polarizedmoments (cid:104) x (cid:105) q and (cid:104) x (cid:105) ∆ q of the nucleon are alsoobtained from phenomenological global fits while acomputation of the nucleon transversity (cid:104) x (cid:105) δq usingLattice QCD is still a prediction due to lack of suf-ficient experimental data [6].We are interested in extracting the forward nu-cleon matrix elements (cid:104) N ( p ) |O| N ( p ) (cid:105) , with the nu-cleon initial and final momenta, p , taken to be zeroin this work. The complete set of one-derivative vec-tor, axial-vector, and tensor operators are: O µνV a = qγ { µ ←→ D ν } τ a q , O µνA a = qγ { µ ←→ D ν } γ τ a q , O µνρT a = qσ [ µ { ν ] ←→ D ρ } τ a q , (4)where q = { u, d } is the isodoublet of light quarksand σ µν = ( γ µ γ ν − γ ν γ µ ) /
2. The derivative ←→ D ν ≡ ( −→ D ν − ←− D ν ) consists of four terms: ψ (Γ −→ D ν − Γ ←− D ν ) ψ ( x ) ≡ (cid:2) ψ ( x )Γ U ν ( x ) ψ ( x + ν ) − ψ ( x )Γ U † ν ( x − ν ) ψ ( x − ν )+ ψ ( x − ν )Γ U ν ( x − ν ) ψ ( x ) − ψ ( x + ν )Γ U † ν ( x ) ψ ( x ) (cid:3) . (5)Lorentz indices within { } in Eq. (4) are symmetrizedand within [ ] are antisymmetrized. It is also implicitthat, where relevant, the traceless part of the aboveoperators is taken. Their renormalization is carriedout nonperturbatively in the regularization indepen-dent RI (cid:48) -MOM scheme as discussed in Appendix B.A more detailed discussion of these twist-2 operators and their renormalization can be found in Refs. [19]and [20].In this work, we consider only isovector quantities.These are obtained from Eq. (4) by choosing τ a = τ for the Pauli matrix. The decomposition of thematrix elements of these operators in terms of thegeneralized form factors at zero momentum transferis: (cid:104) N ( p, s (cid:48) ) |O µνV a | N ( p, s ) (cid:105) = u pN ( s (cid:48) ) A (0) γ { µ p ν } u pN ( s ) (6) (cid:104) N ( p, s (cid:48) ) |O µνA a | N ( p, s ) (cid:105) = iu pN ( s (cid:48) ) ˜ A (0) γ { µ p ν } γ u pN ( s ) (7) (cid:104) N ( p, s (cid:48) ) |O µνρT a | N ( p, s ) (cid:105) = iu pN ( s (cid:48) ) A T (0) σ [ µ { ν ] p ρ } u pN ( s ) (8)The relation between the momentum fraction, helic-ity moment, and the transversity moment, and thegeneralized form factors is (cid:104) x (cid:105) q = A (0), (cid:104) x (cid:105) ∆ q =˜ A (0) and (cid:104) x (cid:105) δq = A T (0) respectively.We end this discussion by mentioning that otherapproaches have been proposed to calculate themoments of PDFs from Lattice QCD in recentyears [21–23]. IV. CORRELATION FUNCTIONS ANDMOMENTS
We use the following interpolating operator N tocreate/annihilate the nucleon state N = (cid:15) abc (cid:104) q aT ( x ) Cγ (1 ± γ )2 q b ( x ) (cid:105) q c ( x ) , (9)where { a, b, c } are color indices, q , q ∈ { u, d } and C = γ γ is the charge conjugation matrix. Thenonrelativistic projection (1 ± γ ) / state. The zero momentum 2-point and 3-point nu-cleon correlation functions are defined as C ptαβ ( τ ) = (cid:88) x (cid:104) |N α ( τ, x ) N β (0 , ) | (cid:105) (10) C pt O ,αβ ( τ, t ) = (cid:88) x (cid:48) , x (cid:104) |N α ( τ, x ) O ( t, x (cid:48) ) N β (0 , ) | (cid:105) (11)where α , β are spin indices. The source is placed attime slice 0, the sink is at τ and the one-derivativeoperators, defined in Sec. III, are inserted at timeslice t . Data have been accumulated for the valuesof τ specified in Table I, and in each case for allintermediate times 0 ≤ t ≤ τ .To isolate the various operators, projected 2- and3-point functions are constructed as C pt = Tr (cid:0) P pt C pt (cid:1) (12) C pt O = Tr (cid:0) P pt C pt O (cid:1) . (13)The projector P pt = (1 + γ ) in the nucleon cor-relator gives the positive parity contribution for thenucleon propagating in the forward direction. Forthe connected 3-point contributions P pt = (1 + γ )(1 + iγ γ ) is used. With these spin projections,the explicit operators used to calculate the forwardmatrix elements are: (cid:104) x (cid:105) u − d : O V = q ( γ ←→ D − γ · ←→ D ) τ q (14) (cid:104) x (cid:105) ∆ u − ∆ d : O A = qγ { ←→ D } γ τ q (15) (cid:104) x (cid:105) δu − δd : O T = qσ [1 { ←→ D } τ q . (16)Our goal is to obtain the matrix elements ( M E ),of these operators within the ground state of thenucleon. These
M E are related to the moments as: (cid:104) |O V | (cid:105) = − M N (cid:104) x (cid:105) u − d , (17) (cid:104) |O A | (cid:105) = − iM N (cid:104) x (cid:105) ∆ u − ∆ d , (18) (cid:104) |O T | (cid:105) = − iM N (cid:104) x (cid:105) δu − δd , (19)where M N is the nucleon mass. The three momentsare dimensionless, and their extraction on a givenensemble does not require knowing the value of thelattice scale a . It enters only when performing thechiral-continuum extrapolation to the physical pointas discussed in Sec. VI. V. CONTROLLING EXCITED STATECONTAMINATION
To calculate the matrix elements of the opera-tors defined in Sec. III between ground-state nu-cleons, contributions of all possible excited statesneed to be removed. The lattice nucleon interpo-lating operator N given in Eq. (9), however, couplesto the nucleon, all its excitations and multiparticlestates with the same quantum numbers. Previouslattice calculations have shown that the ESC can belarge [24–26]. In our earlier works [12–14, 27], wehave shown that this can be controlled to within afew percent. We use the same strategy here. Inparticular, we use HYP smearing of the gauge linksbefore calculating Wilson-clover quark propagatorswith optimized Gaussian smeared sources using themultigrid algorithm [28, 29]. Correlation functions constructed from these smeared source propagatorshave smaller excited state contamination [27]. To ex-tract the ground state matrix elements from these,we fit the three-point data at several τ -values (listedin Table I) simultaneously using the spectral decom-position given in Eq. (21).Fits to the zero-momentum two-point functions,C , were carried out keeping up to four states inthe spectral decomposition: C pt ( τ ) = (cid:88) i =0 |A i | e − M i τ . (20)Fits are made over a range { τ min − τ max } to ex-tract M i and A i , the masses and the amplitudes forthe creation/annihilation of these states by the in-terpolating operator N . In fits with more than twostates, estimates of the amplitudes A i and masses M i for i ≥ τ min . We used the largest timeinterval allowed by statistics, i.e., by the stabilityof the covariance matrix. We perform two types of4-state fits. In the fit denoted { } , we use the empir-ical Bayesian technique described in the Ref. [30] tostabilize the three excited-state parameters. In thesecond fit, denoted { Nπ } , we use as prior for M either the non-interacting energy of N ( − ) π ( ) orthe N ( ) π ( ) π ( ) state, which are both lower thanthe M obtained from the { } fit, and roughly equalfor the six ensembles. The lower energy N ( − ) π ( )state has been shown to contribute in the axialchannel [31], whereas for the vector channel the N ( ) π ( ) π ( ) state is expected to be the relevantone. We find that these two fits to the two-pointfunction cannot be distinguished on the basis of the χ / DOF, in fact the full range of M between thetwo estimates from { } and { Nπ } are viable first-excited-state masses on the basis of χ / DOF alone.The same is true of the values for M . We there-fore, investigate the dependence of the results formoments on the exited-state spectra by doing thefull analysis with multiple strategies as discussed be-low. The ground-state nucleon mass obtained fromthe various fits is denoted by the common symbol M N ≡ M and the mass gaps by ∆ M i ≡ M i − M i − .The analysis of the zero-momentum three-pointfunctions, C O , is performed retaining upto threestates | i (cid:105) in the spectral decomposition: C O ( τ ; t ) = (cid:88) i,j =0 |A i ||A j |(cid:104) i |O| j (cid:105) e − M i t − M j ( τ − t ) . (21)The operators, O , are defined in Eqs. (14), (15)and (16). By fixing the momentum at the sink tozero and inserting the operator at zero momentum a (fm) x u d / dof = 0.74 a m Wa m a m a m a m a m a m L a m a m M (GeV ) x u d / dof = 0.74 a m Wa m a m a m a m a m a m L a m a m Figure 1. Data for (cid:104) x (cid:105) u − d , renormalized in the MS scheme at µ = 2 GeV, for all nine ensembles. The blue band inthe left panel shows the CC fit result evaluated at M π = 135 MeV and plotted versus a , while in the right panel itshows the result versus M π evaluated at a = 0. a (fm) x u d / dof = 0.56 a m Wa m a m a m a m a m a m L a m a m M (GeV ) x u d / dof = 0.56 a m Wa m a m a m a m a m a m L a m a m Figure 2. Data for (cid:104) x (cid:105) ∆ u − ∆ d , renormalized in the MS scheme at µ = 2 GeV, for all nine ensembles plotted as afunction of a (left panel) and M π (right panel). The rest is the same as in Fig. 1. a (fm) x u d / dof = 0.88 a m Wa m a m a m a m a m a m L a m a m M (GeV ) x u d / dof = 0.88 a m Wa m a m a m a m a m a m L a m a m Figure 3. Data for (cid:104) x (cid:105) δu − δd , renormalized in the MS scheme at µ = 2 GeV, for all nine ensembles plotted as afunction of a (left panel) and M π (right panel). The rest is the same as in Fig. 1. { , ∗ } { , free } ensemble τ /a t skip Observable
M E (cid:104) x (cid:105) χ /dof M E (cid:104) x (cid:105) χ /dof a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m W { , , }∗ { , } (cid:104) x (cid:105) u − d − . . . − . . . a m W { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m W { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m L { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m L { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m L { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . a m { , , }† { , } (cid:104) x (cid:105) u − d − . . . − . . . a m { , , } { , } (cid:104) x (cid:105) ∆ u − ∆ d − . . . − . . . a m { , } { , } (cid:104) x (cid:105) δu − δd − . . . − . . . Table II. Our best estimates of the unrenormalized moments from the two fit strategies, { , ∗ } and { , free } , usedto analyze the two- and three-point functions. The second column gives the values of τ used in the fits and thethird column lists t skip = { i, j } , the number of time-slices from the source and sink omitted for each τ for the twofit types to the three-point functions. For each fit type we give the result for the ground state matrix element, ME ,the moment (cid:104) x (cid:105) obtained from it using Eqs. (17)–(19), and the χ /DOF of the fit to the three-point function. In twocases, the values of τ /a included are different: the ∗ in the second column denotes τ /a = { , , } and † denotes τ /a = { , } were used for the { , free } fits. (cid:104) x (cid:105) u − d Ensemble fit-type a ∆ M a ∆ M (cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) χ /dof a m { , } . . . .
34) 0 . . a m { Nπ , } . . . .
81) 0 . . a m { , ∗ } . . . . .
66) 0 . − . − . . a m { Nπ , ∗ } . . . . .
0) 0 . − . . .
2) 1 . a m { , free } . . . .
44) 0 . . a m { , } . . . .
5) 1 . .
31) 0 . a m { Nπ , } . a m { , ∗ } . . . . .
9) 0 . . − . a m { Nπ , ∗ } . . . . .
8) 0 . . − . .
2) 0 . a m { , free } . . . . Table III. Comparison of fits using five strategies, { , } , { Nπ , } , { , ∗ } , { Nπ , ∗ } and { , free } , for the momentumfraction (cid:104) x (cid:105) u − d on two ensembles a m
310 (highest statistics and M π ∼
310 MeV) and a m
135 (physical M π ∼
135 MeV). In the { , free } fit, the excited state mass gap, ∆ M , is left as a free parameter that is determined fromthe fit to the three-point function. The values of τ /a and t skip used are same as listed in Table II. We could not finda { Nπ , } fit to the a m
135 data that gave reasonable values. (cid:104) x (cid:105) ∆ u − ∆ d Ensemble fit-type a ∆ M a ∆ M (cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) χ /dof a m { , } . . . .
6) 0 . . a m { Nπ , } . . . .
5) 0 . . a m { , ∗ } . . . . .
6) 0 . . − . a m { Nπ , ∗ } . . . . .
3) 0 . . . a m { , free } . . . . a m { , } . . . .
09) 1 . . a m { Nπ , } . . − . a m { , ∗ } . . . . .
84) 0 . . − . a m { Nπ , ∗ } . . . − . . − . .
5) 1 . a m { , free } . . . . Table IV. Comparison of fits using five strategies, { , } , { Nπ , } , { , ∗ } , { Nπ , ∗ } and { , free } , for the helicitymoment (cid:104) x (cid:105) ∆ u − ∆ d . The rest is the same as in Table III. (cid:104) x (cid:105) δu − δd Ensemble fit-type a ∆ M a ∆ M (cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) (cid:104) |O| (cid:105)(cid:104) |O| (cid:105) χ /dof a m { , } . . . .
1) 0 . . a m { Nπ , } . . − . .
9) 1 . . a m { , ∗ } . . . . .
0) 0 . . . a m { Nπ , ∗ } . . . − . .
8) 0 . . . a m { , free } . . . . a m { , } . . . .
75) 2 . . a m { Nπ , } . a m { , ∗ } . . . . .
2) 0 . . − . a m { Nπ , ∗ } . . . − . . − . .
7) 1 . a m { , free } . . . . Table V. Comparison of fits using five strategies, { , } , { Nπ , } , { , ∗ } , { Nπ , ∗ } and { , free } , for the transversitymoment (cid:104) x (cid:105) δu − δd . The rest is the same as in Table III. transfer we get the forward matrix element. Thepractical challenge discussed above is determiningthe relevant M and M to use, and failing that, toinvestigate the sensitivity of (cid:104) |O| (cid:105) to possible val-ues of M and M and including that variation as asystematic uncertainty.For a given strategy for determining M and M ,we extract the desired ground state matrix ele-ment (cid:104) |O| (cid:105) by fitting the three-point correlators C O ( t ; τ ) for a subset of values of t and τ simulta-neously. This subset is chosen to reduce ESC—weselect the largest values of τ and discard t skip num-ber of points next to the source and sink for each τ .These values of τ and of t skip are given in Table II.The data for the ratio C O ( τ ; t ) /C ( τ ) areshown in Figs. 5 and 6 in the Appendix A for allnine ensembles. The signal in the three-point corre-lators decreases somewhat from momentum fractionto helicity moment to transversity moment. Never-theless, we are able to make 3 ∗ state (3-state with (cid:104) |O| (cid:105) = 0) fits in all cases. The spectral decompo-sition predicts that the data for all three quantities issymmetric about t = τ /
2, however, on some of theensembles, and for some of the larger values of τ ,the data show some asymmetry, which is indicativeof the size of statistical fluctuations that are present.The fits to C ( τ ) and C O ( τ ; t ) are carried outwithin a single-elimination jackknife process, whichis used to get both the central values and the errors.We have investigated five fit types, { , } , { Nπ , } , { , ∗ } , { Nπ , ∗ } and { , free } , based onthe spectral decomposition to understand and con-trol ESC. The labels { m, n } denote an m -state fitto the two-point function and an n -state fit to thethree-point function. In the 2 free -fit to the three-point function, M is left as a free parameter, whilea 3 ∗ -fit is a 3-state fit with (cid:104) | O | (cid:105) = 0. The re-sults from the five strategies for the momentum frac-tion, (cid:104) x (cid:105) u − d , in Table III, for the helicity moment, (cid:104) x (cid:105) ∆ u − ∆ d , in Table IV, and for the transversity mo-ment, (cid:104) x (cid:105) δu − δd , in Table V illustrate the observedbehavior for the a m
310 ensemble, which has thehighest statistics, and the physical mass ensemble a m
135 at the smallest value of a .For all three observables, the five results in Ta-bles III, and IV, and V for the ground state ma-trix element, (cid:104) | O | (cid:105) , are consistent within 2 σ onthe a m
310 ensemble. On the a m
135 ensem-ble, the difference in ∆ M ≡ M − M between { } and { Nπ } analyses becomes roughly a factorof two, and ∆ M from the { free } fit is larger thaneven the { } value, i.e., the { free } fit does not pre-fer the small ∆ M given by { Nπ } . On the otherhand, the ∆ M from a two-state fit is expected tobe larger since it is an effective combination of themass gaps of the full tower of excited states. Due to a small ∆ M , fits with the spectrum from { Nπ } failon a m { , ∗ } and { , free } fits gives results consistent within 2 σ .The estimates from these two fit types are given inTable II. To summarize, our overall strategy is tokeep as many excited states as possible without over-parameterization of the fits. We, therefore, choose,for the central values, the { , ∗ } results, and to takeinto account the spread due to the fit type, we adda second, systematic, uncertainty to the final resultsin Table VII. This is taken to be the difference be-tween the results obtained by doing the full analysiswith the { , ∗ } and { , free } strategies.The renormalization of the matrix elements is car-ried out using estimates of Z V D , Z AD , and Z T D calculated on the lattice in the RI (cid:48) − MOM schemeand then converted to the MS scheme at 2 GeV asdescribed in the Appendix B. The final values of Z V D , Z AD , and Z T D used in the analysis are givenin Table IX.
VI. CHIRAL, CONTINUUM AND INFINITEVOLUME EXTRAPOLATION
To obtain the final, physical results at M π =135 MeV, M π L → ∞ and a = 0, we make a simulta-neous CCFV fit keeping only the leading correctionterm in each variable: (cid:104) x (cid:105) ( M π ; a ; L ) = c + c a + c M π + c M π e − M π L √ M π L . (22)Note that, since the operators are not O ( a ) improvedand we used the Clover-on-HISQ formulation, wetake the discretization errors to start with a termlinear in a . The fits to the { , ∗ } data from thenine ensembles are shown in Figs. 1, 2 and 3. Thefit parameters are summarized in Table VI.The results of the CCFV fits show that the finitevolume correction term, c , is not constrained at all.We therefore, also present results from a CC fit with c = 0 in Eq. (22). Results for c from the two fitansatz overlap and there is a small positive slope inboth a and M π for all three quantities. The data forboth { , ∗ } and { , free } , given in Table II, are verysimilar, but with a systematic shift of about 0.01–0.02 in all three cases. This difference arises because∆ M for { , free } is larger (except in a m τ is fromabove as shown in Figs. 5 and 6, i.e., a larger ∆ M implies a smaller extrapolation and a larger τ → ∞ value.For our final results we quote the CC fit values asthe coefficient c of the finite-volume corrections infit-type Observable c c c c χ /dofCC (cid:104) x (cid:105) u − d . . . . (cid:104) x (cid:105) ∆ u − ∆ d . . . . (cid:104) x (cid:105) δu − δd . . . . (cid:104) x (cid:105) u − d . . . − . (cid:104) x (cid:105) ∆ u − ∆ d . . . − . (cid:104) x (cid:105) δu − δd . . . . Table VI. Results for the fit parameters in the CCFV ansatz given in Eq. (22) and used for the chiral, continuum andfinite volume (CCFV) extrapolation of the { , ∗ } data. The CC and CCFV fit-types correspond to fits with c = 0or c (cid:54) = 0. Observable { , ∗ } { , free } Best Estimate (cid:104) x (cid:105) MS u − d (cid:104) x (cid:105) MS∆ u − ∆ d (cid:104) x (cid:105) MS δu − δd Table VII. Results for the three moments from the twostrategies { , ∗ } and { , free } . For our best estimates,we take the { , ∗ } values and assign a second, system-atic, error that is the difference between the two results.The results are in the MS scheme at scale 2 GeV. the CCFV fits is undetermined. The CC results withthe two strategies, { , ∗ } and { , free } , are summa-rized in Table VII. For our best estimates, we takethe { , ∗ } results and add a second, systematic, er-ror that is the difference between these two strate-gies and represents the uncertainty in controlling theexcited-state contamination.A comparison of these results with other lat-tice QCD calculations on ensembles with dynam-ical fermions is presented in the top half of Ta-ble VIII. Our results agree with those from theMainz group [20] that have also been obtained usingdata on a comparable number of ensembles, but allwith M π >
200 MeV, to perform a chiral and contin-uum extrapolation. The one difference is the slope c of the chiral correction. For our clover-on-HISQformulation, we find a small positive value while theMainz data show a small negative value [20]. Our re-sults are also consistent within 1 σ with the ETMC20 [32] and ETMC 19 [33] values that are from asingle physical mass ensemble. Results for momen-tum fraction and the helicity moment from RQCD18 [34] are taken from their Set A with the differ-ence between Set A and B values quoted as a secondsystematic uncertainty. The result for the transver-sity moment is from a single 150 MeV ensemble.These values are larger, especially for the helicity and transversity moment. Other earlier lattice re-sults show a spread, however, in each of these cal-culations, the systematics listed in the last columnof Table VIII have not been addressed or controlledand could, therefore, account for the differences.Estimates from phenomenological global fits, mostof which have also been reviewed in Ref. [6], are sum-marized in the bottom of Table VIII. We find that re-sults for the momentum fraction from global fits are,in most cases, 1–2 σ smaller and have much smallererrors. Results for the helicity moment are consis-tent and the size of the errors comparable. Latticeestimates of the transversity moment are a predic-tion. VII. CONCLUSIONS
In this paper, we have presented results for theisovector quark momentum fraction, (cid:104) x (cid:105) MS u − d , helic-ity moment, (cid:104) x (cid:105) MS∆ u − ∆ d , and transversity moment, (cid:104) x (cid:105) MS δu − δd , from a high statistics lattice QCD calcula-tion. Attention has been paid to the systematic un-certainty associated with excited-state contamina-tion. We have carried out the full analysis with dif-ferent estimates of the mass gaps of possible excitedstates, and use the difference in results between thetwo strategies that give stable fits on all ensemblesas an additional systematic uncertainty to accountfor possible residual excited-state contamination.The behavior versus M π , the lattice spacing a and finite volume parameter M π L have been inves-tigated using a simultaneous fit that includes theleading correction in all three variables as given inEq. (22). The nine data points cover the range0 . < a < .
15 fm, 135 < M π <
320 MeV and3 . < M π L < .
5. Over this range, all three mo-ments, (cid:104) x (cid:105) MS u − d , (cid:104) x (cid:105) MS∆ u − ∆ d and (cid:104) x (cid:105) MS δu − δd , do not showa large dependence on a or M π or M π L . As shown0Collaboration Ref. (cid:104) x (cid:105) u − d (cid:104) x (cid:105) ∆ u − ∆ d (cid:104) x (cid:105) δu − δd RemarksPNDME 20 0 . . . N f = 2 + 1 + 1(this work) clover-on-HISQETMC 20 [32] 0 . N f = 2 + 1 + 1 Twisted MassN-DIS, N-FVETMC 19 [33] 0 . . . N f = 2 + 1 + 1 Twisted MassN-DIS, N-FVMainz 19 [20] 0 . stat . stat . stat N f = 2 + 1 Clover(+14 , − sys (+10 , − sys (+16 , − sys RQCD 18 [34] 0 . . . N f = 2 CloverETMC 17 [35] 0 . N f = 2 Twisted MassN-DIS, N-FVETMC 15 [36] 0 . . . N f = 2 Twisted MassN-DIS, N-FVRQCD 14 [25] 0 . N f = 2 CloverN-DIS, N-CE, N-FVLHPC 14 [37] 0 . N f = 2 + 1 CloverN-DIS ( a ∼ .
12 fm)RBC/ [38] 0.124–0.237 0.146–0.279 N f = 2 + 1 Domain WallUKQCD 10 N-DIS, N-CE, N-ESLHPC 10 [39] 0 . . N f = 2 + 1Domain Wall-on-AsqtadN-DIS, N-CE, N-NR, N-ESCT18 [40] 0 . † [6, 41] 0.241(26)NNPDF3.1 [42] 0 . . . . . . . . Table VIII. Our Lattice QCD results are compared with other lattice calculations with N f flavors of dynamicalfermions in rows 2–9, and with results from phenomenological global fits in the remainder of the table. In bothcases, the results are arranged in reverse chronological order. All results are in the MS scheme at scale 2 GeV. Fora discussion and comparison of lattice and global fit results up to 2017, see Ref. [6] and a more recent comparisonin [40] for (cid:104) x (cid:105) u − d . The JAM17 † estimate for (cid:104) x (cid:105) ∆ u − ∆ d is obtained from [6], where, as part of the review, an analysiswas carried out using the data in [41]. The following abbreviations are used in the remarks column for various sourcesof systematic uncertainties in lattice calculations—DIS: Discretization effects, CE: Chiral extrapolation, FV: Finitevolume effects, NR: Nonperturbative renormalization, ES: Excited state contaminations. A prefix “N-” means thatthe systematic uncertainty was not adequately controlled or not estimated. G l oba l F i t La tt i c e Q CD MMHT2014 CT14 HERAPDF2.0CJ15 ABMP2016 NNPDF3.1 CT18 LHPC 10 RBC/UKQCD 10LHPC 14 RQCD 14 ETMC 15 ETMC 17 ETMC 19 ETMC 20 RQCD 18 Mainz 19 PNDME 20 〈 x 〉 u − d G l oba l F i t La tt i c e Q CD DSSV08 NNPDFpol1.1JAM17 LHPC 10 RBC/UKQCD 10ETMC 15 ETMC 19 RQCD 18 Mainz 19 PNDME 20 〈 x 〉 ∆ u − ∆ d La tt i c e Q CD ETMC 15 ETMC 19 RQCD 18 Mainz 19PNDME 20 〈 x 〉 δ u − δ d Figure 4. A comparison of results from lattice QCD calculations with dynamical fermions and global fits (below theblack line) summarized in Table VIII. The left panel compares results for the momentum fraction, the middle for thehelicity moment, and the right for the transversity moment. The PNDME 20 result is also shown as the blue bandto facilitate comparison.
10 5 0 5 10 t a m
135 : x = 0.155(14), /34 = 0.87 = =22 =20 =18 =16
10 5 0 5 10 t a m
135 : x = 0.191(12), /34 = 1.00 = =22 =20 =18 =16
10 5 0 5 10 t a m
135 : x = 0.185(16), /34 = 1.32 = =22 =20 =18 =16
10 5 0 5 10 t a m W : x = 0.170(13), /28 = 1.02 = =24 =22 =20 =18
10 5 0 5 10 t a m W : x = 0.223(15), /28 = 1.00 = =24 =22 =20 =18
10 5 0 5 10 t a m W : x = 0.213(18), /28 = 0.80 = =24 =22 =20 =18 t a m
130 : x = 0.177(8), /22 = 0.93 = =16 =14 =12 =10 t a m
130 : x = 0.218(7), /22 = 1.30 = =16 =14 =12 =10 t a m
130 : x = 0.212(11), /22 = 1.30 = =16 =14 =12 =10 t a m
220 : x = 0.184(5), /22 = 0.89 = =16 =14 =12 =10 t a m
220 : x = 0.227(4), /22 = 0.92 = =16 =14 =12 =10 t a m
220 : x = 0.234(6), /22 = 1.29 = =16 =14 =12 =10 Figure 5. Data and fits for a m
135 (top row), a m W (second row), a m
130 (third row) and a m
220 (lastrow). In each row, the three panels shows the ratio C O ( τ ; t ) /C ( τ ) scaled according to Eq. (17)–(19) to give (cid:104) x (cid:105) u − d (left), (cid:104) x (cid:105) ∆ u − ∆ d (middle), and (cid:104) x (cid:105) δu − δd (right). For each τ , the line in the same color as the data points isthe result of the fit used (see Sec. V) to obtain the ground state matrix element. The ensemble ID, the final result (cid:104) x (cid:105) (also shown by the blue band and summarized in Table II), the values of τ , and χ / DOF of the fit are also givenin the legends. The interval of the y-axis is selected to be the same for all the panels to facilitate comparison. t a m x =0.196(4), /22=1.25 = =16 =14 =12 =10 t a m x =0.233(3), /22=1.24 = =16 =14 =12 =10 t a m x =0.239(4), /22=0.78 = =16 =14 =12 =10 t a m x =0.199(8), /16=1.32 = =14 =12 =10 =8 t a m x =0.234(9), /16=0.92 = =14 =12 =10 =8 t a m x =0.246(11), /16=1.24 = =14 =12 =10 =8 t a m L : x =0.191(9), /16=1.44 = =14 =12 =10 =8 t a m L : x =0.240(7), /16=1.33 = =14 =12 =10 =8 t a m L : x =0.237(10), /16=1.25 = =14 =12 =10 =8 t a m x =0.195(11), /16=1.66 = =14 =12 =10 =8 t a m x =0.238(16), /16=0.76 = =14 =12 =10 =8 t a m x =0.251(16), /16=0.69 = =14 =12 =10 =8 t a m x =0.214(6), /10=1.94 = =9 =8 =7 =6 t a m x =0.266(7), /10=0.76 = =9 =8 =7 =6 t a m x =0.293(12), /6=0.66 = =9 =8 Figure 6. Data and fits for a m
310 (top row), a m
220 (second row), a m L (third row), a m
310 (fourth row)and a m
310 (bottom row) ensembles. The rest is the same as in Fig. 5. p p (GeV) Z a fm Z VD Z AD Z TD p p (GeV) Z a fm Z VD Z AD Z TD p p (GeV) Z a fm Z VD Z AD Z TD p p (GeV) Z a fm Z VD Z AD Z TD Figure 7. Nonperturbative renormalization factors for (cid:104) x (cid:105) u − d , ( Z V D ), (cid:104) x (cid:105) ∆ u − ∆ d , ( Z AD ), and (cid:104) x (cid:105) δu − δd , ( Z TD ), at thefour lattice spacings in the MS scheme at µ = 2 GeV. The shaded bands mark the region in (cid:112) p that is averagedand the error in the estimate. in Table VI, possible dependence on the lattice size,characterized by M π L , is not resolved by the data,i.e., the coefficient c is unconstrained. We, there-fore, take for our final results those obtained fromjust the chiral-continuum fit.The small increase with a and M π , evident inFigs 1–3, is well fit by the leading correction termsthat are linear in these variables. Also, for all threeobservables, the chirally extrapolated value is con-sistent with the data from the two physical massensembles. In short, the observed small dependencein all three variables, and having two data pointsat M π ∼
135 MeV to anchor the chiral fit, al-lows a controlled extrapolation to the physical point, M π = 135 MeV and a = 0.Our final results, given in Table VII, are comparedwith other lattice calculations and phenomenolog-ical global fit estimates in Table VIII. Estimates of all three quantities are in good agreement withthose from the Mainz collaboration [20], also ob-tained using a chiral-continuum extrapolation, andthe ETMC collaboration [32, 33] that are from asingle physical mass ensemble. On the other hand,most global fit estimates for the momentum fractionare about 10% smaller and have much smaller er-rors, while those for the helicity moment are consis-tent within 1 σ . Lattice estimates for the transversitymoment are a prediction. The overall consistencyof results suggests that lattice QCD calculations ofthese isovector moments are now mature and futurecalculations will steadily reduce the statistical andsystematic uncertainties in them.5 p p (GeV) Z V D a m a m a m a m Figure 8. Nonperturbative renormalization factor Z V D for (cid:104) x (cid:105) u − d is calculated on four ensembles, one from eachlattice spacing. The shaded bands give the interval in (cid:112) p over which the data are averaged to get the result andthe error in the estimate. ACKNOWLEDGMENTS
We thank the MILC collaboration for sharing theHISQ ensembles, and Martha Constantinou, Gian-nis Koutsou, Emanuele Nocera and Juan Rojo fordiscussions. The calculations used the Chroma soft-ware suite [51]. Simulations were carried out on com-puter facilities of (i) the National Energy ResearchScientific Computing Center, a DOE Office of Sci-ence User Facility supported by the Office of Sci-ence of the U.S. Department of Energy under Con-tract No. DE-AC02-05CH11231; and, (ii) the OakRidge Leadership Computing Facility at the OakRidge National Laboratory, which is supported bythe Office of Science of the U.S. Department of En-ergy under Contract No. DE-AC05-00OR22725; (iii)the USQCD Collaboration, which are funded by theOffice of Science of the U.S. Department of Energy,and (iv) Institutional Computing at Los Alamos Na-tional Laboratory. T. Bhattacharya and R. Guptawere partly supported by the U.S. Department ofEnergy, Office of Science, Office of High EnergyPhysics under Contract No. DE-AC52-06NA25396. T. Bhattacharya, R. Gupta, S. Mondal, S. Park andB.Yoon were partly supported by the LANL LDRDprogram. The work of H.-W. Lin is partly supportedby the US National Science Foundation under grantPHY 1653405 “CAREER: Constraining Parton Dis-tribution Functions for New-Physics Searches”, andby the Research Corporation for Science Advance-ment through the Cottrell Scholar Award.
Appendix A: Plots of the Ratio C pt O ( τ ; t ) /C pt ( τ ) In this Appendix, we show in Figs. 5 and 6, plotsof the unrenormalized isovector momentum fraction, (cid:104) x (cid:105) u − d , the helicity moment, (cid:104) x (cid:105) ∆ u − ∆ d , and thetransversity moment, (cid:104) x (cid:105) δu − δd , for the nine ensem-bles. The data show the ratio C O ( τ ; t ) /C ( τ )multiplied by the appropriate factor given inEqs. (17)–(19) to get (cid:104) x (cid:105) . The lines with thesame color as the data are the result of the fit to C O ( τ ; t ) using Eq. (21). In all cases, to extract theground state matrix element, the fits to C ( τ ) and C O ( τ ; t ) are done within a single jackknife loop.6 Appendix B: Renormalization a fit range Z V D Z AD Z T D (fm) (GeV )0 .
06 9 . − . . . . .
09 5 . − . . . . .
12 4 . − . . . . .
15 2 . − . . . . Table IX. Results for the renormalization factors, Z V D,AD,TD , in the MS scheme at 2 GeV. These are cal-culated in the RI (cid:48) − MOM scheme as a function of scale p = √ p µ p µ on the lattice, matched to the MS scheme atthe same scale µ = p , and then run in the continuum MSscheme from µ to 2 GeV. The results are the average ofvalues over the range of | p | specified in the second col-umn. The final error estimate is taken to be twice thatshown in Figs. 7 and 8. In this appendix, we describe the calculationof the renormalization factors, Z V D,AD,T D , forthe three one-derivative operators. These are de-termined nonperturbatively on the lattice in theRI (cid:48) − MOM scheme [52, 53] as a function of thelattice scale p = p µ p µ , and then converted to MSscheme using 3-loop perturbative factors calculatedin the continuum in Ref. [54]. For data at each p , we perform horizontal matching by choosing theMS scale µ = | p | . These numbers are then run inthe continuum MS scheme from scale µ to 2 GeVusing three-loop anomalous dimensions [54]. Thiscalculation of Z V D,AD,T D is done for one value of M π at each a . Based on our experience with local operators [13], where we found insignificant depen-dence of results on M π , we assume that these values,within errors, give the mass-independent renormal-ization factors at each a . Also, the decomposition ofthe three operators into irreducible representationsgiven in Refs. [19, 20], show that they can only mixwith higher dimensional operators. In our CCFVfits, such O ( a ) effects would also be taken into ac-count and removed by the continuum extrapolation.In Fig. 7, we show the behavior of the renormaliza-tion factors Z V D,AD,T D in the MS scheme at µ = 2GeV for the four ensembles as a function of | p | —thescale of the RI (cid:48) − MOM scheme on the lattice. InFig. 8 we compare Z V D used to renormalize (cid:104) x (cid:105) u − d for the four ensembles, one at each lattice spacing.For all three operators, the data do not show awindow in | p | where the results are independent of | p | . Thus, for the final estimates we use the strategylabeled “Method B” in Ref. [13]. This correspondsto taking an average over the data points in an inter-val of 2 GeV about p = Λ /a , where the scale Λ = 3GeV is chosen to be large enough to avoid nonper-turbative effects and at which perturbation theory isexpected to be reasonably well behaved. This choicesatisfies both pa → /p → p , we take twice the error (full height of the band)for the error estimate for all three Z (cid:48) s in the finalanalysis.These final estimates of Z V D , Z AD and Z T D usedto renormalize the momentum fraction, the helicitymoment and the transversity moment, respectively,are given in Table IX. [1] R. Brock et al. (CTEQ), Rev. Mod. Phys. , 157(1995).[2] B. Yoon, M. Engelhardt, R. Gupta, T. Bhat-tacharya, J. R. Green, B. U. Musch, J. W. Negele,A. V. Pochinsky, A. Schafer, and S. N. Syritsyn,Phys. Rev. D , 094508 (2017), arXiv:1706.03406[hep-lat].[3] M. Diehl, Generalized parton distributions , Ph.D.thesis (2003), arXiv:hep-ph/0307382.[4] A. Accardi et al. , Eur. Phys. J. A , 268 (2016),arXiv:1212.1701 [nucl-ex].[5] D. Boer et al. , (2011), arXiv:1108.1713 [nucl-th].[6] H.-W. Lin et al. , Prog. Part. Nucl. Phys. , 107(2018), arXiv:1711.07916 [hep-ph].[7] S. Aoki et al. (Flavour Lattice Averaging Group),Eur. Phys. J. C , 113 (2020), arXiv:1902.08191[hep-lat].[8] K. Cichy and M. Constantinou, Adv. High Energy Phys. , 3036904 (2019), arXiv:1811.07248 [hep-lat].[9] N. Karthik, “Lattice computations of PDF:Challenges and progress,” https://indico.cern.ch/event/764552/contributions/3420535/attachments/1864018/3064443/LatticeTalk19_NK.pdf (2019), accessed: 2020-05-20.[10] E. Follana et al. (HPQCD Collaboration, UKQCDCollaboration), Phys.Rev. D75 , 054502 (2007),arXiv:hep-lat/0610092 [hep-lat].[11] A. Bazavov et al. (MILC Collaboration), Phys.Rev.
D87 , 054505 (2013), arXiv:1212.4768 [hep-lat].[12] 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].[13] T. Bhattacharya, V. Cirigliano, S. Cohen, R. Gupta,H.-W. Lin, and B. Yoon, Phys. Rev.
D94 , 054508 (2016), arXiv:1606.07049 [hep-lat].[14] R. Gupta, Y.-C. Jang, B. Yoon, H.-W. Lin,V. Cirigliano, and T. Bhattacharya, Phys. Rev. D98 , 034503 (2018), arXiv:1806.09006 [hep-lat].[15] B. Sheikholeslami and R. Wohlert, Nucl. Phys.
B259 , 572 (1985).[16] A. Hasenfratz and F. Knechtli, Phys.Rev.
D64 ,034504 (2001), arXiv:hep-lat/0103029 [hep-lat].[17] G. S. Bali, S. Collins, and A. Schafer,Comput.Phys.Commun. , 1570 (2010),arXiv:0910.3970 [hep-lat].[18] T. Blum, T. Izubuchi, and E. Shintani, Phys.Rev.
D88 , 094503 (2013), arXiv:1208.4349 [hep-lat].[19] M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt,P. E. L. Rakow, G. Schierholz, and A. Schiller,Phys. Rev.
D53 , 2317 (1996), arXiv:hep-lat/9508004 [hep-lat].[20] T. Harris, G. von Hippel, P. Junnarkar, H. B. Meyer,K. Ottnad, J. Wilhelm, H. Wittig, and L. Wrang,Phys. Rev.
D100 , 034513 (2019), arXiv:1905.01291[hep-lat].[21] Z. Davoudi and M. J. Savage, Phys. Rev.
D86 ,054505 (2012), arXiv:1204.4146 [hep-lat].[22] W. Detmold and C. J. D. Lin, Phys. Rev.
D73 ,014501 (2006), arXiv:hep-lat/0507007 [hep-lat].[23] W. Detmold, I. Kanamori, C. J. D. Lin, S. Mon-dal, and Y. Zhao,
Proceedings, 36th InternationalSymposium on Lattice Field Theory (Lattice 2018):East Lansing, MI, United States, July 22-28, 2018 ,PoS
LATTICE2018 , 106 (2018), arXiv:1810.12194[hep-lat].[24] 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].[25] G. S. Bali, S. Collins, B. Glle, M. Gckeler, J. Najjar,R. H. Rdl, A. Schfer, R. W. Schiel, A. Sternbeck,and W. Sldner, Phys. Rev.
D90 , 074510 (2014),arXiv:1408.6850 [hep-lat].[26] G. S. Bali, S. Collins, M. Deka, B. Glassle, M. Gock-eler, J. Najjar, A. Nobile, D. Pleiter, A. Schafer,and A. Sternbeck, Phys. Rev.
D86 , 054504 (2012),arXiv:1207.1110 [hep-lat].[27] B. Yoon et al. , Phys. Rev.
D93 , 114506 (2016),arXiv:1602.07737 [hep-lat].[28] R. Babich, J. Brannick, R. Brower, M. Clark,T. Manteuffel, et al. , Phys.Rev.Lett. , 201602(2010), arXiv:1005.3043 [hep-lat].[29] M. Clark, R. Babich, K. Barros, R. Brower,and C. Rebbi, Comput.Phys.Commun. , 1517(2010), arXiv:0911.3191 [hep-lat].[30] B. Yoon et al. , Phys. Rev. D , 074508 (2017),arXiv:1611.07452 [hep-lat].[31] Y.-C. Jang, R. Gupta, B. Yoon, and T. Bhat-tacharya, Phys. Rev. Lett. , 072002 (2020),arXiv:1905.06470 [hep-lat].[32] C. Alexandrou, S. Bacchio, M. Constantinou,J. Finkenrath, K. Hadjiyiannakou, K. Jansen,G. Koutsou, H. Panagopoulos, and G. Spanoudes,(2020), arXiv:2003.08486 [hep-lat].[33] C. Alexandrou et al. , Phys. Rev. D , 034519(2020), arXiv:1908.10706 [hep-lat].[34] G. S. Bali, S. Collins, M. Gckeler, R. Rdl, A. Schfer, and A. Sternbeck, Phys. Rev. D , 014507 (2019),arXiv:1812.08256 [hep-lat].[35] C. Alexandrou, M. Constantinou, K. Hadjiyian-nakou, K. Jansen, C. Kallidonis, G. Koutsou, A. Va-quero Avils-Casco, and C. Wiese, Phys. Rev. Lett. , 142002 (2017), arXiv:1706.02973 [hep-lat].[36] A. Abdel-Rehim et al. , Phys. Rev. D92 , 114513 (2015), [Erratum: Phys.Rev.D93,no.3,039904(2016)], arXiv:1507.04936[hep-lat].[37] J. R. Green, M. Engelhardt, S. Krieg, J. W. Negele,A. V. Pochinsky, and S. N. Syritsyn, Phys. Lett.
B734 , 290 (2014), arXiv:1209.1687 [hep-lat].[38] Y. Aoki, T. Blum, H.-W. Lin, S. Ohta,S. Sasaki, et al. , Phys.Rev.
D82 , 014501 (2010),arXiv:1003.3387 [hep-lat].[39] J. Bratt et al. (LHPC Collaboration), Phys.Rev.
D82 , 094502 (2010), arXiv:1001.3620 [hep-lat].[40] T.-J. Hou et al. , (2019), arXiv:1912.10053 [hep-ph].[41] J. Ethier, N. Sato, and W. Melnitchouk, Phys. Rev.Lett. , 132001 (2017), arXiv:1705.05889 [hep-ph].[42] R. D. Ball et al. (NNPDF), Eur. Phys. J.
C77 , 663(2017), arXiv:1706.00428 [hep-ph].[43] S. Alekhin, J. Blmlein, S. Moch, and R. Placakyte,Phys. Rev.
D96 , 014011 (2017), arXiv:1701.05838[hep-ph].[44] A. Accardi, L. T. Brady, W. Melnitchouk, J. F.Owens, and N. Sato, Phys. Rev.
D93 , 114017(2016), arXiv:1602.03154 [hep-ph].[45] H. Abramowicz et al. (H1, ZEUS), Eur. Phys. J.
C75 , 580 (2015), arXiv:1506.06042 [hep-ex].[46] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston,P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump,and C. P. Yuan, Phys. Rev.
D93 , 033006 (2016),arXiv:1506.07443 [hep-ph].[47] L. A. Harland-Lang, A. D. Martin, P. Motylinski,and R. S. Thorne, Eur. Phys. J.
C75 , 204 (2015),arXiv:1412.3989 [hep-ph].[48] E. R. Nocera, R. D. Ball, S. Forte, G. Ridolfi, andJ. Rojo (NNPDF), Nucl. Phys.
B887 , 276 (2014),arXiv:1406.5539 [hep-ph].[49] D. de Florian, R. Sassot, M. Stratmann, andW. Vogelsang, Phys. Rev.
D80 , 034030 (2009),arXiv:0904.3821 [hep-ph].[50] D. de Florian, R. Sassot, M. Stratmann, andW. Vogelsang, Phys. Rev. Lett. , 072001 (2008),arXiv:0804.0422 [hep-ph].[51] R. G. Edwards and B. Joo (SciDAC, LHPC,UKQCD), Nucl. Phys. B Proc. Suppl. , 832(2005), arXiv:hep-lat/0409003.[52] M. Gockeler et al. , Phys. Rev.
D82 , 114511(2010), [Erratum: Phys. Rev.D86,099903(2012)],arXiv:1003.5756 [hep-lat].[53] M. Constantinou, M. Costa, M. Gckeler, R. Hors-ley, H. Panagopoulos, H. Perlt, P. E. L. Rakow,G. Schierholz, and A. Schiller, Phys. Rev.
D87 ,096019 (2013), arXiv:1303.6776 [hep-lat].[54] J. A. Gracey, Nucl. Phys.