A complete non-perturbative renormalization prescription for quasi-PDFs

In this work we present, for the first time, the non-perturbative renormalization for the unpolarized, helicity and transversity quasi-PDFs, in an RI' scheme. The proposed prescription addresses simultaneously all aspects of renormalization: logarithmic divergences, finite renormalization as well as the linear divergence which is present in the matrix elements of fermion operators with Wilson lines. Furthermore, for the case of the unpolarized quasi-PDFs, we describe how to eliminate the unwanted mixing with the twist-3 scalar operator. We utilize perturbation theory for the one-loop conversion factor that brings the renormalization functions to the MS-scheme at a scale of 2 GeV. We also explain how to improve the estimates on the renormalization functions by eliminating lattice artifacts. The latter can be computed in one-loop perturbation theory and to all orders in the lattice spacing. We apply the methodology for the renormalization to an ensemble of twisted mass fermions with Nf=2+1+1 dynamical light quarks, and a pion mass of around 375 MeV.


Introduction
Parton distribution functions (PDFs) describe the inner dynamics of partons inside a hadron [1]. They have a non-pertubative nature and, thus, they can not be computed in perturbation theory. Lattice QCD is an ideal formulation to study the PDFs from first principles, in large scale simulations. However, PDFs are usually defined on the light cone, which poses a problem for the standard Euclidean formulation. Hence, hadron structure calculations in lattice QCD are related to other quantities that are accessible in a Euclidean spacetime. This led to a long history of investigations of Mellin moments of PDFs and nucleon form factors (see [2][3][4][5][6] for recent reviews). In practice, there are severe limitations in the reconstruction of the PDFs mainly due to the small signal-to-noise ratio for the high moments. In addition, there are inevitable problems with power divergent mixings with lower dimensional operators. Therefore, the task of reconstructing the PDFs from their moments is practically unfeasible.
Ab initio evaluations of PDFs are of high importance as they would be a stringent test of non-perturbative aspects of QCD. The fact that a calculation of the PDFs from first principles is missing is a pressing problem that prevents a deeper understanding of the nucleon structure. Our current knowledge on the PDFs relies on phenomenological fits from experimental data using perturbation theory. These parameterized PDFs serve, for example, as input for the computation of cross sections used for presently running colliders, most notably the LHC, and also to plan future collider experiments which are themselves tests of the Standard Model. The parameterizations are, however, not without ambiguities [7]. In addition, there are kinematical regions that are not experimentally accessible. The large Bjorken x region is one of them, with large uncertainties most dramatically seen in the down quark distributions [8,9]. The transversity PDF is yet another example of a PDF that is only poorly constrained by phenomenology.
A pioneering method for a direct computation has been suggested by X. Ji [10]. In this approach, instead of matrix elements defined on the light cone, one calculates matrix elements of fermion operators including a finite-length Wilson line, and whose Fourier transform defines the so-called quasi-PDFs. This is achieved by taking the Wilson line in a purely spatial direction, conventionally chosen to be the z-direction, instead of the +-direction on the light cone. In terms of hadron kinematics, the nucleon momentum is usually taken along this spatial direction. At large, but finite momenta, the quasi-PDFs can then be related to the light-front PDFs via a matching procedure [11,12]. Ji's approach has already been tested in Refs. [13][14][15][16][17] and all these results are promising, i.e. they give a correct shape of the PDFs after the matching procedure. Certain properties of quasi-PDFs, like the nucleon mass dependence and target mass effects, have also been analyzed via their relation with transverse momentum dependent distribution functions (TMDs) [18,19]. Note also the appearance of a related approach, pseudo-PDFs, which is a differ-ent generalization of light-cone PDFs to finite nucleon momenta [20,21]. Refs. [22,23] discussed the role of the Euclidean signature in the computation of quasi-PDFs, as compared to light-front PDFs, and the latter reference proved that matrix elements obtained from Euclidean lattice QCD are identical to those obtained using the LSZ reduction formula in Minkowski space. Nevertheless, the matrix elements of quasi-PDFs contain divergences that need to be eliminated via renormalization in order to obtain meaningful results that can be compared to the physical PDFs.
To date, all works on the quasi-PDFs only considered the bare matrix elements, as the renormalization process is highly non-trivial and was not addressed until recently. In particular, new complications arise in the renormalization of the Wilson line operators, compared to the local operators. For one, in addition to the logarithmic divergences there is a linear divergence [24] with respect to the lattice regulator, a, that prohibits one to take the continuum limit prior to its elimination. To one-loop level in perturbation theory the divergence is manifesting itself as a linear divergence, computed in Ref. [25] for a variety of fermion and gluon actions. However, it is of utmost importance to extract the power divergence non-perturbatively, which is one of the goals of this work. Another feature of these operators that brings in new complications is the fact that certain choices of the Dirac structure exhibit mixing [25].
To show that these matrix elements can be renormalized, in particular that the linear divergence associated with the Wilson line can be eliminated, is of paramount importance. Without renormalization, the whole quasi-PDF strategy is incomplete and unable to provide any useful information to the theoretical and experimental community. Some suggestions for the elimination of the linear divergence via the static potential were proposed in Refs. [26,27]. In Ref. [12] a one loop calculation of the linear divergence has been made, and that motivated the definition of an improved quasi-PDF that is free of power divergences. One has, nevertheless, to show that such a procedure can be done non-perturbatively. Another method to extract the coefficient of the linear divergence using the nucleon matrix elements of the quasi-PDFs was also presented in Ref. [25]. An alternative technique to suppress the linear divergence was discussed in Ref. [28] utilizing the gradient flow. Very recently, two papers discussed the employment of the auxiliary field formalism for the renormalization of quasi-PDFs [29,30], where the Wilson line is replaced by a Green's function of the introduced auxiliary field. The renormalizability of quasi-PDFs to all orders in perturbation theory was addressed in Ref. [31].
In this paper we propose, for the first time, 1 a concrete renormalization method of the quasi-PDFs in a fully non-perturbative manner. We provide the prescription of the method and show examples of the renormalized matrix elements. We also discuss the elimination of the mixing between the unpolarized quasi-PDF and the twist-3 scalar operator. We employ the RI renormalization scheme [33] and we convert the results to the MS scheme at a reference scale, μ = 2 GeV, using the one-loop conversion factor computed in Ref. [25]. As a test case, we focus on the helicity quasi-PDF to demonstrate results, as it is free of mixing.
The outline of the paper is as follows: In Section 2 we provide the theoretical setup related to the nucleon matrix elements and the renormalization prescription in the presence and absence of mixing, as well as the data on the conversion factor computed perturbatively. Section 3 includes results on the renormalization functions, a discussion on the systematic uncertainties in the Z-factors, renormalized matrix elements for the helicity case, as well as results of the matching to light-front PDFs. Finally we conclude and give our future directions.

