Nucleon Helicity Parton Distribution Function in the Continuum Limit with Self-Renormalization

We present the first lattice calculation of the nucleon isovector helicity parton distribution function (PDF) in the framework of large-momentum effective theory (LaMET) that uses the hybrid scheme with self-renormalization. We use ensembles generated by the MILC collaboration at lattice spacings $a=\{0.1207,0.0888,0.0582\}$ fm, with $N_f=2+1+1$ flavors of highly improved staggered quarks at sea pion mass of $M_{\pi}\approx 315$ MeV. We use clover-improved action for our valence quarks with nucleon boost momentum $P_z\approx 1.75$ GeV and high-statistics measurements for the LaMET matrix elements. We perform an extrapolation to the continuum limit and improve the handling of systematic errors using renormalization-group resummation (RGR) and leading-renormalon resummation (LRR). Our final nucleon helicity PDF is renormalized in the $\overline{\text{MS}}$ scheme at energy scale $\mu=2.0$ GeV. We compare our results with and without the two systematic improvements of RGR and LRR at each lattice spacing as well as the continuum limit, and we see that the application of RGR and LRR greatly reduces the systematic errors across the whole $x$ range. Our continuum results with both RGR and LRR show a small positive antiquark region for the nucleon helicity PDF as well as a change of as much as a factor of two in the central values compared to results with neither RGR or LRR. By contrast, the application of RGR and LRR only changes the central values by about 5\% in the quark region. We compare our lattice results with the global fits by the JAM, NNPDF and DSSV collaborations, and we observe some tension between our results.


