Parton Distribution Function with Non-perturbative Renormalization from Lattice QCD

We present lattice results for the isovector unpolarized parton distribution with nonperturbative RI/MOM-scheme renormalization on the lattice. In the framework of large-momentum effective field theory (LaMET), the full Bjorken-$x$ dependence of a momentum-dependent quasi-distribution is calculated on the lattice and matched to the ordinary lightcone parton distribution at one-loop order, with power corrections included. The important step of RI/MOM renormalization that connects the lattice and continuum matrix elements is detailed in this paper. A few consequences of the results are also addressed here.


I. INTRODUCTION
Parton distribution functions (PDFs) are probability densities of quarks and gluons seen by an observer moving at the speed of light relative to the hadron. They are universal nonperturbative properties of the hadron. In a global analysis the hard-scattering cross sections can be factorized into the PDFs and the short-distance matrix elements calculable in perturbation theory. Once known, the PDFs can be used as inputs to predict cross sections in high-energy scattering experiments, one of the major successes of QCD. Today, multiple collaborations provide regular updates concerning the phenomenological determination of the PDFs [1][2][3][4][5][6] using the latest experimental results, either focusing on medium-energy QCD experiments or high-energy ones such as those at the LHC. After the past half century of theoretical and experimental efforts, the precision needed in PDFs to further test the Standard Model has increased significantly, and experiments are planned to push further into unexplored or less known regions, such as sea-quark and gluonic structure.
In this work, we continue the first-principles calculation of the PDFs using lattice QCD. In parton physics, the PDFs are defined as the nucleon matrix elements of quark or gluon correlation operators along the lightcone direction. For example, the unpolarized quark distribution is defined as where µ is the renormalization scale in a given renormalization scheme, such as the MS scheme, the nucleon momentum P µ = (P 0 , 0, 0, P z ), ξ ± = (t ± z)/ √ 2 are the lightcone coordinates, and the Wilson line is inserted to ensure gauge invariance. The main obstacle in directly computing PDFs from lattice QCD is their lightcone dependence. Since lattice QCD is formulated in Euclidean space which maps the whole Minkowski lightcone to a single point, the ξ − dependence is completely lost. Early lattice studies of PDFs used the operator product expansion (OPE) to calculate their moments, which are matrix elements of local gauge-invariant operators [7][8][9][10]. However, discretization error and operator mixing due to the breaking of rotational symmetry on the lattice make it hard to go beyond the first few moments. There exist proposals for obtaining higher moments by using smeared sources [11] or computing current-current correlators in Euclidean space [12][13][14][15]. Recently, there is a lattice study on the feasibility of computing pion DA with Euclidean current correlators [16]. Recently, Ji [17,18] proposed a new approach for the direct computation of parton physics, large-momentum effective field theory (LaMET). According to this approach, in order to get the normal PDF, one can start by calculating a "quasi-PDF", which is defined as a spatial correlation of partons along, say the z direction, in a moving nucleon, h(z, P z ,μ) = P |O Γ |P , arXiv:1706.01295v4 [hep-lat] 23 Jan 2018 where O Γ =ψ(z)ΓW z (z, 0)ψ(0) with Γ = γ z or Pz Pt γ t [19][20][21][22],μ is the renormalization scale in a particular scheme, and the spacelike Wilson line is Unlike the definition in Eq. 1, which is invariant under a Lorentz boost along the z direction, the quasi-PDF changes dynamically under such a boost and depends nontrivially on the nucleon momentum P z . For a nucleon of mass M N moving with finite but large momentum P z M N , Λ QCD , LaMET allows us to match the quasi-PDF to the PDF through a factorization formula [17,18]: where C is the matching kernel, and the O(M 2 N /P 2 z , Λ 2 QCD /P 2 z ) terms are power corrections suppressed by the nucleon momentum. Here q(y, µ) for negative y corresponds to the antiquark contribution. Theq and q have the same infrared (IR) divergences, so the matching kernel C depends on ultraviolet (UV) physics only and, thus, can be calculated in perturbative QCD.
There has been rapid development following Ji's proposal. Lattice-QCD calculations of the proton isovector quark distribution f u−d [23][24][25][26], including the unpolarized, polarized and transversity cases, as well as the pion distribution amplitude [27], have been carried out within the LaMET approach. The one-loop matching kernels were calculated in the continuum theory for the isovector quark distributions in a transverse-momentum cutoff scheme in Ref. [19] and reproduced in Refs. [24,28]; the matching for GPDs was addressed in Refs. [29,30]. Recently also studies in lattice perturbation theory are available [31][32][33]. The nucleon-mass corrections to all orders in M 2 N /P 2 z have already been derived in Refs. [23,25] and included in the lattice calculations [23,25], while the higher-twist O(Λ 2 QCD /P 2 z ) correction was numerically removed by fitting the results at different P z with a polynomial of 1/P 2 z and extrapolating to infinite momentum [23,25].
Despite many promising features in the previous lattice calculation of PDFs [23][24][25][26][27], one important piece was missing until recently to form a complete image: the lattice renormalization of the quasi-PDFs. The UV transverse-momentum cutoff scheme used in the one-loop matching computation [19,28] is not the same regularization as used on the discretized lattice. To reduce systematic uncertainties from this mismatch, a proper renormalization of the bare lattice matrix elements is required. An alternative approach is to replace the lattice regularization by the gradient flow and match the continuum extrapolated results to the MS PDF [34], where the latter will be rather complicated due to the new vertices introduced by the gradient flow [35]. With larger statistics and the momentum-smearing technique, which allows high momentum with small statistical errors [36], the uncertainty of lattice simulations will soon be dominated by the renormalization, which, therefore, need to be properly addressed.
The renormalization of the quasi-PDF has been closely studied from the perturbative point of view [19,28,[37][38][39][40][41][42]. The bare quasi-PDF suffers from both logarithmic and linear UV divergences [19,28]. The linear divergence originates from the self-energy of the spacelike Wilson line W z (z, 0) and can be absorbed into an exponential factor exp(δm|z|) where δm has a mass dimension [43][44][45]. This linear divergence is not affected when the Wilson line is inserted between two separated quark fields, so the exponential factor is capable of removing the same divergence in the quasi-PDF [38,39]. Since the remaining divergences are logarithmic and can be subtracted by a renormalization factor that only depends on the endpoints, the quasi-PDF has been proven to be multiplicatively renormalizable forh(z, P z ,μ) at any individual z [41,42].
This property makes it possible to carry out a nonperturbative renormalization of the quasi-PDF in the regularization-invariant momentum-subtraction scheme (RI/MOM) [46] that has been widely used for quark operators on the lattice. In the RI/MOM scheme, the UV divergence in the quasi-PDF can be removed to all orders in perturbation theory by the renormalization constant determined non-perturbatively, leaving the theoretical uncertainty to how precisely one can match the renormalized quasi-PDF onto the MS-renormalized PDF. Since the RI/MOM scheme is regularization independent, the matching kernel can be calculated analytically in the continuum theory with dimensional regularization (d = 4 − 2 ), which is free of linear divergence. The oneloop result has already been obtained, and shows nicely convergent features for Eq. 5, compared to the matching in the transverse-momentum cutoff scheme [47].
In Ref. [40], the additional mixing due to the chiral symmetry breaking on the lattice was first identified in the 1-loop perturbative calculation of O γz , where O γt can be free of this problem. We will come back to this mixing effect in the following sections. Note that O γt is also a feasible choice for the quasi-PDF in LaMET as it belongs to the same universality class of O γz [48]. In a recent new proposal [22,49] O γt is used to define an Ioffe-time or pseudo distribution, which requires the same lattice setting and similar factorization formula to extract out the PDF [50].
In this work, we present lattice results for the nonperturbatively renormalized quasi-PDF in RI/MOM scheme for the isovector unpolarized case 1 , and match it to the MS PDF at one-loop order in perturbative QCD.
We demonstrate the procedure with the previously calculated lattice quasi-PDF [23,25] using clover valence fermions on N f = 2+1+1 (degenerate up/down, strange and charm) flavors of highly improved staggered quarks (HISQ) [52] generated by MILC Collaboration [53] with lattice spacing a = 0.12 fm, box size L ≈ 3 fm and pion mass m π ≈ 310 MeV. The presentation of the paper is organized as follows: In Sec. II, we provide the theoretical setup of the RI/MOM renormalization and explain how to implement it on the lattice. In Sec. III, we show the result of the renormalization factor, and use it to renormalize our previous quasi-PDF [23,25] obtained on the same lattice. We then match the renormalized quasi-PDF in the RI/MOM scheme to the PDF in MS scheme following the procedure elaborated in Ref. [47]. In Sec. IV, we summarize our results and discuss possible directions for further studies.

II. RENORMALIZATION OF WILSON-LINK OPERATORS
For continuum QCD, the renormalization of nonlocal quark bilinear operators has been discussed since the 1980s [43][44][45]54], and the multiplicative renormalizability of the operator has been proven [41,42]. Recent studies based on one-and two-loop perturbative analysis [37][38][39] also indicate that this property might be valid to all orders. Under this assumption, operator mixing does not appear for nonsinglet operators in renormalization. Here, we address the situation in the lattice case, when certain symmetries are broken.