Theoretical setup
In this section we briefly introduce the nucleon matrix elements for the quasi-PDFs that we aim to renormalize. We also explain the renormalization prescription for the three types of PDFs: unpolarized, helicity and transversity. For details on the computation of the nucleon matrix elements, we refer to Refs. [15,17].

Nucleon matrix elements
We consider matrix elements of non-local fermion operators that contain a straight Wilson line, denoted by h (P 3 , z). The variable z is the length of the Wilson line, and P 3 is the nucleon momentum, which is taken in the same direction as the Wilson line. The quasi-PDFs can be computed from the Fourier transform of the following local matrix elements: where |N is a nucleon state with spatial momentum P = (0, 0, P 3 ) along the 3-direction and W 3 (z) is a Wilson line of length z in the same direction. denotes the Dirac structure of the operator insertion, which is γ μ (unpolarized), γ μ · γ 5 (helicity), σ μν (transversity). In the works appearing in the literature γ μ is taken along the Wilson line. In principle, one may choose γ μ orthogonal to the direction of the Wilson line. In this case, the unpolarized operator is free of mixing, while the helicity and transverity do mix. For example, choosing the γ -matrix in the temporal direction is important for a faster convergence to the physical PDFs, as discussed in Ref. [19]. Our recent work [25] indicates that P 3 ⊥ z ( P 3 z) is ideal for the unpolarized (helicity and transversity) case.
To calculate the bare matrix elements, we use the setup of Ref. [17]. We consider one ensemble of dynamical N f = 2+1+1 twisted mass fermions produced by ETMC [34], with volume 32 3 ×64, lattice spacing a≈0.082 fm [35] and a bare twisted mass of aμ = 0.0055, which corresponds to a pion mass of around 375 MeV. We performed our calculations on 1000 gauge configurations with 15 forward propagators and 2 stochastic propagators, i.e. 30000 measurements in total. We will present results for momentum P 3 = 6π L , which is around 1.4 GeV in physical units. Gaussian smearing has been employed on the nucleon interpolating fields in the calculation of the matrix elements [17].

Renormalization scheme
Here we discuss a fully non-perturbative renormalization prescription that will remove all divergences inherited in the matrix elements of the quasi-PDFs, as well as the mixing, as indicated in the perturbative analysis of Ref. [25]. In a nutshell, the proposed renormalization program: 1 removes the linear divergence that resums into a multiplicative exponential factor, e −δm|z|/a+c|z| . The coefficient δm represents the strength of the divergence and is expected to be operator independent, as it is related only to the Wilson line. c is an arbitrary scale [36] that can be fixed by such a renormalization prescription; 2 takes away the logarithmic divergence with respect to the regulator, log(aμ 0 ), where μ 0 is the RI renormalization scale; 3 applies the necessary finite renormalization related to the lattice regularization; 4 eliminates the mixing that appears in the unpolarized operator, as the bare matrix element is a linear combination of the unpolarized quasi-PDF and the twist-3 scalar operator. The two may be disentangled by the construction of a 2×2 mixing matrix.
We adopt a renormalization scheme which is applicable non-perturbatively, that is, the RI scheme [33]. We compute vertex functions of the operators under study, between external quark states, with the setup being in momentum space, and the operator defined as: where = γ μ , γ μ · γ 5 , σ μν (ν = μ). The path ordering of the exponential appearing in the above expression becomes, on the lattice, a series of path ordered gauge links. The renormalization functions (Z-factors) depend on the length of the Wilson line and, thus, we perform a separate calculation for each value of z. Typically, z goes up to half of the spatial extent of the lattice. The renormalization prescription is along the lines of the program developed for local operators and the construction of the vertex functions is described in Ref. [37]. The difference between the renormalization of the local operators and the Wilson-line operators is the linear divergence that appears in the latter case. However, there is no need to separate this divergence from the multiplicative renormalization and, therefore, the technique described below may successfully extract both contributions at once.