Introduction
Parton distribution functions (PDFs) describe non perturbatively the probability distribution of specific longitudinal momentum fractions, x, of a hadron's constituent quarks and gluons.Among them, helicity PDFs provide information on the difference between the parton having its spin aligned and opposite to the hadron's spin.Experimentally, great progress has been made in the study of nucleon helicity PDFs through semi-inclusive deep-inelastic scattering and high-energy muon scattering experiments, such as by the HERMES [1] and COMPASS [2,3] collaborations at DESY and CERN, respectively.Other experiments, such as the STAR [4,5,6,7] and PHENIX [8,9,10] collaborations which study proton-proton collisions at the Relativistic Heavy Ion Collider (RHIC) are able to probe the antiquark contribution to the helicity PDF.Future experiments such as those that have been proposed at the Electron Ion Collider (EIC) [11,12,13,14,15] in the United States as well as the Electron Ion Collider in China (EicC) [16] and the AMBER experiment at CERN [17] will probe the spin structure of the nucleon with even greater accuracy and help refine the global fits of the helicity PDFs.
Progress in global fits has been made in extracting helicity PDFs from experimental data by the DSSV [18], NNPDF [19] and JAM [20] collaborations.DSSV [18] showed that there is a positive antiquark region for the proton helicity PDF (∆u(x) − ∆d(x) > 0) as well as that quarks and antiquarks only contribute up to one third of the proton's spin.NNPDF [19] extended the proton helicity PDF calculation to the small-x region for polarized gluons as well as performed a calculation of the individual unpolarized light quark and antiquark PDFs.JAM [20] determined that the contribution to the proton's spin from strange quarks is compatible with zero, which provided more information about the total spin content of the proton.However, assumptions often have to be made in the global analysis of helicity PDFs, such as the antiquark region being zero, the functional form of the PDF, or exact SU(3) flavor symmetry.Both the DSSV [18] and NNPDF [19] collaborations made the assumption ∆s = ∆s in their calculations.JAM relaxed this assumption by including semi-inclusive deep inelastic scattering (SIDIS) data [20] but found that any deviation from zero was unable to rise above background noise.While global fitting has the advantage of using experimental data from real events, it does not perform a direct PDF calculation from the first principles of quantum chromodynamics (QCD) and is limited by the quality of experimental data that is available.
A direct calculation of PDFs became possible with the advent of large-momentum effective theory (LaMET) in 2013 [21,22,23], which allowed parton physics to be studied on the Euclidean lattice.By studying the behavior of spatially separated correlators boosted to large momentum, parton physics can then be recovered through perturbative matching.The bare matrix elements extracted from the lattice are h B (z, a) = P z ψ − z 2 ΓW − z 2 , z 2 ψ z 2 P z , where |P z ⟩ is a hadron state at boost momentum P z in the z direction, ) is the Wilson line connecting the two spacetime coordinates (0, 0, −z/2, 0) and (0, 0, z/2, 0), ψ is the fermion field, and Γ is a Dirac structure which for helicity matrix elements is Γ = γ z γ 5 .To extract the ground-state matrix element, we use a two-state fit on the two-point correlators and a two-sim fit on the three-point correlators, following, for example, Ref. [24].The extraction of the matrix element can be visualized by plotting the ratio of the threepoint and two-point correlation functions at multiple values of source-sink separation t sep .The ratio is observed to approach the computed ground-state matrix element as t sep increases.In addition, we vary the minimum and maximum values of t sep (denoted by t min sep and t max sep ), again following the example of Ref. [24].We find that the ground-state matrix elements are compatible for these different values.This shows that excited-state contamination is under control and not a significant source of systematic error.The LaMET method has been used for a numerical determination of the nucleon helicity PDF in Refs.[25,26,27,28,29,30,31,32,33].The very first helicity PDF calculation was done in 2016 [25] using 310 MeV pion mass at lattice spacing a = 0.12 fm with boost momenta P z = {0.43,0.86, 1.29} GeV.Mass correction and one-loop matching was used to better align the behavior of the quasi-PDF with that of the lightcone.This was the first time the antiquark section had been studied and showed an asymmetry u(x) > d(x).Soon after, there were many followup LaMET calculations by different collaborations with different fermion actions and even physical pion masses [27,28,29,30,31,32,33,34]. Most LaMET-method helicity PDFs have been renormalized nonperturbatively only in the RI/MOM and similar renormalization schemes [27,28,30,31,32,33].The latest hybrid-or selfrenormalization schemes [35,36] have not been applied to the helicity PDF yet.
In this work, we make the first calculation of the isovector nucleon helicity PDF using matrix elements renormalized in the hybrid scheme with self-renormalization (HSR).The hybrid scheme was introduced in Ref. [35] and the selfrenormalization in Ref. [36].We use this method to address two sources of divergence in the matrix elements: the linear divergence which arises from the self-energy in the Wilson line of the bare matrix elements and the renormalon divergence which arises from the asymptotic nature of the perturbation series.This is done by fitting the matrix elements to a functional form derived from perturbative QCD which has the added advantage of accounting for discretization effects.It extracts the renormalization factor and the remaining non-perturbative physics directly from the matrix elements.It also reduces the dependence on lattice spacing a during the renormalization process, which allows for a reliable extrapolation to the continuum limit.These two sources of divergence can also be addressed using the hybrid scheme introduced in Ref. [35].However, this method does not account directly for discretization effects, and the removal of the linear divergence in this scheme is a very delicate ex-ercise and a significant source of systematic error, as shown in Ref. [36].The HSR method has been used in the calculation of nucleon transversity PDFs [37] as well as the calculation of pion and kaon distribution amplitudes [38,39].In both cases, the HSR process greatly reduces the dependence on lattice spacing compared with pure RI/MOM scheme.
We can further improve our calculation by augmenting both the renormalization scheme and the lightcone matching with the methods of renormalization-group resummation (RGR) [38,40,41] and leading-renormalon resummation (LRR) [41].RGR involves resumming the large logarithms that arise from the difference in renormalization scale and the intrinsic physical scale.The method involves setting the renormalization scale such that the logarithms vanish and then evolving to the desired energy scale using the renormalization-group equation (RGE).The perturbation series also contains a renormalon divergence [42], which is enhanced by the use of RGR on its own.We account for this effect by including LRR in the calculations which resums the series to account for the renormalon divergence.We present the first application of these methods to the nucleon helicity PDF.

