Self-Renormalization of Quasi-Light-Front Correlators on the Lattice

In applying large-momentum effective theory, renormalization of the Euclidean correlators in lattice regularization is a challenge due to linear divergences in the self-energy of Wilson lines. Based on lattice QCD matrix elements of the quasi-PDF operator at lattice spacing $a$= 0.03 fm $\sim$ 0.12 fm with clover and overlap valence quarks on staggered and domain-wall sea, we design a strategy to disentangle the divergent renormalization factors from finite physics matrix elements, which can be matched to a continuum scheme at short distance such as dimensional regularization and minimal subtraction. Our results indicate that the renormalization factors are universal in the hadron state matrix elements. Moreover, the physical matrix elements appear independent of the valence fermion formulations. These conclusions remain valid even with HYP smearing which reduces the statistical errors albeit reducing control of the renormalization procedure. Moreover, we find a large non-perturbative effect in the popular RI/MOM and ratio renormalization scheme, suggesting favor of the hybrid renormalization procedure proposed recently.


I. INTRODUCTION
Parton distribution functions (PDFs) provide an effective description of quarks and gluons inside a lighttravelling nucleon [1,2]. They play an essential role in calculating hardronic cross sections involving the nucleons [3], and their uncertainties from phenomenological extractions have been one of the major sources of errors in the theoretical predictions at Large Hadron Collider (LHC) and other hadron facilities. Thus, a precise knowledge of PDFs is very important both for the accurate tests of the Standard Model (SM) and for the search of new physics beyond the SM. On the other hand, the densities of partons in the nucleon provide direct information on its intrinsic properties such as the origin of the * yh3285@columbia.edu † ysu12345@umd.edu nucleon spin and mass, as well as the role of sea quarks for various physical quantities [4]. However, directly calculating parton physics from first principles of quantum chromodynamics (QCD) has been a difficult task. A review on various efforts of doing so can be found in Ref. [5]. The proposal of the large-momentum effective field theory (LaMET) has been an important step toward meeting the challenge [6][7][8]. So far, LaMET has been widely used in calculating quark isovector distribution functions [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], generalized parton distributions [25,26], distribution amplitudes (DAs) [27][28][29], and transverse-momentum-dependent distributions [30][31][32].
LaMET suggests to calculate time-independent physical distributions in a finite momentum nucleon, which are Euclidean observables. Such finite-momentum quantities can then be matched to light-front (LF) parton properties using effective field theory techniques [33]. For example, for the collinear quark distributions, it has been suggested to first compute a Euclidean space correlation arXiv:2103.02965v1 [hep-lat] 4 Mar 2021 (or quasi-light-front (quasi-LF) correlation) on the lattice,h (z, P z ) = P |O γt (z)|P , (1) where |P is the nucleon state with a large momentum P and the non-local operator is where ψ,ψ denote the bare quark field, Γ is a Dirac structure, and U (0, z) = exp(−ig z 0 dz A z (z )) is the Wilson link along the direction z, where A µ is the gluon gauge potential. After renormalizingh properly and matching it to some continuum scheme such as dimensional regularization and (modified) minimal subtraction (MS), one obtains the so-called quark quasi-PDF via a Fourier transform, from which the quark PDF can be extracted through perturbative QCD matching.
Renormalizing the quasi-LF correlationh under lattice regularization is a non-trivial task. To see this, let us take, as a simple example, the matrix element of O Γ (z) in an off-shell quark state |q with large Euclidean momentum p 2 (actually an off-shell truncated Green's function). It has the following form at 1-loop level [34], where g is the bare gauge coupling and a is the lattice spacing. Unlike the coefficient of the logarithm γ, the coefficient of the linear divergence m −1 can be sensitive to the details of the fermion and gauge actions on the lattice as those actions themselves define the lattice regularization. The linear-divergence term is proportional to the Wilson link length z/a in lattice units, and it can be exponentially large at large z or small a, when higher-order corrections are included. Thus, the linear divergence effect has to be removed before the lattice calculated quasi-PDF can be extrapolated to the continuum limit a → 0 and eventually matched to the continuum-scheme PDF. Theoretical studies so far have shown that the quasi-PDF operator is multiplicatively renormalizable [35][36][37][38] in a continuum theory. On the lattice, due to noncommutativity of the limit z → 0 and a → 0, it has been suggested recently [39] to use a hybrid scheme to renormalize the large and short distance correlations separately, which has the advantage of avoiding certain discretization effects and undesired infrared effects introduced in the renormalization stage. At small z where a short distance expansion is valid, one can use various ratio schemes, such as dividing by an offshell quark matrix element of the quasi-PDF operator in the regularization-independent momentum subtraction (RI/MOM) scheme [34,38,40] or by a hadron matrix element in the rest frameh(z, P z = 0) [41,42]. At large distance, one can directly remove the Wilson line self-energy effect by subtracting the power divergence [36][37][38]43]. Apart from using the RI/MOM factor or the rest-frame hadron matrix element, other methods have been suggested to extract the linear divergence, including Wilson loop [27,43], vacuum expectation values O Γ (z) [44], and gauge fixed Wilson link [8], and so on [39].
However, numerically subtracting the linear divergence is an extremely delicate exercise. First, linear divergence could be sensitive to computational systematics in lattice calculations. In data, there may be slight differences between two matrix elements with the same linear divergence. These differences may lead to a failure of the suggested methods, especially for small lattice spacings. Second, some of the renormalization factors have other problems. For example, the leading contribution to the vacuum expectation value of O Γ (z) at short distance vanishes and therefore, it is numerically very challenging to obtain the relevant linear divergence. Finally, as we shall see, linearly-divergent chiral symmetry breaking effects for Wilson fermions may render the linear divergence non-universal [45].
To understand better the effects of the linear divergence in LaMET applications, we study systematically the linear divergences of matrix elements for sets of lattice data, generated with different lattice actions. We propose a self-renormalization method to eliminate all divergences and discretization errors when data for several different lattice spacings are available. The idea of this method is to extract the renormalization factor and the residual contribution directly from the matrix element we want to renormalize, without using an additional matrix element for renormalization. We then match the empirically renormalized matrix elements to those in the continuum. Our method is largely model-independent, and each term of our fitting functions is motivated by physics considerations. Tests show that our method works well for all data sets considered, which include matrix elements calculated with the lattice spacings from 0.03 fm to 0.12 fm using the valence clover and overlap actions on the MILC [46] N f = 2 + 1 + 1 and RBC [47] N f = 2 + 1 configurations. We also find that the method works for matrix elements after applying smearing which is needed to improve the statistical precision.
The rest of the paper is organized as follows. We present the theory of perturbative renormalization in Sec. II. The definition of the matrix elements we calculate and the simulation setup can be found in Sec. III. In Sec. IV the linear divergence is analysed for the different ensembles. Sec. V presents our strategy to selfrenormalize a matrix element and collects the results for the O γt (z) matrix elements in different states and calculated with different valence fermion actions. Sec. VI discusses other fitting options, which mainly differ in the treatment of higher-order terms. In Sec. IV, V, and VI, we only analyze matrix elements without HYP smearing. In Sec. VII, we extend our method to HYP smeared cases. In Sec. VIII, we test the linear divergence in some vacuum state matrix elements, including Wilson loop, quasi-PDF operator, and gauge fixed Wilson link, for the HYP smeared cases. In Sec. IX, we summarize the results and discuss issues to be addressed in further studies.

II. RENORMALIZATION IN PERTURBATION THEORY
According to the standard renormalization in local quantum field theories [48], the renormalization of the matrix elements of an operator does not depend on the external states but is only related to the short-distance property of the operator itself. Therefore, one can study the renormalization property of the operator in perturbative Green's functions, or off-shell quark and gluon matrix elements. For LaMET applications to PDFs, we are interested in the non-local operator O Γ (z) = ψ(0)ΓU (0, z)ψ(z) in Eq. (2). There are two types of divergences: The linear divergence associated with the Wilson link (its self-energy) and the logarithmic divergences associated with the renormalization of the vertices involving the Wilson line and light quark. The renormalization is multiplicative, and only the linear divergence has a (linear) z dependence. Therefore the renormalized operator is [36][37][38][39] where e δmz contains the linear divergence and Z O the logarithmic ones. δm is not uniquely defined apart from the linear divergence, introducing a subtraction scheme dependence which affects the z-dependence of the renormalized operator [39]. The linearly-divergent mass renormalization can be calculated in perturbation theory. At one-loop order, the result is independent of lattice action [43], The energy scale of α s can be chosen as 1/a to match the lattice results, where b 0 = 11 − 2 3 n f (we will take n f =3) is the QCD β function at one-loop order with n f species of fermions. Λ QCD is the non-perturbative QCD scale. Including higher-orders in the β function will lead to a more complicated expression. Different choices of energy scale amount to including high-order corrections. Λ QCD and higher-order α s terms in δm also depend on the lattice action.
The perturbation theory does not converge due to infrared renormalons, see from example [49][50][51][52]. The existence of the renormalons signals a non-perturbative term in δm which is independent of a, where the minus sign is just a convention. The uncertainty in summing the perturbation series to get m −1 (a) is compensated by the same uncertainty in the nonperturbative m 0 , leaving the total independent of the renormalon uncertainty. Additional uncertainty in m 0 comes from the subtraction scheme, or equivalently from the matching to the continuum scheme. To reduce the subtraction scheme dependence, we require the renormalized lattice correlation function to be consistent with the MS result from continuum perturbation theory at short-distances z. To be more concrete, we determine m 0 by matching the lattice result to the MS one within a window a z < z S , where z S < 1/µ is the point beyond which perturbation theory ceases to work and µ is a perturbative scale. The condition a z ensures that the discretization effect on lattice results is small. For this special choice m 0c , the LaMET expansion does not have a linear term in 1/P z [39].
At one-loop order and in dimensional regularization, the renormalization factor of the logarithimic divergence is [8,53], where C F = 4/3 is the Casimir operator for the fundamental representation of SU (3) and d is the space-time dimension. One can resum the logarithimic divergence through solving the renormalization group equation, where = (4 − d)/2 and γ = − 3C F 2π α s (µ, Λ QCD ) is the leading anomalous dimension, which is independent of the regularization method. The form of the leadingorder solution of Eq. (9) is independent of regularization scheme and, on the lattice, is where the bare ultraviolet (UV) cut-off is 1/a, and the renormalization scale is µ. If one considers the contributions from sub-leading logarithms to the anomalous dimension γ when solving the renormalization group equation Eq. (9), we obtain [54] where Λ QCD and the constant d depend on the specific lattice action.

III. LATTICE MATRIX ELEMENTS AND SIMULATION SETUP
The standard non-perturbative renormalization of lattice matrix elements is through calculating some auxiliary matrix elements of the operator (which can also be computed in perturbation theory) and using them as the renormalization factor to approach the continuum limit. In this paper, we focus on the study of the auxiliary matrix elements potentially useful for the renormalization of quasi-LF correlations. We mainly consider the matrix elements of the quasi-PDF operators in the following cases: 1) in a large Euclidean momentum quark state, 2) in the physical zero-momentum state of a hadron such as the pion or nucleon. We will also study the matrix element in the vacuum [44] (case 3), and show that it is not a good choice for renormalization because at small z such a matrix element is suppressed as can be shown by operator product expansion (OPE). Moreover, since the zdependent renormalization is mainly associated with the Wilson link, we shall also consider the matrix elements of Wilson loop as well as the Wilson link in a fixed gauge (case 4).
We calculate the pion and nucleon matrix elements h 0 H (z) = H|O γt (z)|H P =0 in the rest frame evaluating the following three-point function, where O γt (z, ( x, t)) =ψ( x, t)γ t U (( x, t), ( x+n z z, t))ψ( x+ n z z, t),n z is the unitary vector along the z direction, and O H is the interpolation field of a hadron such as a pion and a nucleon. t 2 corresponds to the source-sink separation.
To obtainh 0 H (z) accurately, we need both t and t 2 − t to be large enough to suppress excited-state contaminations, or fit the three point function with a proper parametrization on the excited-state contaminations. We can use the first strategy for the pion case since the signalto-noise ratio will not decay for the pion in the rest frame; while the second strategy is essential for the other cases where the statistical uncertainty increases exponentially with t and t 2 . Due to (anti)periodic boundary conditions we can use at most t 2 = T /2 which is larger than 2 fm on most modern lattice ensembles. Since the mass gap δm ∼ 1 GeV between the ground and first exited state of pion h0 π (z) can be extracted with sufficient accuracy. At small perturbative z, h0 π (z) −1 acts as a renormalization factor up to certain O(Λ 2 QCD z 2 ) power corrections [41]. A common non-perturbative renormalization method for lattice QCD matrix elements with local operators is the RI/MOM scheme [55] and its modified versions [56]. One can calculate the bare matrix element with given operators in an off-shell quark state, for both the lattice and dimensional regularizations, and renormalize their difference in the MS scheme to get the renormalization factor for the bare quantities. For the quasi-PDF, such a renormalization constant (the RI/MOM renormalization factor) is defined through the O γt (z) matrix element in the off-shell quark state with given momentum p 2 [34,38,40,57]: where N c = 3 is the number of colors, µ R is the RI/MOM renormalization scale. We will take µ R = 3 GeV since the RI/MOM factor has little dependence on the renormalizaiton scale in a certain range [45]. We work with Landau gauge fixed configurations, and choose p z = p t = 0 to eliminate the difference between the projections [16]. Factor Z RI can be calculated non-perturbatively for any lattice regularization, or equivalently for any quark and gluon action defined on the lattice. However, the current lattice QCD calculation is limited to the Landau gauge. The perturbative matching with off-shell states in Landau gauge can be complicated beyond one-loop level [58].
The z-dependent part of renormalization comes from the Wilson line, therefore one may extract the linear divergences from matrix elements of pure Wilson lines. First one can consider a Wilson loop [27,28,38,43,59], Ref. [43] proposed this scheme to renormalize the pion DA, and the linear divergence coefficient was extracted precisely based on the calculation with several lattice spacings in the following Kaon DA study [27]. But most of subsequent lattice QCD calculations have switched to the RI/MOM scheme or used a hadronic matrix element.
One can also consider the matrix element of a single Wilson Link U (0, z) in a fixed gauge such as the Landau gauge. As the simplest choice without any external state, U (0, z) can provide a reference to identify whether the linearly divergent behavior is sensitive to the existence of the external state, as suggested in the multiplicative renormalizability studies of the quasi-PDF operator [35][36][37][38]. Our lattice calculations are performed using the Chroma software suite [60] and QUDA [61][62][63] in the HIP programming model [64]. We use 2+1+1 flavors (degenerate up and down, strange, and charm degrees of freedom) of highly improved staggered quarks (HISQ) [65] ensembles generated by the MILC Collaboration [46] at five lattice spacings, and 2+1 flavor domain wall (DW) quarks and Iwasaki gauge ensembles from the symbol 6/g 2 L T a(fm) DW11 2. 13  RBC/UKQCD collaboration [47] at three lattice spacings. The pion mass for all ensembles are tuned to be roughly 310 MeV based on the pion mass of the light sea quark mass on the corresponding ensemble. The lattice spacings for the MILC ensembles are determined using Wilson flow based on the parameters determined by Ref. [66].
We use the matrix elements without hyper-cubic (HYP) smearing in order to test the perturbative renormalization analysis (Sec. IV, V, and VI). However, since in practical calculations the results without smearing can be rather noisy, it is standard to use some type of smearing in lattice simulations. It is unclear to us how strongly the smearing will interfere with renormalization. Therefore, we regard our method as a purely phenomenological approach for the time being to analyze data with one step HYP smearing [67] in Sec. VII and VIII, hoping that the gain in statistical precision outweighs the additional systematic uncertainty from moderate smearing.
We use two kinds of fermion actions for the valence quarks: clover and overlap fermions. Clover fermions break chiral symmetry and the action is defined by wheren µ is the unit vector along the µ direction, the clover coefficient c sw is the tadpole improved tree level value, and m 0 q a is the bare quark mass which is determined by requiring the pion mass to be roughly 310 MeV. Parameters c sw and m 0 q a should approach 1 and 0 respectively in the continuum, while the lattice spacing dependence is weaker than the O(a) discretization effect in the present range of a, and closer to that of the gauge coupling g 2 as predicted by lattice perturbative theory.
The overlap action preserves chiral symmetry but the simulation is rather expensive. We use it to test the dependence of renormalization on the fermion action. The overlap fermion action is defined by [68,69] S ov q = x,yψ , where and −ρ should be smaller than the bare quark mass for which vanishes the pion mass to make D ov to be the same as the standard Dirac operator in the continuum limit. We choose −ρ = 1.5.
Although different active and sea fermion formulations will in general introduce non-unitarity, the effect shall vanish in the continuum limit. However, at finite a, the difference will generate systematic uncertainties which can affect the final results.
We start with nucleon and pion matrix elements for the MILC ensembles and clover action. For the pion matrix element, We average the R π (t 2 , t, z) data with t 2 = T /2 and t ∈ [T /8, 3T /8] to get a precise estimate of the ground state matrix element, and it should be also precise since T is between 7.7 and 9.1 fm in the ensembles we used, as shown in the upper panel of Fig. 1. Note that R π (T /2, T /4, z) equals to π|O γt (z)|π /2 in such a case. The additional factor 1/2 corrects for the fact that there are forward and backward propagating states for (anti)periodic boundary conditions. The nucleon case is known to be much noisier especially at large t 2 . Thus, we just consider t 2 = 2t 0.72 fm in all cases, and plot R N (t 2 , t, z) in the lower panel of Fig. 1.
The normalized bareh 0 N (z)/h 0 N (0) whereh 0 N (z) is approximated by R N (t 2 , t, z) with t 2 = 2t 0.72 fm is plotted in Fig. 2, with a normalization factor 1/h 0 N (0) using jackknife resampling to make it to be exactly one at z = 0 at all lattice spacings. As in the figure, the linear divergence makesh 0 N (z)/h 0 N (0) decay exponentially with both z and a, and thus one cannot obtain any meaningful continuum limit when a → 0. A renormalization is necessary to recover a good continuum limit.