Helicity and transversity quasi-PDFs
We first provide the methodology for a general operator with a Wilson line in the absence of any mixing. This is applicable for the helicity and transversity quasi-PDFs, provided that their Dirac structure is chosen along the Wilson line. The renormalization functions of the Wilson-line operators, Z O , are extracted by imposing the following conditions: where Z q is the renormalization function of the quark field obtained via The trace is taken over spin and color indices, and the momentum p entering the vertex function is set to the RI renormalization scale μ 0 . In Eq. (3) V(p, z) is the amputated vertex function of the operator and V Born is its tree-level value, i.e. V Born (p, z) = iγ 3 γ 5 e ipz for the helicity operator. Also, S(p) is the fermion propagator and S Born (p) is its tree-level. The RI scale μ 0 is chosen such that its z-component is the same as the momentum of the nucleon. Such a choice serves as a suppression of discretization effects, as different classes of spatial components have different discretization effects, and scales of the form (n t , 3, 3, 3) have small discretization effects [38]. We test both diagonal (democratic) and parallel momenta to the Wilson line (in the spatial direction), that is a μ 0 = 2π L (P 3 , P 3 , P 3 ) and a μ 0 = 2π L (0, 0, P 3 ), respectively. We will refer to these choices as "diagonal" and "parallel". The latter are expected to have larger lattice artifacts, as the ratiô is higher than for diagonal momenta. Using renormalization scales leading to a small value for such a ratio has been successful for the local fermion operators [39,37]. The vertex functions V(p) contain the same linear divergence as the nucleon matrix elements. This is crucial as it allows the extraction of the exponential together with the multiplicative Z-factor, through the renormalization condition of Eq. (3). Regarding the renormalizability of quasi-PDFs, it was proven to be multiplicative to all orders in Refs. [27,31]. The authors consider quasi-PDFs defined in coordinate space and prove the multiplicative renormalizability, and the same holds for the definition in Eq. (1). This is due to the fact that the former require, in addition, a consistent subtraction scheme to remove all terms divergent in the limit ξ z → 0. Moreover, the renormalization of bilocal composite operators is also studied in Ref. [27,31]. Based on the above, Z O can be factorized as where Z O is the multiplicative Z-factor of the operator and δm is the strength of the linear divergence. The exponential includes a term with an arbitrary scale c that could be of the form c|z| [36]. Note that the exponential comes with a different sign compared to the nucleon matrix element, as Z O is related to the inverse of the vertex function, Such a construction of the Z-factor justifies the reason why the elimination of the power divergence in the nucleon matrix elements is successful by multiplying with Z O , provided that it has been calculated on the same ensemble. Consequently, the above renormalization condition handles all the divergences which are present in the matrix element under consideration. Note that in the absence of a Wilson line (z = 0), the Z-factors reduce to the ones for local currents, which are free of any power divergence. For example, for the helicity operator, Z A (z = 0) is the standard Z-factor of the axial current. We would like to point out that knowledge of the coefficient δm provides insight on the strength of the power divergence. One can pursue this direction via the static potential [27] or the technique proposed in Ref. [25]. If the linear divergence is extracted, one may apply the matching of Ref. [12], which includes the coefficient δm. However, there is still a necessity to compute the multiplicative Z-factor to cure any logarithmic divergences, apply finite renormalization, as well as fix the arbitrary scale c.

Unpolarized quasi-PDF
The case of the unpolarized quasi-PDF requires special treatment, if the Dirac structure is in the same direction as the Wilson line. As demonstrated in Ref. [25], for such a choice there is a mixing with the twist-3 scalar operator, 2 that must be eliminated. This mixing appears in some lattice regularizations (in particular, Wilson and twisted mass fermions) due to the breaking of chiral symmetry, and is found to be finite. We establish notation by using the subscripts S and V for the scalar and vector (unpolarized) operators, respectively. The corresponding operators are: and we represent their nucleon matrix elements as h S (P 3 , z) and h V (P 3 , z). To disentangle the two contributions from their bare matrix elements, one must compute the multiplicative renormalization and mixing coefficients from the following 2×2 matrix: where Ẑ is the matrix of the mixing between the scalar and vector operators, that iŝ According to the above mixing, the renormalized unpolarized quasi-PDF, h R V (P 3 , z), is related to the bare scalar and unpolarized via: where Z V V and Z V S are computed in a certain scheme. In the present work we employ the RI scheme and then convert to the MS scheme, at an energy scale μ. The Z ii factors can be computed following a prescription similar to Eq. (3), that is: where the elements of the vertex function matrix V are given by the trace In the above equation V Born i is the tree-level expression of the operator O i . Thus, all matrix elements of Ẑ can be extracted by a set of linear equations, which can be written in the following matrix form: As can be seen from Eq. (12), a major component in the renormalization of the unpolarized quasi-PDF is knowledge of the scalar nucleon matrix element h S (P 3 , z). To date, no lattice calculation is available for h S (P 3 , z), as a mixing was not anticipated prior the work of Ref. [25]. Therefore, a proper renormalization of the unpolarized quasi-PDF is still pending. However, the mixing appears to be greatly suppressed in the presence of a clover term in the fermionic action. This is of high importance, as we are currently computing the quasi-PDFs on an ensemble of twisted mass clover improved fermions at the physical pion mass [40,41].