Self-Renormalization in Hybrid Scheme
In this work, we use clover lattice fermion action for the valence quarks on top of 2+1+1 flavors (degenerate up and down quarks plus strange and charm quarks at their physical masses in the QCD vacuum) of hypercubic (HYP)-smeared [43] HISQ [44,45], generated by MILC Collaboration.The lattice parameters include lattice spacings a ∈ [0.06, 0.12] fm, pion mass M π ≈ 310 MeV and M π L ≈ 4.5.The quark masses for the clover action are tuned to reproduce the lightest sea HISQ pseudoscalar meson masses, and the clover parameters are set to the tree-level tadpole-improved values.We carefully monitor for signatures of exceptional configurations due to the non-unitary lattice-QCD formulation of mixed-action approach and found exceptional configurations to be absent for these three MILC ensembles [44].This mixed-action setup is the same as the one used in works done by PNDME Collaboration in many studies of nucleon structure [46,47,48,49].
On the lattice, we calculate the time-independent, nonlocal matrix elements stated in Sec. 1.We use Gaussian momentum smearing [50] for the quark field.Such a momentum source is designed to increase the overlap with nucleons of the desired boost momentum, and we are able to reach higher boost momentum for the nucleon states.We use multigrid algorithm [51,52] in the Chroma software package [53] to speed up the inversion of the quark propagator for the clover fermions.To make sure excited-state contamination is under control, we measure at least four nucleon three-point source-sink separations, and we perform simultaneously two-state extraction of ground-state nucleon matrix elements.Details of our calculation parameters can be found in Table 1.
To renormalize the nucleon matrix elements for helicity PDFs, the renormalization factors in the RI/MOM scheme [54] are calculated in Ref. [55].The real renormalization factors as a function of z are shown in the leftmost panel of Fig. 1  the imaginary ones are consistent with zero for p R z = 0 with lattice spacings a ≈ {0.12, 0.09, 0.06} fm denoted by red circles, green squares and blue triangles, respectively.We can see that the RI/MOM factor decays more quickly as a decreases across the whole z range; this is due to the linear divergence becoming more severe as lattice spacing decreases.The same behavior was also observed by ETMC [28], in which the inverse renormalization factors in the RI ′ -MOM scheme for the helicity operator were calculated as a function of z at lattice spacings 0.09, 0.08 and 0.06 fm.This linear divergence can be quantified as the kz/a ln aΛ QCD term mentioned in Eqs. 1 and 3 in Sec. 2, which becomes larger as a decreases at fixed z.Later in this work, we will remove this divergence with the HSR procedure.In Fig. 1, we show the matrix elements renormalized in the pure-RI/MOM scheme with lattice spacings a ≈ {0.12, 0.09, 0.06} fm shown as red circles, green squares and blue triangles, respectively, with the real (imaginary) part in the middle (rightmost) panel.Our renormalized matrix elements are also normalized to 1 at z = 0.This is equivalent to dividing the renormalized matrix elements by the axial charge g A which can be determined both on the lattice and experimentally.
LaMET calculations use matrix elements boosted to large momentum, and we renormalize them in HSR [36].The method involves fitting bare matrix elements at P z = 0 to a functional form dictated by perturbative QCD and the fitting parameters determined in this process can then be used to divide out the sources of divergence in the large-momentum matrix elements relevant for LaMET.The case of NLO can be supplemented with the systematic improvements of RGR and LRR, which we denote by NLR.The renormalization factor Z(z, a) fitted to the functional form can either be an RI/MOM factor Z RI (z, a, µ R , p R z = 0) or a hadron matrix element evaluated at zero boost momentum.The two functional forms are equally valid, as was demonstrated in Ref. [56], where the hybrid-RI/MOM and hybrid-ratio schemes were compared.The m 0 parameter was extracted in the same way for both schemes by demanding the short-distance behavior (z ≲ 0.3 fm) agree with perturbation theory.It was demonstrated that the renormalized matrix elements are compatible well within error bars between these two schemes.This shows that we can choose either an RI/MOM factor or a zero-boost-momentum matrix element without introducing additional systematic uncertainties.In our calculation, we choose the RI/MOM factor.In addition, the matching kernels coincide at NLO for the hybrid-ratio scheme and the hybrid-RI/MOM scheme for RI/MOM matrix elements evaluated at zero-momentum, as was shown in Ref. [57].We follow the procedure of Refs.[36,37,39] for NLO and of Ref. [38] for NLR.
We describe, first, the self-renormalization procedure at NLO.The logarithm of the renormalization factors is fitted to the functional form [36] where z is the Wilson length; k is the linear divergence arising from the self-energy of the Wilson link.The linear divergence contains a renormalon ambiguity [36] which is accounted for later in the HSR procedure.a is the lattice spacing; Λ QCD is the cutoff scale for QCD; f 1 (z) is a function describing the discretization effects of the lattice.The ensembles generated in Ref. [44] use clover-improved action for the valence quarks and highly improved staggered quarks (HISQ) for the sea, so the discretization terms are O(a).µ is the final desired renormalization scale; d is a parameter determined by demanding that the short-distance behavior (z ≲ 0.3 fm as suggested by Ref. [35]) of the renormalized matrix element agrees with perturbation theory; C F is the quadratic Casimir for the fundamental representation of SU(3) and β 0 is the first coefficient of the QCD beta-function.The expression g 1 (z) is the residual term that describes the non perturbative physics after the divergences and discretization effects are removed.We, thus, have an expression for the bare matrix element from which the different sources of divergence can be divided out.We interpolate the renormalization factors of Wilson length z ∈ [0, 1.20] fm in uniform steps of 0.06 fm for all three ensembles.Next, we determine the linear divergence and QCD scale (k and Λ QCD ) by setting the parameter d to zero initially and fitting to Eq. 1 for multiple pairs of k and Λ QCD .For this work, we choose Λ QCD = 0.2 GeV at minimum χ 2 = 1.04 at k = 0.795 when d = 0.The next step is to remove the logarithmic divergence and renormalon divergence by setting the short-distance renormalized matrix elements to agree with the Wilson coefficient C NLO 0 (z, µ).The helicity Wilson coefficients have been computed to NLO [58,59]: + 7 , where α s (µ) is the strong coupling at energy scale µ, and γ E is the Euler-Mascheroni constant.Since we tabulated the RI/MOM factors in steps of z = 0.06 fm, we fit to perturbation theory in the interval z ∈ [0.06, 0.18] fm.We demand that g 1,2 (z) − ln C NLO 0 (z, µ) = m 0 z + c (a linear function) whose y-intercept is less than 10 −3 to ensure good matching to perturbation theory.We tune the d parameter to match ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲   C NLO 0 (z, µ) and we find the value d = 0.497 gives a y-intercept of O(10 −5 ).We now have a minimum χ 2 = 1.02 at k = 0.798, which is a change of less than 0.5% in the value of k compared to d = 0. We then construct a full renormalization factor To perform HSR at NLR, we first modify the fitting function to [38] where the µ dependence has been removed since it will be resummed with the RGE and we include the conversion constant ∆I.We use the same values for k, Λ QCD and d from the NLO case as was done in Ref. [60] since the first two are global parameters that depend on the bare matrix elements regardless of whether we use RGR or LRR.This time, the conversion constant is tuned so as to match the short-distance behavior with C NLR 0 (z, µ) Wilson coefficient which is defined in Ref. [41].We find ∆I = −0.330has a corresponding y-intercept of O(10 −3 ).The full self-renormalization factor for NLR is We perform the continuum extrapolation on the quantity Z X self (z, a)h B (z, a), where X is NLO or NLR and finally convert to the hybrid scheme with self-renormalization with the factor where z s ≈ 0.3 fm is the maximum distance at which perturbation theory is valid.Thus, the full hybrid-renormalized matrix element is h R,X HSR (z, a) = Z X HSR (z, a)h B (z, a).
The matrix elements renormalized in the selfrenormalization scheme, h R,X self (z), for both X being NLO and NLR are shown in Fig. 2 to demonstrate the effectiveness of self-renormalization.We show the renormalization factors at a ≈ {0.12, 0.09, 0.06} fm in red circles, green squares and blue triangles, respectively, and the corresponding Wilson coefficient as a black-dotted line.Note that the dependence on lattice spacing a is almost completely removed by the selfrenormalization process.The results at the smallest and largest lattice spacings differ by no more than 6% for z ∈ [0, 1.0] fm.In addition, the matrix elements agree with the corresponding Wilson coefficient for short distances, plotted as a black dashed line.The systematic errors are estimated using the method of "scale variation" as in Refs.[41,60].When we use the RGR for the matrix elements, we set the initial scale µ = z −1 so as to eliminate the logarithms and then evolve to the final desired energy scale µ = 2.0 GeV.Scale variation involves setting the initial scale to c ′ × z −1 for c ′ = 0.75 and c ′ = 1.5, as was used in Ref. [41]; the central value corresponds to c ′ = 1.0.This corresponds to a variation of approximately 15% on either side of α s (µ = 2.0 GeV) in the strong coupling.
The first thing we notice is that there is a significant decrease in systematic errors when going from NLO to NLR for all lattice spacings in the real and imaginary parts.The relative systematic errors for the real part decrease by up to a factor of eight for a = 0.1207 fm and a = 0.0582 fm, and up to a factor of seventeen at a = 0.0888 fm.The effect on the imaginary part is even greater with a decrease of as much as a factor of twenty for all lattice spacings.
The continuum extrapolation is performed on the quantity h R,X self (λ, a) by fitting to a linear function h R,X self (λ, a) = c(λ) × a + h R,X self (λ, a = 0) where λ = zP z is Ioffe time.In the continuum case, the relative systematic errors for the real part decrease by as much as a factor of six when NLO is supplemented with both RGR and LRR; the same quantities for the imaginary part decrease by as much as a factor of eight.This is the same effect observed in Ref. [60], in which the systematic errors are drastically reduced with the RGR and LRR improvements, as well as the fact that the large logarithms and renormalon divergence are significant sources of systematic errors.In Fig. 3 we show the renormalized matrix elements at fixed X and variable a.The top (bottom) row shows the hybrid-renormalized matrix elements at NLO (NLR) while the left (right) column shows the real (imaginary) part with those from lattice spacings a = {0.1207,0.0888, 0.0582} fm shown in red circles, green squares and blue points and continuum-extrapolated ones shown as a purple band, respectively.We show both statistical errors (solid inner lines) and statistical and systematic errors combined with quadrature (dashed outer lines).We also see that the relative systematic errors vary by as much as a factor of two between the largest and smallest lattice spacings for both NLR and NLO in the real and imaginary parts.The scale of the systematics at NLO and NLR affect the systematics of the corresponding continuum extrapolation.The smaller systematic errors in the a 0 data at NLR compared to NLO results in smaller systematics in the continuum extrapolation.Comparing both NLO and NLR renormalized matrix elements from a = 0.0888 fm and 0.0582 fm (the results with the smallest statistical errors), we see that the real parts differ by at most one sigma across the full range of λ.The imaginary parts show some tension in the region λ ∈ [4,6] but elsewhere are also compatible within two sigma.By contrast, the matrix elements renormalized in the pure-RI/MOM scheme shown in Fig. 1 show compatibility for both real and imaginary parts within two sigma but have much larger error bars than the matrix elements renormalized with the HSR.In the continuum limit, the real part of the central values in the region λ ∈ [0, 4] differ by up to 10% with those at a = 0.0582 fm and up to 20% with those at a = 0.1207 fm.The imaginary part of the central values in the same λ range differ by no more than 5% with those at a = 0.0582 fm and 7% with those at a = 0.1207 fm.Both of these demonstrate good convergence in the continuum limit.