IV. SIMPLE TEST OF LINEAR DIVERGENCE
To study the renormalization properties of the quasi-PDF operator, we need to calculate matrix elements at several different lattice spacings. To ensure the data to be useful for a refined high-precision analysis, we start by testing whether they approximately show the linear divergence predicted by perturbation theory. In particular, one needs to show that the z-dependence of the linearly divergent term is linear.
To achieve this, we first extract the factor e −δmz from the bare matrix element to test the 1/a dependence in δm in Eq. (5). Based on Eqs. (4), (5), and (6), if we take the natural log on a bare matrix element M, we have,   where the first term is linearly divergent and g(z) is the residual. We have ignored the logarithmic divergence and discretization error here because these terms are much smaller and thus have little influence on this simple test of the linear divergence.
We use Eq. (18) to fit the bare matrix elements as a function of a. We treat e(z) and g(z) as unknown functions of z so that our fitting is performed at each value of z at which they are treated as free parameters. We need to do linear interpolation on the data ln M(z, a) with respect to z. Here are the steps we perform in detail: 1. We do the linear interpolation with neighbourhood two data points to obtain the central value and uncertainty at z = 0.06 × n fm (n = 3, 4, 5, ..., 16).
2. For each z, we fit the dependence on a, treating e(z), g(z) as fitted parameters. For this initial test, we fix Λ QCD at 0.2 GeV, which is a reasonable guess to start with [70].
3. We plot the fitted e(z) with respect to z to see if e(z) has a linear z dependence.   If we obtain a reasonable fit in step 2, we then see approximately the 1/a dependence of the divergences. In step 3, we can test the linear z dependence. The fitting results for the bare O γt (z) matrix element in the offshell quark state for valence overlap and clover actions without HYP smearing are shown in Fig. 3 and Fig. 4, respectively. There are eight different lattice spacings. The range of z is taken from 0.18 fm to 0.96 fm and we analyse 14 different z values in this range.
The 1/a dependence in both cases is fitted very well, although due to the high precision of the data, chi-square is large. However, this is not our concern at this initial step. The z dependence of the divergent coefficients shows nicely the linear feature, approximately going through zero at z = 0. This is a strong indication that the linear divergence follows roughly the prediction of perturbation theory which gives us confidence to perform more refined analyses in the next sections.