Conversion to the MS-scheme
The fermion operators with a finite-length Wilson line are scheme and scale dependent. As a result, having obtained the Z-factors in the RI scheme as depicted in Eqs. (3) and (13) at the RI scale μ 0 , we must convert them to the MS scheme at a scale μ. This conversion factor has been computed in one-loop continuum perturbation theory in Ref. [25]. In addition, comparison with phenomenological estimates is done typically at μ = 2 GeV. This defines the value of the MS renormalization scale chosen in the expressions for the conversion factors. Alternatively, one can choose μ =μ 0 , and then evolve the scale to 2 GeV, which requires the anomalous dimension of the operator: Note that to one-loop level, the anomalous dimension does not dependent on the Dirac structure of the operator and is the same in the RI and MS schemes. The evolution to the MS renormalization scale μ = 2 GeV can be performed using the intermediate Renormalization Group Invariant scheme (RGI), as employed in Refs. [42,38]. In the framework of this paper we tested both methods and their difference was found to be negligible. As discussed in the previous section, we employ two types of the RI renormalization scale regarding the spatial direction: one the same as in the nucleon matrix elements, a μ 0 = 2π L (0, 0, P 3 ), (parallel to the Wilson line direction), and one which is diagonal and each direction has a value of P 3 , that is a μ 0 = 2π L (P 3 , P 3 , P 3 ). The conversion factor depends on both the RI and MS renormalization scales. While μ is typically fixed to 2 GeV, μ 0 can change affecting the numerical values of the conversion factor. Such a case is illustrated in Fig. 1 for the unpolarized and helicity operators, which share the same conversion factor. 3 We focus on the renormalization of the matrix elements with the nucleon boosted by aP 3 = 6π L , and we apply the same to the conversion factor, for the two cases of the renormalization scale. For the specific value of P 3 , we choose the temporal direction of μ 0 such that the ratio P defined above, is as small as possible in order to suppress lattice artifacts. Nevertheless, for the "parallel" case the minimum value for the ratio is P = 0.54 which is already very high. The chosen values for aμ 0 are: 2π 32 ( 7 2 + 1 4 , 3, 3, 3) and 2π 32 ( 4 2 + 1 4 , 0, 0, 3) for the "diagonal" and "parallel" case, respectively.
The real and imaginary parts of the conversion factor are plotted as a function of the length of the Wilson line, z. We stress that the dependence of the conversion factor on the renormalization scales and the length of the Wilson line is highly non-trivial and is expressed in terms of integrals of modified Bessel functions. Consequently, the data points shown in Fig. 1 have been computed numerically for the specific scales, at each value of z separately. We observe that the real part of the conversion is an order of magnitude larger than the imaginary part. The real part consistently increases for increasing values of z, while the imaginary part almost immediately stabilizes when z becomes non-zero. Also, the conversion factor at z = 0 is equal to unity, as it corresponds to the local vector and axial currents. Note that the case z = 0 is not extracted from the calculation of Ref. [25] as it is strictly for z =0: the appearance of contract terms beyond tree level renders the limit z → 0 nonanalytic. On the contrary, the non-perturbative prescription of the previous section can be applied for any value of z, as the calculation is performed on each z separately. The values of Fig. 1 will be used in the following section to bring the RI non-perturbative Z-factors to the MS scheme. To reliably extract the Z-factors we have extend the calculation including several values of μ 0 as explained in Section 3. The conversion factor was found to have the same qualitative behavior for all values of μ 0 .

Results
In this section we apply the prescription suggested above and we present our results for the non-perturbative Z-factors both in the RI and the MS schemes. For demonstration purposes we focus on the data with 5 steps of Hypercubic (HYP) smearing that suppress the power divergence and bring the results closer to renormalized nucleon matrix elements. We focus on the multiplicative renormalization for the helicity quasi-PDF, and only briefly discuss the case of the unpolarized operator.