Helicity Parton Distribution Function
To obtain the nucleon isovector helicity PDF, we first need to Fourier transform the HSR matrix elements h R,X HSR (z, a) to momentum space to obtain the quasi-PDF ∆ qX (x, P z ) = ∞ −∞ P z dz 2π e ixzP z h R,X HSR (z, a).However, simply truncating the integral at the limit of the Wilson length extent on the lattice will cause unphysical oscillations in the quasi-PDF [35,61,62].To prevent this, we extrapolate the renormalized matrix elements to infinite distance.The large-distance behavior of the renormalized matrix element determines the small-x behavior of the lightcone PDF; the extrapolation model for infinite distance is derived from a model assumption x k 1 , commonly used in the global-fit community, for naive small-x region of PDFs.We adopt the large-z extrapolation model previously used in Refs.[35,60,61,63]: h R,X HSR (z, a)/g A → Ae −mz |zP z | n as z → ∞, where A, m and n are fitting parameters.The range of z values used must be large enough to realistically capture the long-distance behavior but not so large that the matrix elements are too noisy.In all cases, the χ 2 /dof value is less than 1 which demonstrates that the large distance behavior is well captured by the fitting model.We can then use combined renormalized matrix elements directly from the lattice and those large-z extrapolation to Fourier transform to obtain quasi-PDFs.
The final stage in the calculation is the perturbative matching which is used to align the UV behavior of the quasi-PDF with the lightcone.The lightcone PDF is related to the quasi-PDF via ∆q where K −1,X is the matching kernel for case X, which has been derived up to one-loop for helicity PDFs renormalized in the hybrid scheme [57,58] with finite-P z corrections [61,64].To make the matching kernel at NLR, we first add the LRR modification term ∆K LRR defined in Refs.[38,41] to the matching kernel to make K NLO+LRR = K NLO + ∆K LRR .The process of RGR in the matching follows the same philosophy as in the case of renormalization.The matching is performed with the kernel K NLO+LRR at energy scale µ = 2xP z (such that the scale dependence vanishes from matching formula).We then evolve the matched PDF to the final desired energy scale (in our case, µ = 2.0 GeV), this time using the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation: d∆q X (x,µ) dln(µ 2 ) = 1 x dz |z| P(z)∆q X x z , µ , where P(z) is the DGLAP kernel, which has been calculated up to three loops [65].It should be noted that the DGLAP evolution formula begins to break down at |x| ≲ 0.2 where α s (µ = 2xP z ) becomes nonperturbative.We use the same algorithm for RGR matching as was detailed in Appendix D Ref. [40].
We compare the lattice-spacing dependence of NLO (left) and NLR (right) helicity PDFs in Fig. 4 with a = {0.1207,0.0888, 0.0582, 0.0} fm bands shown in solid red, solid green, horizontally-hatched blue and vertically-hatched purple, respectively.We can see that at the NLO antiquark region has the greatest deviation from zero at a = 0.1207 fm, the largest lattice spacing, but tends towards zero as lattice spacing decreases, including in the continuum limit.All NLO antiquark distributions are compatible with zero with full error bars for x = [−1, −0.15].By contrast, in the NLR case, the antiquark region does not show dependence on lattice spacing and shows a small positive value.In the quark region, the relative systematic errors of the NLO helicity PDF are not correlated with lattice spacing: the largest systematic error occurs at a = 0.0888 fm in the interval x = [0.3,0.5] and at a = 0.0582 fm in the interval x = [0.5, 0.85].For x = [0.15,0.4] the smallest systematic error occurs at the largest lattice spacing.However, with NLR PDF, both the upper and lower systematic errors decrease with lattice spacing in the range x = [0.3,0.7].The correlation disappears at small x and large x where the RGR matching and LaMET calculations begin to break down, and the results become unreliable.This suggests that the discretization effects are the dominant source of systematics in the NLR case.We also compare the PDFs at a = 0.1207 fm and a = 0.0582 fm (the largest and smallest lattice spacings): in the interval x = [0.2,0.8], the central values differ by no more than 12% for both NLO and NLR.The lattice-spacing dependence of the helicity PDF by ETMC [28] with RI ′ -MOM and RI-xMOM renormalization schemes shows a greater difference across a narrower range of lattice spacings: a ∈ {0.06, 0.08, 0.09} fm.Our own relatively small variation between our largest and smallest lattice spacings compared to the above shows that more of the lattice-spacing dependence is removed by the HSR procedure compared to RI ′ -MOM and RI-xMOM.self (z, a) at lattice spacings a ≈ {0.12, 0.09, 0.06} fm are shown as red circles, green squares and blue triangles, respectively; the left (right) plot shows NLO (NLR).All but the largest a have been offset slightly to the right from their true z value to allow for readability.We also plot the Wilson coefficient C X 0 (z, µ) as a dotted black line, which agrees with the renormalized data for short distances.The purple band shows the term exp(g(z) − m 0 z).
▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲  Figure 3: Real (left column) and imaginary (right column) matrix elements renormalized in the HSR at lattice spacings a = {0.1207,0.0888, 0.0582} fm and the continuum limit plotted in red, green, blue and purple, respectively.Data are shown at both NLO (top row) and NLR (bottom row).The continuum extrapolation is computed by performing a weighted fit to a function linear in a of the discrete data at fixed z-value.The solid error bars are statistical and the dotted error bars are combined statistical and systematic errors the latter of which are derived from scale variation as described earlier in this section.