A. The Strategy
The non-perturbative matrix elements contain the following important contributions: a) the linear divergence, b) a finite term coming from non-perturbative renormalon physics, and c) the residual contribution encoding intrinsic non-perturbative physics. It is part c) that is required to extract the partonic structure information. Thus, it is important to separate out the latter from the former in a systematic way and with high precision.
Here we develop a self-renormalization method to extract the residual intrinsic non-perturbative physics and renormalization factor from a matrix element itself without using a different matrix element. The extraction of the residual is much harder than the extraction of the linear divergence term because after we take the logarithm of the matrix element (Eq. (18)), the linear divergence term is dominant and the residual is very sensitive to the subtraction. We need to properly take into account the fine features of the data, including logarithmic divergences and discretization errors.
Based on Eq. (11), we modify the fitting function Eq. (18) to be, where the first term is linearly divergent. g(z) is the residual, which contains the non-perturbative m 0 effect and the intrinsic non-perturbative physics. f 1,2 (z)a takes into account the discretization errors, which we allow to be different for the data calculated from MILC and RBC ensembles (f 1 (z) for MILC and f 2 (z) for RBC). The last two terms come from the resummation of leading and sub-leading logarithmic divergences, which only affect the overall normalization at different lattice spacings. We use Eq. (19) to fit the logarithms to the data for the bare matrix elements for each z. The z points can be chosen in a large range where the lattice discretization error is small and, at the same time, the statistical error is limited. We fit the dependence on a, treating g(z), f 1 (z), f 2 (z) as free parameters and fixing µ at 2 GeV, which is the renormalization scale we choose to use. We will discuss k and Λ QCD in detail in Sec. V B and d in Sec. V C. We obtain the residual g(z) from the fit.
The g(z) obtained this way does not correspond to what is in the perturbative MS scheme because it contains a non-perturbative m 0 effect. We need to eliminate this effect through matching the lattice result to the continuum scheme at short distance z. Based on Eq. (4) and Eq. (7), we use the following equation within a window a z < 1/µ to extract m 0 , where Z MS is the perturbative matrix element in the MS scheme. For the O γt (z) matrix elements in the pion, nucleon, and off-shell quark state, we take Z MS at one loop [42], from operator product expansion, where we fix µ = 2 GeV and Λ MS = 0.3 GeV. The renormalized physical matrix element is then obtained from and is expected to be consistent with Z MS at small z.
One could consider higher-order MS matrix elements for matching, including summing over large logarithms. We choose not to do this here as it does not affect the procedure for the subsequent analysis, and therefore does not change the main conclusions of the paper. However, to obtain physical results at high precisions, one might need to examine this more carefully.
We thus define a renormalization factor which includes the discretization error as well. All the parameters here are either fixed, fitted or fine-tuned through the renormalization procedure from the matrix element M(z, a) itself. Dividing the bare matrix element M(z, a) by Z(z, a) R , we get If our procedure eliminates all the divergences and discretization errors, we should expect that M R (z) is aindependent and the same as M(z) R (Eq. (22)). The degree to which this is the case can be used as a test of the self-renormalization procedure.