RI scheme and conversion to the MS scheme
As a starting point, we have applied the two values of the RI scale μ 0 used in Fig. 1 ("parallel" and "diagonal"). After converting both cases to the MS scheme at μ = 2 GeV, we can quantify the systematic uncertainties related to lattice artifacts and truncation of the conversion factor, as explained in the next subsection.
In Fig. 2, we show the extracted values for the helicity Z-factor, Z h , that renormalizes the bare matrix element h(P 3 , z). In each plot we overlay the results for the RI (open symbols) and the MS (filled symbols) schemes, for the real and imaginary part of the Z-factor. We employ the momentum source technique [49,38] that offers high statistical accuracy with a small number of measurements, typically of O (10). As can be seen in the plots, the statistical uncertainties are almost invisible. The left (right) plot corresponds to the "parallel" ("diagonal") choices for μ 0 . We find that the imaginary part of Z MS h is reduced compared to its counterpart in Z RI h , and is also rather small for low values of z. This is more pronounced in the right plot of Fig. 2 for "diagonal" μ 0 , for which Im [Z MS h ] (blue circles) is smaller than the corresponding data from the "parallel" scale, especially for large values of z. It is worth mentioning that the perturbative Z-factor in dimensional regularization and in the MS scheme is real to all orders in perturbation theory, as it is extracted only from the poles. Therefore, it is expected that the imaginary part of the non-perturbative estimates should be highly suppressed. The behavior of the "diagonal" scale is encouraging, as the imaginary part is very close to zero for |z| up to ∼10 a.  Table 1 Values for the RI scale defined as aμ 0 = 2π L ( n t 2 + π 4 , n x , n y , n z ), where L is the spatial extent of the lattice. The values are given in lattice and physical units. The last column corresponds to the ratio P defined in Eq. (5).

Label
(n t , n x , n y , n z ) ( ā μ 0 ) 2μ (GeV)P As we have done in previous applications of the RI scheme to extract the renormalization functions of ultra-local operators (see, e.g., Ref. [38]), here we also use a range of values for the RI renormalization scale, μ 0 . This allows us to study the scale dependence of the vertex functions with the external momentum and identify the window in which renormalization factors can be extracted as reliably as possible. Unlike the case of the local currents, the computation of the lattice artifacts to O(g 2 a ∞ ) is extremely laborious and not available yet. Thus, one has to be careful with the choices for μ 0 , as we want to avoid non-perturbative contaminations, and satisfy the criterion that P (Eq. (5)) is small, preferably below 0.3, to avoid enhanced cut-off effects. We compute the renormalization functions using the values reported in Table 1, from which we choose an optimal range for the fit using the diagonal choices. The parallel μ 0 have only been used to explore the systematic uncertainties in Subsection 3.2.
Since we present the renormalized matrix elements for the helicity PDF, we focus on Z h for this analysis. In Fig. 3 (aμ 0 ) 2 , using the diagonal values labeled by m1-m8. The real (imaginary) part is shown in the left (right) panel. We find a residual dependence on (aμ 0 ) 2 , mostly affecting the imaginary part. This is due to the fact that the vertex function depends not only on the magnitude and direction of the renormalization scale, especially the z-direction, as this is parallel to the Wilson line. Upon evolving to the same scale, the Z-factors should not depend on the initial scale if the evolution is known to higher loops in perturbation theory and discretization effects are sufficiently small. The nonzero slope of the plots shown in Fig. 3 is an indication of non-negligible lattice artifacts, as well as, a consequence of the truncation of the conversion factor to one-loop level. However, further investigation is required in order to attempt disentangling the two effects. We address this in Subsection 3.2.
To remove the residual dependence on (aμ 0 ) 2 we fit Z MS h with the function and we take Z MS 0, h as our final result. This process is done for each value of the length of the Wilson line z, and repeated for several fit ranges for the diagonal choices m1−m8 given in Table 1. For the extrapolated value Z MS 0, h at 2 GeV, we choose the one obtained in the range (aμ 0 ) 2 [1.4, 2.0]. By excluding (aμ 0 ) 2 close to 1, we avoid contamination from nonperturbative effects. The choice for the above fit range also excludes the two higher scales that have P >0.3, which may carry sizable lattice artifacts. As an additional check of the quality of the fit, we find that the χ 2 /d.o.f is small for the chosen range. In Fig. 4, we plot the extrapolated value Z MS 0, h (for 5 steps of HYP smearing), which is applied on the bare nucleon matrix element of the helicity quasi-PDF in Subsection 3.2.
The prescription for obtaining Z V V and Z V S has also been applied on the same ensemble. We find that both the HYP smearing and the choice of "diagonal" scales suppresses the mixing coefficients Z V S and Z SV . In the left panel of Fig. 5, we plot Z RI V V and Z RI V S for the scale "(7,3,3,3)" and 5 steps of HYP smearing. The real and imaginary parts of Z RI V S are of the same magnitude, but at least an order of magnitude smaller than the multiplicative factor Z RI V V . This can be seen from the right panel of Fig. 5   both Z RI V V and Z RI V S is compatible with zero for |z/a| < 4. Therefore, the large values of the blue points in the region |z/a| < 4 in the right plot are no indication of significant mixing.