NLR
In Fig. 5 we compare our continuum-limit NLR x-dependent PDFs (green) with the global fits of the NNPDFpol1.1 [19] (red), JAM [20] (cyan) and DSSV [18] (blue) collaborations with the quark-region (antiquark-region) shown in the left (right) panel.Up until now, our results have been divided by the axial charge g A .The global fits do not use normalization so we multiply our helicity PDFs by the value g A = 1.218(25)(30) computed in Ref. [46].We use this value since the calculation was performed on the same lattices as our LaMET calculation.In the quark region, we note that there is tension between the global fits of JAM'17 and each of DSSV'08 and NNPDFpol1.1 from mid to large x, which suggests that there are other sources of systematic errors that may have been ignored or underestimated, likely due to the difference in the experimental cuts, theory inputs, parametrization, and so on.For example, JAM excludes SIDIS data; DSSV and NNPDFpol1.1 rely on assumptions such as SU(3) symmetry to constrain the analysis and add a very small symmetry-breaking term.These assump-tions are needed due to the difficulties in constraining data from polarized experiments.Future experiments with neutral-and charged-current DIS (such as at the EIC) will provide useful measurements to constrain our understanding of the antiquark helicity distribution.Our lattice PDF result also has significant tension with each of the global fits at large x and do not become compatible with zero as x → 1.This is in contrast to the global fits, which use a (1 − x) b term in the parametrization form, enforcing that the PDF goes to zero as x → 1.Similar behavior was observed in the past LaMET lattice calculations of the nucleon helicity PDF [25,28,29] at heavy quark mass where no parametrization form is used.Turning to the antiquark region, we reflect our x-dependent PDFs in the vertical to allow a direct comparison with global fits.Our NLR result for antiquark helicity favors more polarized up than down flavor (whereas NLO does not).This agrees with the results of the STAR [6] and PHENIX [10] collaborations which measured ∆u(x) > ∆d(x) and the global fits who use these experimental data as inputs.The total antiquark flavor asymmetry from this work is 1 0.2 d x (∆u(x) − ∆d(x)) = 0.037 +0.019 −0.023 .