A. Operator Mixing
On the Euclidean lattice, QCD is invariant under discrete symmetries, which include parity P, time reversal T and charge conjugation C. The parity and time-reversal operation are generalized into any direction in the Euclidean space. Because there is no distinction between time and spatial directions, we call the generalized parity and time-reversal operations P µ and T µ , respectively. We investigate the transformation properties of the nonlocal operator O Γ (z) and as some of the discrete transformation can flip the sign of z, it is convenient to define the combinations [55]  (Hermitian). The transformation properties of C, P µ and T µ prohibit O γz (z) from mixing with other operators except for O I (z), where I is the identity matrix. In the zero quark mass limit, we have chiral symmetry (a continuous symmetry), which eliminates the mixing between O γz (z) and O I (z). Some lattice fermions, such as Wilson-type fermions, explicitly break chiral symmetry and introduce a mixing between O γz (z) and O I (z). The situation for the other vector operators, Γ = γ x , γ y and γ t , is different. Discrete symmetries alone prohibit their mixing with other Γ's even if chiral symmetry is broken.
The same discussion can also be applied to pseudoscalar, axial vector, and tensor operators [55]. Our analysis is consistent with what was found in one-loop lattice perturbation theory [40,51].
Therefore, the renormalization of the nonlocal vector operators for lattice fermions without chiral symmetry can be schematically presented as where all Z's are complex functions. For the diagonal elements, Re[Z 11 (22) In the past, Γ = γ z has been chosen for the unpolarized quark distributions. As we discussed above, the renormalization for this operator involves mixing with the scalar operator, whose signal is generally worse in lattice simulations. Alternatively, as pointed out in Ref. [19], one can choose Γ = γ t instead of Γ = γ z to define the unpolarized quasi-PDF. This choice also approaches the normal PDF in the infinite-momentum limit and has the advantage of avoiding the mixing problem. However, the matching kernel, which involves vectors in the z and t directions, becomes much more complicated in this case. Therefore, we leave it for future investigation, and concentrate in this work on Γ = γ z and will estimate the mixing effect from scalar operator nonperturbatively.
where n µ = (0, 0, 0, 1) is the unit vector along the z direction and the summation is over all lattice sites w. The quark propagators are defined as S(p, x) = y e ipy ψ (x)ψ(y) , S(p) = x e −ipx S(p, x) .
(10) and two γ 5 are inserted to the both side of S † (p, w + zn) in Eq. (9) to get the necessary propagator y e −ipy ψ (y)ψ(w + zn) . By imposing the RI/MOM renormalization condition, where the superscript R denotes a renormalized quantity, and µ R is the renormalization scale. Note that the vertex functions are projected withΓ = / p/p z to avoid the ambiguity arising from additional operator mixing in the offshell matrix elements [47], and the prescription of equating the proton momentum P z to the quark momentum p z is used. The renormalization matrix Z(z, p z , a, µ R ) with lattice spacing a is inverse ofZ in Eq. (7,) which can be extracted via We drop the renormalization of the quark self energy, since it only contributes to the overall constant factor, which can eventually be determined by normalizing q(x, µ)dx to unity. It is equivalent to normalize the vector charge to unity. Such a strategy have been used in many recent nucleon matrix elements calculations like Ref. [56][57][58], and the origin of this idea can be traced to the original reference of the RI/MOM scheme, which used the discrete conserved vector current to define the renormalization constant of the quark self energy [46].
Next, the renormalized proton matrix element of where Det(Z) is the determinant of the renormalization matrixZ. The a dependence on the right-hand side cancels up to discretization errors of order O(aP z , aµ R ). The renormalized quasi-PDFq R (x, P z , µ R ) in the RI/MOM scheme can be obtained by a Fourier transform: In the next section, we apply the RI/MOM scheme to the renormalization of quasi-PDF on the lattice, and eventually extract the MS PDF through a sequence of systematic corrections.