Assessment of systematic uncertainties
In the framework of this study we have performed several investigations on the systematic uncertainties related to the truncation of the one-loop conversion factor, as well as, discretization effects. In this subsection we present the main conclusions of this study, and we give estimates of these effects.
Prior converting to the MS scheme, we compute the Z-factors for different values of the RI scale μ 0 . For example, Z RI h shown in the two plots of Fig. 2 on the scales m4 and m9 is different. Upon evolving to the same scale, the extracted values should agree if the evolution has converged; typically this happens at two or three loops in perturbation theory. Thus, one should be able to compare the data extracted for Z MS h at 2 GeV for the two scales presented in Fig. 2. The difference between the two, Z(m4, m9), is an indication of the presence of lattice artifacts (mainly in the "parallel" case), coupled with the truncation of the conversion factor. This is demonstrated in Fig. 6 for both the real and imaginary parts, and it is interesting to see that the difference increases as z becomes larger. Based on this observation, we expect that such an increase of the lattice artifacts is also present for the "diagonal" case, but less severe. Another evidence of the presence of non-negligible systematics is the fact that the imaginary part of Z MS h is nonzero. In Fig. 2, we find a small imaginary part for the "diagonal" scale, which has an increasing trend at large values of z.
Based on the above, we conclude that the systematic uncertainties related to lattice artifacts and the conversion factor must be addressed, in order to extract reliable estimates on the renormalization functions. It is our intention to reduce both effects in the near future, which will eliminate systematic uncertainties propagated to the estimates of the quasi-PDFs. Understanding the uncertainties dominating the large-z region is crucial, as the matrix element for these values also enters in the Fourier transform that yields the quasi-PDF. Preliminary explorations indicate that a likely magnitude of the two-loop contribution might suppress the imaginary part of Z MS h . Even though we are currently in no position to accurately quantify these systematics, we will estimate upper bounds. In particular, we try to estimate the effect of each one individually.

Lattice artifacts
The Z-factor for the helicity quasi-PDF at z = 0 reduces to the renormalization function of the local axial current Z A . Z A is scheme-and scale-independent, and thus, the case z = 0 allows us to study the lattice artifacts. Z A has been already evaluated for this ensemble in Ref. [38]. In the latter calculation several diagonal scales have been employed and, together with a technique for the removal of the lattice artifacts, the extracted value was found to be Z A = 0.7556 (5). In the present calculation at z = 0 we find a value of Z A = 0.8620 (15) using the scale m9 and Z A = 0.7727(2) for the scale m4. This is yet another indication that lattice artifacts are large for the "parallel" renormalization scale and less severe, but not negligible, for the "diagonal" case.
An additional test that may be performed for any value of z is the comparison of Z MS h between scales with different components, but same (aμ 0 ) 2 . For this, we utilize the values labeled as m10 and m11, which can be compared to m1 and m4, respectively. Since these pairs are approximately at the same value of (aμ 0 ) 2 , we can compare them directly in the RI scheme without any evolution, 4 and the difference can be interpreted as lattice artifacts. We observe similar behavior for the two pairs and as a demonstration we plot in Fig. 7 the following ratio for each value of z In the left panel we show the case where μ 0 and μ 0 correspond to the values m4 and m11, respectively. The difference from zero can be attributed to lattice artifacts. We find that the Z-factor extracted from m11 deviates from the value extracted from m4 by 5-10% in the real part and up to 30% in the imaginary part. This deviation seems to stabilize after z/a∼8 to 5% (10%) for the real (imaginary) part.
Note that the left plot of Fig. 7 gives the estimate of the discretization effects for the "parallel" case compared to the "diagonal" one. It would be interesting to compute D Re(Im) between two neighboring "diagonal" momenta in order to understand the change in the artifacts. We choose 2 pairs (m2, m3) and (m5, m6) and we find that both D Re (m2, m3) and D Re (m5, m6) are less than 1%, while D Im (m2, m3) and D Im (m5, m6) are of the order of 6%. This is a confirmation that excluding "parallel" momenta in the fit of Eq. (17) reduces significantly the discretization effects.