Conclusion
In this paper we perform the first LaMET calculation of the isovector nucleon helicity PDF using the hybrid scheme with self renormalization and the first lattice-QCD helicity PDF results with the RGR and LRR improvements in both the renormalization and the lightcone matching process.We use lattice spacings a = {0.1207,0.0888, 0.0582} fm and perform a continuum extrapolation with a pion mass of M π ≈ 315 MeV.We demonstrate the use of the hybrid-self renormalization scheme in its application to the isovector nucleon helicity PDF.We show the fine-tuning procedure for d and ∆I to match the renormalized matrix elements to perturbation theory at z ≲ 0.2 fm for NLO and NLR, respectively.
Our matrix elements renormalized in the hybrid-self scheme show a reduced dependence on lattice spacing compared to pure-RI/MOM allowing for a more accurate extrapolation to the continuum.Both renormalized matrix elements and quasi-PDFs have their systematic errors due to scale variation reduced by the LRR and RGR procedures by as much as a factor of 10, but the central values by only about 10%.We compared the central values of the NLO and NLR x-dependent helicity PDFs and found that that the application of RGR and LRR has a minor impact on the PDF in the quark region but a significant one in the antiquark region, where it changes sign.We then use the NLR PDF to compare with global fits, since it is more reliable, due to its accounting for large logarithms and the renormalon divergence.Our results in quark region do show some tension with the global fits of DSSV'08 [18], JAM'17 [20] and NNPDFpol1.1 [19], but so are the JAM'17 and NNPDFpol1.1,for example.Our antiquark NLR PDF shows agreement with STAR [6] and PHENIX [10] results in that u(x) > d(x) and contribute to asymmetry of 1 0.2 d x (∆u(x) − ∆d(x)) = 0.037 +0.019 −0.023 .Our calculation can be supplemented with future improvements in lattice-QCD calculation systematics for the LaMET method, such as a higher boost momentum and a physical pion mass.   ) at NLR (green) with global fits computed in Refs.[18] ("DSSV'08", blue), [20] ("JAM'17", cyan) and [19] ("NNPDFpol1.1",red).The left (right) plot shows the quark (antiquark) region.Our results use the axial charge g A = 1.218(25)(30) from Ref. [46].