III. LATTICE CALCULATIONS
The results of our lattice calculations are presented in two parts: The first part is the nonperturbative renormalization constants in RI/MOM scheme, the second part is the result of the isovector unpolarized PDF. The bare quasi-PDF is renormalized using the renormalization in the first part, and then matched to the PDF using the one-loop matching formula after the power corrections in P z are applied. The final result is the isovector unpolarized PDF of the proton in the MS scheme.

A. Renormalization Constants in the RI/MOM Scheme
For the renormalization calculation, we used 33 configurations on the L 3 × L t = 24 3 × 64 lattice [53] used for the previous quasi-PDF calculation [23,25]. We connected the ends of the quasi-PDF operator to the sinks of the momentum-source quark propagators with p = 2π(3/L, 2/L, 3/L, 8/L t ), which enables us to take the volume average of the operator position as in Eq. (9), which improves the signal-to-noise ratio. This treatment allows us to access all the operators with different Wilsonlink lengths, although we must repeat the calculation for different momenta.
The p z = 6π/L computation used p = 2π(3/L, 2/L, 3/L, 8/L t ) while the p z = 4π/L computation uses p = 2π(3/L, 3/L, 2/L, 8/L t ). To make the RI/MOM scheme µ 2 R = p 2 to be the same, we use the same p for the latter case as the former one but change the operator to ψ(y)γ y W y (y, 0)ψ(0). In both cases, µ 2 R = 5.74 GeV 2 . In the renormalization procedure, there is an arbitrariness to define the renormalization constant as long as it subtracts all the UV divergences. In our case, the renormalization constant should depend on how we choose µ R and p z , and p z = P z is just one of the choices that we can make. This P z dependence shall be cancelled in the matching, as the RI/MOM scheme dependence must be cancelled in our matching up to O(α 2 s ) and the the discretization effects. The discretization effect introduced by this choice can be eliminated eventually when the simulation with the same p z (in the physical unit) are repeated at smaller lattice spacings and the continuum extrapolation are applied on the renormalized results with different lattice spacings.
A comparison of the signals between the point source and the momentum source for the p z = 6π/L case is given in Fig. 1. It is obvious that with the same configurations, the signal with the momentum source can be much better than that with the point source, while the central values are consistent with each other.   2 shows both the renormalization factor and the mixing with the scalar quasi-PDF operator O I (z) for the p z = 6π/L case. We would estimate the systematic uncertainty from ignoring the mixing by assuming O I (z) ∼ iO γz (z) with the overall factor i from the wick rotation, then h R (z) will not change at z = 0, but it will change about ∼ 8% at z = 6 and ∼ 15% at z = 12 respectively, which is smaller than the present statistical uncertainties.
Since the matching from quasi-PDF to PDF will be processed in the space of the momentum fractions, we will check the µ R dependence on the final distributions instead of that on the effective MS renormalization constants.