Truncation effects
In order to assess quantitatively the influence of the conversion factor that is known to oneloop perturbation theory, we form the ratios 4 We confirmed numerically that the conversion factors from m1 to m10, and from m4 to m11 deviates from unity by less than 1‰, so we ignore it. both for the real and imaginary parts of the helicity Z-factor. In Fig. 8, we plot Eq. (19) for the RI' and MS case. These have been extracted at different values of (aμ 0 ) 2 (m1−m8) and evolved perturbatively to the same scale of 2 GeV, using the results of Ref. [25]. The ratio is always taken with respect to the smallest "diagonal" scale, m1. R RI Re(Im) depend on the truncation effects in the one-loop evolution to 2 GeV, while R MS Re(Im) is affected by scheme conversion truncation effects. Without contamination from lattice artifacts and truncation effects, Eq. (19) should equal 1 in both schemes and for all values of (aμ 0 ) 2 and z/a. This realization allows for investigation of the truncation effects.
One observes that the imaginary part is more sensitive to the change of the initial RI scale, which is given in the x-axis. To demonstrate this, we keep the range of y-axis the same for all plots of Fig. 8. The difference between the ratios of the real part extracted from different RI scales (μ 0 ) 2 is consistently small and reaches at most approx. 5% in the RI scheme and 3% in the MS scheme, regardless of the Wilson line length, z. As we mentioned above, for z = 0 we can compare to the highly reliable value found in Ref. [38]; our current value from the scale m1 is around 2% higher. Hence, we can estimate the typical size of discretization effects to be of order 2-5% in the real part of the Z-factors.
The differences in the MS ratios are combinations of lattice artifacts and the conversion truncation. We observe that the conversion to the MS scheme decreases the differences in the MS scheme to at most 3%. Hence, truncation alone in the real part seems not to have an effect larger than 2%. Since we know that conversion truncation effects are very small for small values of z (in particular, they are zero for z = 0, where Z A is scheme-independent), we can take 0-2% as our estimate. In the end, given the fact that truncation effects seem to be opposite to the influence of lattice artifacts, we estimate that the total effects present in the real part of Z MS h should not exceed 5%.
The situation is somewhat different in the ratio of the imaginary parts shown in the right panel of Fig. 8. We conclude that the uncertainty from lattice artifacts can be of the order of 10% in the RI scheme for intermediate and large z (|z|/a ≥ 5 5 ). From the right panel of Fig. 8, the total uncertainty in R MS Im can get enhanced to around 40% for large and up to 80% for intermediate Wilson line lengths. In addition, the total uncertainty in the imaginary part of Z MS h is basically of the order of its magnitude, as it is expected to be real. Note that ignoring the imaginary part does not provide any solution, as they need to come out zero from the computation in order to claim that the computation is fully reliable. Thus, conversion truncation effects become the most major source of uncertainty of the imaginary parts of the Z-factors.
Based on the study presented in this subsection, we present in Table 2 a summary of our quantitative estimates on the systematic uncertainties present in Z MS h . The estimate of the total effect is not a simple sum of isolated effects, but takes into account their different signs discussed above. For the imaginary part, the estimates are for |z/a| ≥ 5, since Im [Z MS h ]∼0 for smaller values of z and the relative effects become meaningless.
We would like to stress that due to the complex multiplication of the Z-factors and the bare matrix elements, the large uncertainty in the imaginary part of the Z-factor implies that also real part of the renormalized matrix elements is affected. Furthermore, the uncertainties that we consider in the Z-factors translate differently to the uncertainties of the renormalized matrix elements for different Wilson line lengths and thus, we postpone the discussion of this influence to the next subsection.

Renormalized results
Once the Z-factors are obtained for the unpolarized, helicity and transversity quasi-PDFs, one may proceed with the application of the renormalization in the nucleon matrix elements. Here we mostly focus on the helicity case, as the renormalization is multiplicative. The case of the unpolarized quasi-PDF is briefly mentioned in the end of the subsection.
In Fig. 9, we show the renormalized helicity nucleon matrix elements, and compare them with the bare ones. This is a straightforward procedure as there is no mixing for the axial operator and, therefore, the renormalization is only multiplicative: The above formula involves complex quantities, and thus, h MS (z) is a mixture of the real and imaginary part of the bare matrix element. However, each value of z is renormalized independently. We use Z MS h from the extrapolation of Eq.  The renormalized matrix elements of the helicity quasi-PDF are presented here for the first time, and they demonstrate the enormous progress in the field of the quasi-PDFs. However, before attempting to compare with the physical PDFs, we must understand the uncertainties that are inherited to h MS (z) from its renormalization function. As we argued in the previous subsection, a robust computation needs the subtraction of O(g 2 a ∞ ) lattice artifacts and a significant reduction of truncation uncertainties in the conversion between the RI and MS schemes. We attempt setting bounds on these systematics, starting with the real part of renormalized matrix elements, which reads: The uncertainty in the small-z region of renormalized matrix elements is dominated by uncertainties in the real part of the Z-factor, which, as argued in the previous subsection, should not exceed 5%. The local minimum observed at z = 0 is likely to result from lattice artifacts and truncation effects. In the large-z region, the real part of the renormalized matrix elements receives negative contributions from the imaginary part of the bare matrix element if Z MS h has a non-zero imaginary part and Im [ h bare ] decays to zero more slowly than Re [ h bare ], which is what we observe in the data. This effect is unphysical and is expected to be strongly suppressed or eliminated by extending our calculation to the two-loop order in the conversion and by subtracting lattice artifacts.
The imaginary part of renormalized matrix elements, , which is of the order of 5% and comparable to the currently attained statistical uncertainty. Thus, h MS is rather robust in this region. However, when z is increased, the uncertainty from Im [Z MS h ] starts to dominate and reaches 100% for large z. In this way, the improvements expected by the perturbative subtraction of O(g 2 a ∞ ) lattice artifacts and the extension to the two-loop perturbative conversion to the MS scheme are very important for obtaining meaningful values of renormalized matrix elements and hence, also quasi-PDFs.
The values of the multiplicative vector renormalization factor and the mixing coefficient can be used to properly renormalize the unpolarized quasi-PDF through Eq. (12). The successful renormalization requires the bare nucleon matrix elements for the scalar and vector operators, in order to extract the renormalized h RI V . Then, it must be multiplied by the conversion factor C V to bring the results to the MS scheme. Note that once the mixing is treated, the conversion factor is multiplicative and not a 2×2 matrix. This is due to the fact that no mixing is present in the continuum dimensional regularization.