B. Renormalon Uncertainty
Now let us turn to the parameter k and Λ QCD in Eq. (19). In Eq. (19), we have only considered the leading order perturbative term for the linear divergence and neglected higher-order ones. These can be included partially by a proper choice of Λ QCD . For example, the coupling constants with different choices of Λ are related to each other by perturbation theory [70],  In principle, Λ QCD depends on the specific lattice action, and one needs to perform multi-loop calculations to relate them. In our analysis, we treat Λ QCD as a fitting parameter which, therefore, can effectively take into account some higher-order effects.
Next we want to discuss the coefficient k of the linear divergence. QCD perturbation theory not only predicts the 1/a dependence and linear z dependence (which we have tested in Sec. IV), but also the value of k. Comparing Eq. (5) and Eq. (19), we get the one-loop perturbative value of k: ] for different sets of (k, ΛQCD) along the "small-χ 2 band".
Next, we check whether this value is consistent with our data.
We use Eq. (19) to fit the bare O γt (z) matrix element in the pion state for the clover action without HYP smearing, treating g(z), f 1 (z), f 2 (z) as fitted parameters. If we take (k, Λ QCD ) to be one set of values, e.g., (7.4 GeV −1 fm −1 , 0.09 GeV), we can perform a fitting for each z, as seen in Fig. 5. There is a χ 2 /d.o.f. for the fit at each z and we calculate the average of them, denoted as χ 2 /d.o.f. z . Varying (k, Λ QCD ) generates a χ 2 map, as seen in Fig. 6.
The sets of (k, Λ QCD ) of small χ 2 lie in a band, which we call the "small-χ 2 band". That means that k and Λ QCD are strongly correlated, while there is a large uncertainty for k and Λ QCD separately. If we choose some sets of values from the "small-χ 2 band" to do the fitting, we find the fitted residuals Exp[g(z)] are very different, as seen in Fig. 7. However, after eliminating the nonperturbative m 0 effect, the renormalized matrix elements Exp[g(z) − m 0 z] are all the same (see Fig. 8). Therefore, the uncertainty in k and Λ QCD will not significantly influence the final physical result.
The reason for this uncertainty lies in the behavior of k a ln[aΛQCD] (the term related to the linear divergence in Eq. (19)) in the a range of our data. Varying along the "small-χ 2 band" effectively shifts k a ln[aΛQCD] by a constant C, see Fig. 9. The Cz term is absorbed into g(z) automatically during fitting so there is a large uncertainty for k and Λ QCD . However, when we eliminate the m 0 effect, the Cz term has been taken into account, and the final result does not change.
Another way to interpret this is that the perturbative value of k is consistent with the "small-χ 2 " band and therefore is consistent with the fit. As a consequence, we can fix k at the one-loop perturbative value, for different sets of (k, ΛQCD) along the "small-χ 2 band".
FIG. 9. The linear divergent term in the a range of our data for different sets of (k, ΛQCD) along the "small-χ 2 band". still vary within a reasonable range of χ 2 while the physical result is stabilized by the choice of m 0 . We call this correlation between Λ QCD and m 0 the "phenomenological renormalon uncertainty", in the sense that Λ QCD is characteristic for the different truncation/resummation schemes for the perturbation series, but that the nonperturbative m 0 compensates the effects of the differences [49][50][51][52].