B. From Quasi-PDF to PDF: Numerical Results and Discussion
In this subsection, we present our results for the unpolarized isovector quark distribution. We first calculate the time-independent, nonlocal correlator of a nucleon with finite P z h Γ (z, µ, P z ) = P ψ (z)Γ n U z (nẑ) ψ(0) P , (16) where U z is the gauge link pointing from nẑ to (n + 1)ẑ, and P = (0, 0, P z ) is the momentum of the nucleon. We calculate the bare lattice nucleon matrix elements h γz and h I at P z = {1, 2, 3}2π/L, which are 0.43, 0.86 and 1.29 GeV, respectively. As observed in Refs. [25,27], the correction terms for the smallest-momentum distribution is less well-behaved; thus, we drop it in the rest of this work. We then renormalize the bare matrix elements with the RI/MOM renormalization factors defined in the previous section: where the mixing with h I turns out to be numerically negligible because h I /h γz M/P z and |Z SV /Z V V | 1. In Fig. 3, we show the bare (h γz ) and renormalized (h R ) matrix elements for P z = {2, 3}2π/L. In the renormalized matrix elements the mixing effect is temporarily ignored. We note that in both cases, the bare matrix elements vanish within error bands when the link length reaches 10a-12a. After renormalization, the error bands become much broader at large z due to an exponential increase of the renormalization factor, and consistent with 0 within error bands.
Next, we Fourier transform Eq. (15) to convert the lattice matrix elements as functions of spatial link length z into the quasi-PDF with µ R the RI/MOM renormalization scale. Then we take the one-loop RI/MOM-to-MS  matching calculated in Ref. [47] and mass corrections for the renormalized quasi-PDF. We invert Eq. (5) to obtain the PDF in the MS scheme, where C (1) has been computed in Ref. [47], where and with r R = µ 2 R /(p z R ) 2 , and the Feynman and Landau gauges correspond to τ =1 and 0 respectively. p z R is the momentum used in the RI/MOM renormalization condition, whereas p z is the parton momentum for the matching. In our calculation, we choose p z R = P z , and p z has to be |y|P z in the factorization theorem of Eq. (18). "⊕, +, and " denote plus functions whose defintions are given in [47].
Theq M (x, P z , µ R ) in the above equation is the quasi-PDF in the RI/MOM scheme with the nucleon mass correction removed [23,25], where c = M 2 N /P 2 z , f + = √ 1 + c + 1 and c ≡ c/f 2 + < 1 for any P z . The remaining Λ 2 QCD /P 2 z correction will be removed by a parametrization, as was done in Ref. [25]. The µ R dependence on the right-hand side should cancel modulo residual O(a 2 µ 2 R , α 2 s ) corrections. The case with P z = 6π/L and three difference µ R is illustrated in Fig. 4. As in the figure, the residual O(a 2 µ 2 R , α 2 s ) dependence are smaller than the statistical uncertainties.
The final results are shown in Fig. 5. In contrast to the previous result in Ref. [25], the sea flavor asymmetry is hardly visible, mainly due to the rapid increase of the renormalization factor with distance, which amplifies the error. The peak in the positive-x region is shifted slightly to the left. This is expected since the renormalization enhances the long-range correlation, and thereby enhancing the contribution in the small x region when Fourier transformed to momentum space. After renormalization the unphysical dip near x = 0 in the previous result also vanishes. This is because the linear divergence is removed and, therefore, the RI/MOM matching kernel has a smoother form than the matching used to relate bare PDFs.
Another observation concerning our renormalized distribution is an oscillating behavior in negative-x (antiquark) region, which is absent from the previous bare-PDF results. This is likely because the bare matrix elementh(z) decays very fast with the distance z, so the long-range correlation plays a less important role. However, the long-range correlation becomes more important in the renormalized distributions due to the exponential increase in the renormalization factor at large distance. In such a case, the cut-off on z will introduce large truncation errors by Fourier transforming theh(z) intoq(x) [59]. To examine this hypothesis, we apply the reverse matching and mass corrections procedure to the central values of one of the global fitted of PDF, "CJ15" (from the CTEQ-JLab collaboration [60]), to make a direct comparison with our renormalized function h(z). Note that the PDF community fits xq(x) and has larger uncertainty near x = 0, so their h(z) will also have large uncertainties at very large z. The result is shown in Fig. 6. The renormalizedh(z) at 310-MeV pion mass and a relatively small lattice (L=2.9 fm) has a narrower peak around zp z = 0 and differs significantly from the Fourier transform of the CJ15 result at large values of zp z . Further studies on removing the higher-twist contribution at large z in the RI/MOM renormalization will be carried out in the future.
Finally, we have several comments regarding our RI/MOM treatment. The first one is the possible gauge dependence induced by taking the external quarks offshell in the nonperturbative renormalization. The gauge dependence should be canceled by the matching kernel, but the cancellation is not complete, since the kernel is only computed at one loop. It is encouraging that the one-loop matching effect is numerically small in Landau gauge that we employed. Whether the higher-loop contributions will remain small requires further study. The second one is treating p z of the off-shell quark the same as the proton P z . Numerically, the renormalization factor is rather insensitive to p z and µ R , so we do not expect this treatment to cause a big error.