Matching to light-front PDFs
Having the renormalized matrix elements, one can perform a Fourier transform and obtain the renormalized quasi-PDF, which represents the distribution of quark momenta for a finitemomentum nucleon moving in the chosen spatial (z) direction: where x is the quark momentum fraction. The quasi-PDF, expressed in our case in the MS scheme at μ = 2 GeV, can then be connected to the light-front PDF, in the same scheme and at the same scale, using one-loop perturbative matching. The matching procedure uses the fact that only the ultraviolet physics is different in quasi-(q) and light-front (q) PDFs [11]. Hence, the one-loop difference between them is expressed as the difference between vertex corrections (denoted Z (1) below) and wave function corrections (δZ (1) ) for the finite momentum and infinite momentum cases. The generic matching formula for quasi-PDFs is: [50,12] q (x, μ) =q(x, μ, where α s (μ) is the strong coupling constant at the scale μ. In the integral, we exclude a small region around ξ = 0, such that the argument of q is not too large, as we would then pick up contributions from the mirror images of q resulting from the periodicity of the Fourier transform in Eq. (23). The linearly divergent terms ∝ /P 3 ( : transverse momentum cutoff) of the matching formulae from Ref. [11] are not present in our renormalized results. For purely orientational purposes, we also plot phenomenological PDFs (DSSV08 [51] and JAM15 [52]). However, we emphasize that no quantitative comparison with our results is aimed at, since careful consideration of a number of systematic effects is still needed. These include: cut-off effects, non-physical pion mass, finite volume effects, possible contamination by excited states, extrapolation to infinite nucleon boost, as well as the improvements in the computation of MS renormalization functions, postulated in the previous subsection.
In Fig. 10, we show the matched helicity PDF computed with either fully renormalized matrix elements obtained in this work (blue) or with bare matrix elements multiplied by the local (z = 0) axial vector current renormalization function Z A (magenta), corresponding to our results from Ref. [17]. We observe that the renormalized matrix elements from this work move towards the phenomenological PDFs, which is promising. In particular, the antiquark asymmetry is not overestimated any longer and actually this asymmetry becomes compatible with zero under current uncertainties. In the quark part, there is an enhancement of the matched PDF for all values of x. We emphasize again that the comparison with the phenomenological PDFs should be understood as only qualitative. For quantitative comparison, a careful investigation of a number of systematic effects is still needed. These include: cut-off effects, non-physical pion mass, finite volume effects, possible contamination by excited states, extrapolation to infinite nucleon boost, as well as the improvements in the computation of MS renormalization functions, postulated in the previous subsection: subtraction of lattice artifacts computed in lattice perturbation theory and reduction of truncation effects in the perturbative conversion to the MS scheme.

Conclusions and discussion
In this work we have presented a concrete prescription to renormalize non-perturbatively the matrix elements needed for the computation of quasi-PDFs. The employed scheme is RI , which is then converted to the MS scheme and evolved to 2 GeV; this is done perturbatively to one-loop. We have argued that the renormalization condition properly handles both kinds of divergences present in the matrix elements: the standard logarithmic divergence and the power divergence specific to non-local operators containing a Wilson line. Furthermore, we provide the renormal-ization conditions to eliminate the mixing in the case of the unpolarized quasi-PDF that mixes with the twist-3 scalar operator.
We have also demonstrated the implementation of the proposed prescription to the helicity quasi-PDF and presented the corresponding renormalized matrix elements. This has allowed us to draw conclusions how to make the computation more robust, which is the main outcome of this work.
• First, an essential ingredient of a computation with controlled systematic uncertainties is the subtraction of one-loop lattice artifacts in the framework of lattice perturbation theory. Following the ideas of Ref. [38], we are currently computing the O(g 2 a ∞ ) artifacts that will be subtracted from the non-perturbative estimates for the Z-factors. In this way, the presence of large cut-off effects in the renormalized functions (especially for "parallel" momenta) will be avoided to a large extent. • Second, the conversion factor from the RI renormalization scheme to the MS scheme is likely to have sizable higher order corrections that, among others, are responsible for the unphysical feature of the real part of the renormalized matrix element becoming negative for large Wilson line lengths. A two-loop computation of this conversion factor is expected to resolve this issue to a sufficient degree. We have performed numerical experiments that indicate that a natural change of the conversion factor by two-loop contributions, i.e. around 10-20% (which is approximately α s at the considered scale), should be enough to suppress the unwanted effect. A perturbative calculation of the conversion factor to two loops is quite laborious and will be presented separately.
To summarize, the renormalization program presented in this work together with future improvements that are being pursued, will provide reliable estimates for the renormalization functions of the Wilson line fermion operators. In this fashion, the obtained renormalized matrix elements can be used as an input to calculate the quasi-PDFs and match them to light-front PDFs, which is the main aim of the whole approach. Apart from the helicity case discussed in this work, we will address the transversity PDF in an analogous manner. For the unpolarized case, one needs to take into account the mixing with the scalar operator, as explained and numerically demonstrated here. With this work, we have proposed and discussed a complete renormalization program of the quasi-PDFs, which has been a major uncertainty prior to this work and constitutes a crucial milestone in connecting lattice QCD results to the light-cone PDFs.