C. Tuning Parameter-d to Match the Continuum Scheme
In principle, the parameter d in Eq. (19) is determined by the lattice action. In our data, there are two different dynamical lattice ensembles (MILC and RBC) and two different valence fermion actions (overlap and clover). Here we treat d as a fine-tuning parameter to make sure that tional to z within a window a z < 1/µ. As shown in Fig. 10, this is quite effective. We have also tested that d has little influence on χ 2 since the effect of tuning d shifts g(z) simply by a constant. Fig. 11 shows that the renormalized matrix element Exp[g(z)−m 0 z] (blue) for clover quarks is consistent with Z MS (black) at small z. However, at large z, there is a significant discrepancy between Exp[g(z) − m 0 z] and Z MS . That means that there is a significant non-perturbative effect at large z in the popular RI/MOM and ratio renor- malization schemes used previously. To our knowledge, this is the first time that this difference has been studied in the literature. This finding supports the usage of the hybrid renormalization procedure proposed recently [39].

D. Self-Renormalization Procedure
The procedure developed in the previous subsections forms our self-renormalization strategy. This procedure can overcome some of the problems encountered in the renormalization of the linear divergence due to the numerical uncertainties associated with exponentiallyamplified small errors. There are three steps in this process: 1. Use Eq. (19) to fit the data of the bare matrix elements M. For each z (here z can be chosen in a large range and linear interpolation with respect to z might be needed for the data in lnM(z, a) ), fit the dependence on a, treating Λ QCD as a global parameter, namely the same for all z, and g(z), f 1 (z), f 2 (z) as fit parameters, with fixing k=7.4 GeV −1 fm −1 (or 1.46 if dimensionless), and µ=2 GeV. Fine-tune d to make sure that g(z) − ln[Z MS ] is proportional to z within a window a z < 1/µ. We obtain the residual g(z) through fitting; 2. Use Eq. (20) to fit the dependence on z (within a window a z < 1/µ) to extract m 0 . Then calculate the renormalized matrix element Exp[g(z) − m 0 z] (Eq. (22)), which is considered valid for a large range of z; In step 1 and 2, one can use another way to get d and m 0 : Modify Eq. (19) through replacing g(z) with (ln[Z MS (z, µ, Λ MS )] + m 0 z), and then use the modified fitting function to fit the bare matrix element at small z to extract d and m 0 . This way should be equivalent to the procedure outlined above and we will not follow it here.
In step 3, we need to calculate Z(z, a) R as well as estimate its error. We shall not propagate the fitting error of Λ QCD to Z(z, a) R because the effect of slight variation of Λ QCD can be compensated by m 0 and will not influence our result as we have discussed in Sec. V B.
Since the lattice data we use do not provide systematic errors, we will add a dummy systematic error δ sys to the data during our fitting [45] where we assume that there is an aµ dependence as well. This choice is intended to give more weight to the small lattice spacing data in the fitting.