Figure 1 :
Figure 1: RI/MOM factors (leftmost figure) at lattice spacings a ≈ {0.12, 0.09, 0.06} fm denoted by red circles, green squares and blue triangles, respectively, as a function of Wilson length, z; the imaginary parts are compatible with zero.The middle (rightmost) figure shows the real (imaginary) matrix elements renormalized in the pure RI/MOM scheme at the same lattice spacings and the same color scheme as the leftmost plot, as functions of Ioffe time λ = zP z .

Figure 2
Figure 2: h Rself (z, a) at lattice spacings a ≈ {0.12, 0.09, 0.06} fm are shown as red circles, green squares and blue triangles, respectively; the left (right) plot shows NLO (NLR).All but the largest a have been offset slightly to the right from their true z value to allow for readability.We also plot the Wilson coefficient C X 0 (z, µ) as a dotted black line, which agrees with the renormalized data for short distances.The purple band shows the term exp(g(z) − m 0 z).

Figure 4 :
Figure 4: Lightcone PDFs at variable lattice spacing with pion mass M π = 310 MeV renormalized in the HSR scheme at NLO (left) and NLR (right).PDFs from lattice spacings a = {0.1207,0.0888, 0.0582} fm are plotted in solid-red, solid-green and horizontally hatched blue, respectively.The continuum extrapolation is plotted in vertically hatched purple.The inner error bars are statistical, and the outer error bars are combined statistical and systematic errors.

Table 1 :
Ensemble information and parameters used in this work.N meas is the total number of measurements of the three-point correlators for different values of t sep .L indicates the spatial length which is aN s (in fm).