IV. SUMMARY
We have carried out a nonperturbative renormalization of the quasi-PDF in the RI/MOM scheme in lattice QCD. Based on the renormalized quasi-PDF, we have updated the lattice result of the unpolarized isovector quark distribution from previous studies by some of the authors. The RI/MOM renormalization of the quasi-PDF is performed for each individual z, where it has been shown to be multiplicatively renormalizable. All the UV divergences, including the linear and logarithmic divergences, are subtracted nonperturbatively by the renormalization constant. Meanwhile, due to chiral symmetry breaking from the lattice fermion action we used, there is a mixing between the isovector quasi-PDF and a scalar operator. The mixing is estimated to be a ∼ 10% effect and is left for future investigation.
It is possible that an alternative process is needed in the large z region. But the errors there are still too large to conclude. Nevertheless, there is an alternative renormalization method used in Ref. [41,61] which is worth investigating. Also, two methods to reduce the weighting of the long range correlation are also discussed in Ref. [59]. Clearly more efforts on this issue in the future are needed.
Compared to the previous results on bare PDFs, our present result is free of the unphysical dip at x = 0 due to the smooth matching kernel. However, we end up with a large uncertainty band that makes it difficult to evaluate whether an improvement has been achieved. The reason behind the large uncertainty band is that the RI/MOM renormalization constant which grows exponentially at large z significantly amplifies the error in the nucleon matrix element of the quasi-PDF. Future work involving higher momentum, larger volume, lighter pion mass, and finer lattice spacing (such that the higher nucleon boosted momenta P z can be used without the additional (P z a) n systematics) or other renormalization conditions may resolve some of the issues we see in this paper.