E. Comparing Results from Different Actions and Matrix Elements
In this subsection, we apply the self-renormalization method to the O γt (z) matrix elements in the pion, nucleon, and off-shell quark state with valence clover and overlap actions without HYP smearing. We compare the renormalization factors and physics results. We find that apart from the case of the RI/MOM matrix element with clover valence fermion, all other renormalization factors are similar. Despite this difference in renormalization factor, the physical results are independent of fermion actions.
Figs. 12 through 16 show the detailed results of applying the self-renormalization method for different cases. For the RI/MOM factor with overlap and clover fermions, we analyse eight different lattice spacings from two different types of gauge ensembles. For the pion matrix element calculated with the overlap fermion action, we analyse three different lattice spacings from MILC ensembles. In the clover pion case, we use six different lattice spacings from two different types of gauge ensembles. For the nucleon matrix element with clover fermion action, we analyse five different lattice spacings from MILC ensembles.
Subfigure (a) in each figure is used to estimate the best fitted value of the global parameter Λ QCD as well as its error. For each chosen Λ QCD , we can use Eq. (19) to fit the bare matrix element. The fitting for each z gives us a separate χ 2 . We can calculate the averaged χ 2 for different z, denoted as χ 2 z . Subfigure (a) shows gives us the best fitted Λ QCD and χ 2 z − [ χ 2 z ] min =1 gives us its error. If the data points for different z were independent, we should use the total χ 2 for the fits at different z to estimate the error of Λ QCD . However, since we have done linear interpolations of the original data points, some of the points are correlated to each other. So here we just use χ 2 z to estimate the error of Λ QCD . The range of z starts from 0.18 fm here but not 0.06 fm because the χ 2 for z = 0.06 or 0.12 fm is much larger than others. d is fixed to be −1 since it has little influence on χ 2 .
For δ sys = 0.002, the best fitted values of Λ QCD are 0.1086 (17), 0.1350 (21), 0.093(10), 0.086 (14), 0.0926 (61) GeV for the correlations in the overlap quark, clover quark, overlap pion, clover pion, and clover nucleon, respectively. As one can see, they are quite similar except for the correlation in the clover quark state. While the size of the systematic error has some effects for the correlation in the quark case, it has little influence on the correlation in the physical states.
In subfigure (b) in each figure, after taking Λ QCD to be the best fitted value, we use Eq. (19) to fit the bare matrix elements to extract g(z). The range of z is taken from 0.06 fm to 1.02 fm with 17 different z values for correlations in the overlap quark case, from 0.06 fm to 1.14 fm, with 19 different z values for correlations in the overlap pion and clover pion cases, from 0.06 fm to 0.96 fm, with 16 different z in the clover quark case, from 0.06 fm to 0.9 fm, and with 15 different z values in the clover nucleon case. We fine-tune d to make sure that In subfigure (d) in each figure, we show our result for the renormalized matrix element Exp[g(z) − m 0 z] in blue points connected by lines. As z increases, Exp[g(z)−m 0 z] increases first, following the perturbative result, and then decreases from about z = 0.25 fm on. In all the cases, Exp[g(z) − m 0 z] is consistent with Z MS at small z. However, at large z, there is a significant discrepancy between Exp[g(z) − m 0 z] and Z MS . This means that there is a large non-perturbative effect at large z in the popular RI/MOM and ratio renormalization scheme used previously, which supports the hybrid renormalization procedure proposed recently [39]. We have already briefly discussed this in Sec. V C.
Subfigure                (z, a) itself. This subfigure is a consistency check of our method. In all cases, M(z, a)/Z(z, a) R has little dependence on a within error and is consistent with Exp[g(z) − m 0 z]. Therefore, our self-renormalization method can isolate the residual intrinsic physics from the linear divergence, renormalon uncertainty, log divergence and discretization error with a reasonable precision.
Although the ratio M(z, a)/Z(z, a) R has a good behavior for most of the lattice spacings, it may have large errors for small lattice spacings like 0.0318 fm or 0.0574 fm. Therefore we do not recommend using M(z, a)/Z(z, a) R as the renormalized matrix element. Since Exp[g(z) − m 0 z] is achieved through our fitting, the large-error data points have little influence on it. We take this as the renormalized matrix element, while using M(z, a)/Z(z, a) R only to test consistency.
We show in Fig. 17 a comparison of renormalized O γt (z) correlation in the pion state between overlap and clover actions. Before we add the systematic error, both results show appreciable difference, which may indicate that two different valence fermions will lead to different results. After we add a small dummy systematic error during the fitting, the correlations in both cases lead to similar results. That means that the systematic error is important, and the residual intrinsic non-perturbative physics is independent of valence quark formulations if it is properly included.
Finally, for curiosity, we show in Fig. 18 the renormalized correlations as well as parameters related to the divergence for all cases we studied. As mentioned before, the fitted Λ QCD for overlap quark, overlap pion, clover pion and clover nucleon are similar to each other. But for the clover quark, it is markedly different. We do not yet understand the reason, but it could be due to the combination of chiral symmetry breaking and the linearly-divergent quark mass. We plan to investigate this in the future, but for the time being, at least we understand why we cannot eliminate the linear divergence using RI/MOM for the clover action [45]. Quite surprisingly, the residual intrinsic non-perturbative correlations for overlap quark, pion and nucleon are very similar to each other, which we also do not fully understand.

VI. FITS WITH OTHER OPTIONS AND STABILITY
In a phenomenological analysis, the amount of information and the accuracy one can obtain depends, obviously, on the quality of the lattice data. Given the data we have, we would like to push the limit of the analysis by including more physics in the fit until the analysis no longer yields useful information.
The first question we would like to address is whether one can get more accurate information about the linear divergence. We have used so far the one-loop perturbative QCD prediction, but treated Λ QCD as a free parameter to partially take into account higher order corrections. However, one might consider directly including at least the second-order corrections. Our strategy here is to fix k to its perturbative value, and to vary Λ QCD and the explicit second order correction characterized by a new parameter λ where the terms omitted are the same as those we had before in Eq. (19). k is related to k by −2π/b 0 (since −k 2π/b 0 = k). On the other hand, one might use the λ parameter to include some higher order corrections [51], ln M(z, a) = k z aα To be consistent, we now have to use the QCD running coupling constantα s up to the second order β function b 1 = 102 − 38 3 n f , Therefore, again we have two global parameters λ and Λ QCD in the fit.
The fitting process with the above parametrization is similar to what is described in Sec. V D. The results are shown in Figs. 19 and 20. As one can see, the range of λ and Λ QCD are rather large, and strongly correlated. After picking several correlated values of the pair, matching to the perturbative result, the resulting residual matrix elements are shown in the lower panel. One can hardly see any difference between the different choices. Moreover, the results are essentially the same as the fit without the extra λ parameter (Fig. 22). Therefore, we conclude that with the data at hand, our result is stable with respect to the higher order corrections. Only with more accurate data, one might be able to see the difference of a higher-order fit.
Another possibility one can explore is that the higherorder corrections represented by Λ QCD might be different for different ensembles. Therefore, one can discuss different Λ QCD values for MILC and RBC data, where we fix the k parameter and obtain a two parameter fit, Λ QCD1 and Λ QCD2 . The results of using the above equation is shown in Fig. 21. As shown, the two parameters are strongly correlated and proportional to each other. The difference of the two is limited to a very small range. Furthermore, with different choices of the correlated pair (Λ QCD1 , Λ QCD2 ), the final renormalized results do not change.
In Fig. 22, we show the results of renormalized matrix elements based on different fitting functions. All of these fitting functions lead to the same result. Thus, given the data, it is unnecessary to consider other high-order terms in the fitting function and Eq. (19) is good enough. This also supports the argument made in Sec. V B that we can use the fitted Λ QCD in leading order to partially represent higher order corrections.
Finally, we consider O(a 2 ) corrections to our fitting formula. These fits were unstable, which indicates that we introduce too many fit parameters for the given data quality. The conclusion is that higher precision data is needed to constrain O(a 2 ) terms. One can also make fits by including O(a 2 ) instead of the linear correction. The physics for doing it is weak in our cases.

VII. RENORMALIZATION OF LATTICE SMEARED MATRIX ELEMENTS
In many lattice calculations, the data can be very noisy which calls for more statistics. However, large statistics means more resources which are hard to come by sometimes. Therefore, lattice practitioners have invented phenomenological approaches to quench the short-distance fluctuations: smearing. However, smearing might lead to configurations that are not connected with fundamental theory. On the other hand, smearing in some sense just makes the effective UV cutoff (the inverse of the lattice spacing) smaller, and in the limit a → 0, all smearing becomes local effects which can be taken into account properly through some kind of renormalization. In this section, we check if the procedure we discussed in the previous section still works for smeared matrix elements.
To ensure that the HYP smearing data are useful for refined analysis, we can do a simple test on whether the HYP smearing data show the properties of the linear divergence predicted by perturbation theory, as we have done in Sec. IV. Here we use the following function to fit the bare matrix element M to extract the linear divergence factor, ln M(z, a) = e(z) a ln[aΛ QCD ] + g(z) + f 1,2 (z)a, where we allow for a discretization error to get a better fit. We fix Λ QCD at 0.39 GeV in this simple test. Fig. 23 shows the fits for different actions and states which all look reasonable. We have not shown χ 2 which will be considered in the extraction of the renormalization factor. Fig. 24 shows the extracted e(z) with respect to z. The linear z dependence is approximately reproduced but some deviations are seen at large z. We can perform a linear fitting of e(z) and compare the fitted slopes. Overall, the quality of the fits is not as good as in the unsmeared cases, although it is still rather impressive. The slope extracted for the linear divergence for the clover quark is totally different from the others, which explains why we cannot use the RI/MOM factor to eliminate it in the clover case for HYP smearing data [45]. Although the slopes for the overlap quark and pion, clover pion and nucleon are similar to each other, there are some small discrepancies. Using the RI/MOM factor for these cases may leave small residual linear divergences, and the selfrenormalization method provides a better way to renormalize the matrix elements. Next, we extend our self-renormalization method to the HYP smeared data. Our fitting function, Eq. (19) may have no clear physical meanings for the HYP smeared date, but it may be regarded as a phenomenological model.
We start with a fit function with two parameters in the linearly divergent part, k and Λ QCD . Fig. 25 shows the χ 2 map with respect to them. It is clear from the plots that in the "small-χ 2 band", we cannot find a solution for k with the one-loop perturbative value 7.4 GeV −1 fm −1 (1.46 if dimensionless, see Eq. (26)) any more. This is expected because the lattice space is enlarged by a factor of two, and therefore the values of k in the "small-χ 2 band" are expected to be at most half of this, named 3.7 GeV −1 fm −1 (0.73 if dimensionless). In fact, most probable k is now around 3.0 GeV −1 fm −1 (0.59 if dimensionless) or below. Another way of saying this is HYP smearing smooths out the short-range behaviors such as the linear divergence. The value of Λ QCD , on the other hand, has a larger range, possibly corresponding again to a larger effective lattice spacing.
So for the HYP smearing cases, we choose a set of (k, Λ QCD ) near the center of the "small-χ 2 band". We choose k = 2.5 GeV −1 fm −1 (0.49 if dimensionless) for overlap quark and clover quark, and k = 2.6 GeV −1 fm −1 (0.51 if dimensionless) for overlap pion, clover pion and clover nucleon. The corresponding Λ QCD shows a larger dispersion, with Λ QCD = 0.43 and 0.38 GeV for overlap quark and pion correlations respectively, and Λ QCD = 0.61, 0.37 and 0.35 GeV for clover quark, pion and nucleon correlations, respectively.
Once the (k, Λ QCD ) parameters are chosen, we subtract the linear divergences, and match the result to the perturbative matrix elements to determine possible further subtractions for the non-perturbative mass effect. Fig. 26 shows that the renormalized matrix elements for HYP smearing in blue ticks and lines, in comparison with the cases without HYP smearing.
The large difference between smeared and unsmeared matrix elements is seen in the case of clover quark correlation at large z. Again, we speculate that due to the chiral symmetry breaking effects, the linear divergence has some curious behavior there. However, the difference is smaller in the pion matrix element where the error bars are also bigger. For the overlap case, the correlations for quark and pion show no substantial change due to smearing effects. On the other hand, the matrix element in the nucleon case shows a much smaller error bar in the smeared case, which demonstrates the power of smearing.
To conclude, our analysis method to renormalize the linear divergent matrix elements can successfully be used for smeared matrix elements as well. Smearing does not seem to change the behavior of renormalization qualitatively, but only quantitatively. The intrinsic physics does not seem to change under (moderate) smearing, consistent with the expectation that smearing is an alternative method to reduce the statistical noise in lattice calculations.

VIII. LINEAR DIVERGENCE OF CORRELATIONS IN QCD VACUUM
In this section, we analyze the linear divergence of the Wilson line in the matrix elements taken in the QCD vacuum. We try to test if the divergence is universal and if the divergent mass can be extracted successfully from these matrix elements. We consider three new types of matrix elements: 1) large-size Wilson loop that has been used to extract heavy-quark potential, 2) the vacuum matrix element of the quasi-PDF operator, 3) the vacuum matrix element of Wilson link operator in a fixed gauge.
The data are for MILC ensembles and, when applicable, with the clover action. According to the last section, smearing does not change the form of renormalization, and therefore, we consider here the higher-precision data with one step smearing. Since we just want to test the linear divergence, we can use the simple function Eq. (32) in our fitting.

A. Wilson Loop
To extract the linear divergence from lattice data, a Wilson loop is a good choice because it is easy to calcu-    late, gauge invariant, and has high precision. It has been extensively studied in the literature [27,28,38,43,59], particularly in relation to the heavy quark potential.
In our analysis, we will not consider the heavy quark potential interpretation but simply consider a Wilson loop as a matrix element and fit its 1/a divergence by looking at the logarithm of the matrix element, as shown in Fig. 27. The lattice spacing dependence has been fitted well with our formula with a choice of Λ QCD = 0.39 GeV, which is a value favored by our fitting in the previous section. After isolating the linear divergent coefficient, we plot it as a function of r + t, for which half of the slope gives us the divergent coefficient. Our value is about 2.5 GeV −1 fm −1 (0.49 if dimensionless). This value is very similar to the value we found in the previous section (Fig. 24), demonstrating that the linear divergence can be very well extracted from the Wilson loop.

B. Vacuum Matrix Element of Quasi-PDF Operator
Here we consider the vacuum matrix elements of a quasi-PDF operator which were suggested as choices for the renormalization factor in Ref. [44].   One can consider the vacuum expectation value (VEV) of O γt (t). As a gauge invariant choice proposed in Ref. [44]. O Γ (t) can be obtained through the following stochastic estimation: where the second term on the right hand side of the second line vanishes as it is not gauge invariant, and S w ( x, t) = y S( x, t; y, 0) is a wall source quark propagator without gauge fixing. Such a proposal can only be applied for the non-vanishing cases with Γ = γ t or I. We will apply it to O γt . However, it is simple to see that the O γt matrix element in the vacuum vanishes at small t using operator product expansion. Therefore the matrix element is susceptible to large O(a) correction. Only at very large t, one might see the correct slope.  Results of a test of the linear divergence for the vacuum matrix element of the PDF operator is shown in Fig. 28. The slope at large t appears to be consistent with that from the Wilson loop case but only within uninterestingly large errors.

C. Landau-gauge-fixed Wilson link
For the gauge-dependent matrix elements, we consider the Wilson line in Landau gauge as the simplest choice.
The result of our test of the linear divergence is shown in Fig. 29. At small-z, where perturbation theory works well, the slope roughly agrees with the prediction from perturbation theory. For large-z, the matrix element does not have a transfer-matrix interpretation in this gauge.
Finally, Fig. 30  others because the large discretization errors make the extraction of the linear divergence very imprecise.

IX. SUMMARY
The proper renormalization of quasi-LF correlations calculated on the lattice has been a main challenge for the applications of LaMET to studying parton physics. Such correlations contain linear divergences associated with Wilson lines which have to be eliminated with high precision to extract the desired physics information. In this work, we propose a self-renormalization method to eliminate both the linear and the logarithmic divergence, as well as the discretization error, from a quasi-LF matrix element which can be matched to a continuum scheme at short distance. As a paradigmatic example, we show that our method works well for O γt matrix elements in the pion, nucleon, and off-shell quark state. Our analysis shows that the renormalization factors are universal in the hadron state considered. Moreover, the renormalized correlations in the pion state for clover and overlap fermions are similar to each other. Besides, We find a large non-perturbative effect in the popular RI/MOM and ratio renormalization scheme used previously, which has to and can be avoided in the hybrid renormalization scheme proposed recently. It is certainly interesting to see that our method works also for HYP smeared matrix elements. This paves the way to do a phenomenological renormalization for smeared high-precision lattice data.
For the convenience of the reader, we collect the parameters related to renormalization for various matrix elements discussed in the previous sections, see Tables III  and IV. For HYP smearing cases, there is much freedom for the choice of Λ QCD . We choose Λ QCD = 0.39 GeV in Table IV. But we also find that Λ QCD = 0.1 GeV gives us smaller fitted discretization errors.  For comparison purposes, we also show results based on 2+1 flavor clover quark and Luescher-Weisz (equivalent to Symanzik) gauge ensembles from CLS collaboration [71]. We choose valence fermion action to be clover action, the same as sea fermion action. The renormalization parameters for CLS ensembles are shown in Table V. The parameters k and Λ QCD for pion matrix elements for CLS ensembles are similar to those for MILC and RBC ensembles, suggesting consistency between results using mixed actions or not. We would like to stress that the application of our method is not limited to the matrix elements discussed in this paper. Instead, it can be applied to any matrix element of an operator in the form of Eq. (2). It can also be generalized, in principle, to any matrix element with a Wilson line in LaMET applications [8]. Of course, in these cases one also needs to consider the physical interpretations of the matrix elements and the details related to lattice calculations carefully.