Improved renormalization scheme for nonlocal operators

In this paper we present an improved RI-type prescription appropriate for the non-perturbative renormalization of gauge invariant nonlocal operators. In this prescription, the non-perturbative vertex function is improved by subtracting unwanted finite lattice spacing ($a$) effects, calculated in lattice perturbation theory. The method is versatile and can be applied to a wide range of fermion and gluon actions, as well as types of nonlocal operators. The presence of operator mixing can also be accommodated. Compared to the standard RI' prescription, this variant can be recast as a supplementary finite renormalization, whose coefficients bring about corrections of higher order in $a$; consequently, it coincides with standard RI' as $a\to 0$, however it can afford us a smoother and more controlled extrapolation to the continuum limit. In this proof-of-concept calculation we focus on nonlocal fermion bilinear operators containing a straight Wilson line. In the numerical implementation we use Wilson/clover fermions and Iwasaki improved gluons. The finite-$a$ terms were calculated to one-loop level in lattice perturbation theory, and to all orders in $a$, using the same action as the non-perturbative vertex functions. We find that the method leads to significant improvement in the perturbative region indicated by small and intermediate values of the length of the Wilson line. This results in a robust extraction of the renormalization functions in that region. We have also applied the above method to operators with stout-smeared links. We show how to perform the perturbative correction for any number of smearing iterations, and evaluate its effect on the power divergent renormalization coefficients.


I. INTRODUCTION
The past decade has seen a great surge in the study of parton distribution functions (PDFs), by means of calculations formulated on a spacetime lattice. The viability of a Euclidean formulation for these intrinsically Minkowski-space functions was argued in ground-breaking papers by Ji [1,2] with the introduction of the Large Momentum Field Theory (LaMET) and the calculation of the so-called quasi-distributions. A number of approaches are currently being employed in addition to quasi-distributions in order to match experimentally measurable distribution functions to quantities amenable to numerical estimation on a lattice, such as the pseudo-distributions [3,4]. Both of these cases rely on the evaluation, via lattice simulations, of matrix elements which involve gauge-invariant nonlocal composite operators. There is a large variety of nonlocal operators which are of physical relevance, notably path-ordered exponentials of the gluon field along a straight-line or staple-shaped contour ("Wilson line"), with a quark-antiquark pair or two gluon field tensors attached at the endpoints. For an extended list of general references, we refer the reader to Refs [5][6][7][8][9].
The study of nonlocal operators in Minkowski continuum spacetime has a long history, dating back decades before lattice investigations began, with seminal works by Mandelstam, Polyakov, Makeenko and Migdal [10][11][12], and later works, among them Refs. [13][14][15][16][17][18][19]. The lattice formulation introduces several new complications, such as mixing and renormalization functions which exhibit an inverse-power divergence as the regulator a (lattice spacing) approaches its limiting value, a → 0. Power divergences can of course manifest themselves also in the renormalization of local operators, but only in the presence of mixing with lower-dimensional operators which share the same symmetry properties; nonlocal operators, on the other hand, by virtue of the natural length scale(s) appearing in their definition, exhibit power divergences despite the absence of mixing with lower-dimensional operators. Further, the effect of finitea corrections -lattice artifacts -becomes ever more significant as the spatial extent of the operator grows.
A number of approaches are routinely applied to alleviate finite-a artifacts. Many of them apply improved discretizations to terms in the action and to the definitions of composite operators; such improvements amount to adding contributions of higher order in a, which is a legitimate and systematic procedure in the spirit of the Symanzikimprovement program ( [20,21], for one of the first lattice applications see, e.g., [22]). Exploitation of symmetries can also lead, in some cases, to elimination of O(a 1 ) (or even O(a 2 ) ) effects (e.g., [23], and [24,25]). A recent approach developed for nonlocal operators, the hybrid scheme [26], implements a different renormalization scheme for small and large values of the length of the Wilson line, z. At short distances one may use an RI type or a ratio scheme. At large values of z one can use the Wilson-line-mass-subtraction scheme [27,28].
Another approach, which will be employed in the present work, has been quite successful in the renormalization of various local operators on the lattice; it relies on perturbation theory carried out to all orders in a (see, e.g., Refs. [29][30][31][32][33] for comparable implementations). In this approach, the Green's functions which are to be evaluated non-perturbatively, in order to stipulate operator renormalization conditions, are also computed perturbatively, keeping the complete dependence on a at the i-th perturbative order 1 , O(g 2i ). The perturbative result, generically denoted as Λ (i) (a), can be split into contributions which are O(a 0 ), O(a 0 ln(a)) [possibly even O(a −1 )], and remaining contributions which vanish as a → 0. This splitting is most straightforwardly performed to one loop: rest (a) ; the quantity Λ (1) 0 (a) could be used to derive standard perturbative renormalization functions, should this be called for, while Λ (1) rest (a), which consists purely of lattice artifacts, can be subtracted from the non-perturbative value of the corresponding Green's function. The latter thus becomes a smoother function of a, and allows for a more stable extrapolation to the continuum limit. The legitimacy of this procedure can be best exposed by recasting the above subtraction in the form of a complementary renormalization containing only higher orders in a [i.e., supplemental multiplicative renormalization functions of the form: Z(a) = 1 + g 2 z(a), (z(0) = 0) and mixing of the original operator O with higher-dimensional operators O i , thus leading to a renormalized operator Ideally, the elimination of artifacts would circumvent the need to carry out an extrapolation to vanishing lattice spacing; this, unfortunately, is not practicable since, in the absence of such an extrapolation, terms of O(g 2i a n ), (i > 1, n > 0), which are beyond computational feasibility, would still be significant for typical values of a used in simulations.
In this work, we apply this last approach to nonlocal operators; we focus on straight Wilson lines ending on a quark and an antiquark field. It is worth noting that, in handling nonlocal operators, elimination of lattice artifacts is performed separately for different values of the operators' length scale(s), since each value corresponds to an independent operator. One of the advantages of this method is that it can be integrated in the hybrid scheme to improve the estimates for the small-z region.
Another aspect of nonlocal operator improvement which is addressed in this work regards stout smearing [34] and, in particular, the elimination of lattice artifacts from Green's functions of operators containing stout links. Implementation of stout smearing is straightforward and can be iterated an arbitrary number of times ("smearing steps") in the context of a numerical simulation. On the contrary, perturbative calculations in the presence of stout links are notoriously cumbersome, and can rarely be carried out beyond one step [35]. We demonstrate that, for the operators under study, as well as for a rather wide class of local and nonlocal operators, one-loop perturbative expressions can be obtained as closed form integrands over the loop momentum, for an arbitrary number of smearing steps. This allows a perturbative treatment of artifacts, using the procedure outlined above.
The layout of this paper is as follows: Section II provides the general setup and definitions of the objects entering the calculation. In Section III we describe the renormalization scheme and the improvement procedure. Sections IV and V contain our perturbative results to all orders in a, and their effect of non-perturbative renormalization functions, respectively. Section VI discusses stout smearing and the results of its application to one loop, for an arbitrary number of smearing steps. Finally, Section VII summarizes our findings.

A. Lattice Actions
We calculate the Green's functions of nonlocal operators using the clover (Sheikholeslami-Wohlert) fermion action [36] in which the clover parameter, c SW , is treated as a free parameter. The Wilson parameter, r, is set to 1; f is a flavor index, σ µν = [γ µ , γ ν ]/2 andF µν /(ia 2 g) is a clover discretization of the gluon field tensor, defined aŝ with Q µν being 2 The Lagrangian masses for each flavor, m f 0 , are set to their critical value, which to one-loop level is zero.
We present results for the Symanzik-improved gluon action [37] with c 0 = 3.648, c 1 = −0.331, and c 2 defined by the normalization condition c 0 + 8c 1 + 16c 2 = 1. We will focus on the Iwasaki action, for which c 2 = 0. Our calculation is readily applicable for a family of Symanzik improved gluon actions, as defined in Ref. [37].

B. Definition of Operators
We study fermion bilinear nonlocal operators of the form where the path-ordered (P) Wilson line is inserted for gauge invariance. We consider only cases where the Wilson line is straight, and, without loss of generality, we choose it to be along theẑ direction. The length of the Wilson line is z, and our perturbative calculation implements values of z which range between z = 0 and z = 15. This range is sufficient for the lattice sizes used currently in non-perturbative calculations of the renormalization functions. Note that, in the limit z → 0, Eq. (5) reduces to the standard ultra-local fermion bilinear operators. However, the calculation of the Green's functions for O Γ is for strictly z = 0 as the limit z → 0 is nonanalytic.
We perform our calculation for all 16 independent combinations of Dirac matrices, Γ, and we distinguish between the cases in which a Lorentz index of a Dirac matrix is in the direction z of the Wilson line or not. Therefore, the possible choices of Γ are separated into 8 subgroups, defined as where j, k ∈ {t, x, y}. An alternative definition for the tensor operators is O γ 5 σ µν . This is redundant if one employs a 4-dimensional regularization, such as the lattice, since the latter operators are just a renaming of the T operators, and they will, thus, renormalize identically.
Ref. [38] revealed a pairwise finite mixing pattern in lattice regularization for some of the operators of Eq. (6), based on one-loop perturbative calculations in lattice QCD. The mixing pattern is in the pairs {S, V z }, {A t , T xy }, {A x , T ty }, {A y , T tx }, and this holds to all orders in perturbation theory. The improvement method that we propose here can be applied to the cases with or without mixing.

C. Perturbative Calculation
The central focus of this work is the evaluation of the lattice-regularized bare amputated Green's functions to one-loop level in perturbation theory and all orders in the lattice spacing, O(g 2 a ∞ ). Exploiting translational invariance in this equation, the operator O Γ is summed over all space-time positions, and thus the external quark and antiquark Fourier transformed fields have the same momentum (p). Λ 1−loop O explicitly depends on p and on the length of the Wilson line z. Such a calculation requires more sophisticated methods of integration over the loop momentum, but is, overall, less complicated than the calculation to O(g 2 a 0 ), which must also be carried out. Arriving at analytic expressions for the latter necessitates that the divergent terms be isolated; this is a rather delicate procedure involving special techniques [38]. On the contrary, the Green's functions O(g 2 a ∞ ) are evaluated numerically for each value of the external momentum that is used in the numerical calculation of the renormalization functions; this momentum serves as the renormalization scale. Besides these 1-loop calculations, we evaluate non-perturbatively Λ 1−loop O , traced with its tree-level value. Such a projection, described in more detail in Sec. III, eliminates the Dirac structure, leading to a computationally less costly calculation. Other kinds of projectors can also be employed.
The one-loop Feynman diagrams which enter our calculation are shown in Fig. 1. The filled rectangle shows the insertion of any one of the nonlocal operators O Γ . Diagram d1 contains the 0-gluon vertex of the operator, whereas diagrams d2-d3 (d4) contain the corresponding 1-gluon (2-gluon) vertex. The "tadpole" diagram d4 is related to the self-energy of the operator and leads to the well-known linear divergence with respect to the dimensionful ultraviolet cutoff a [15,39]. As noted in Ref. [38], such a divergence is independent of the choice of the operator, and, to O(g 2 a 0 ), has the form c |z|/a. The coefficient c depends, however, on the gluon action.

A. Non-perturbative calculation
For the renormalization condition we employ a regularization independent RI scheme [40], which requires the calculation of the non-perturbative vertex functions, G O (p, z). The latter correspond to a nonlocal operator within The vertex functions with momentum p are amputated using the up-quark and down-quark propagator in momentum space, that is, . This amputated vertex function is matched to its tree-level value, Λ tree O (p, z), using the RI prescription, where the vertex momentum is set equal to the renormalization scale. We define the renormalization functions for the quark field and quark operator as where the quark field renormalization, Z q , is given by matching the propagator to its tree level value Λ tree q (p) Eq. (9) is a generalization of the condition used for local operators at a scale µ 0 ; here it is applied at each value of z separately.

B. Improvement scheme
The improvement procedure is applied on Eq. (9) using the perturbative expressions for the O(g 2 a ∞ ) results obtained through the diagrams of Fig. 1. A similar improvement is applied on Z q following the work of Ref. [33]. More precisely, we redefine the renormalization functions as follows: where ∆Λ is the difference between the Green's functions to O(g 2 a ∞ ) and the Green's functions to O(g 2 a 0 ), defined above as Λ rest (a). With this improvement scheme, one subtracts the unwanted finite-a contamination from the nonperturbative trace of the vertex functions using the results from 1-loop perturbation theory. There are alternative subtraction schemes, such as to subtract the finite-a contamination on Z O . The latter differs from the one in Eq. (11) to O(g 4 a n ) (n ≥ 1). As described in the Introduction, the variations of improvement schemes can also be understood as a finite renormalization whose effects disappear in the limit a → 0.
Once the improvement procedure is applied on the renormalization conditions, the analysis proceeds as outlined in Ref. [41], that is 1. Calculation of the vertex functions using isotropic ("democratic") momenta in the spatial directions, ap = 2π L ( nt 2 + 1 4 , n x , n x , n x ). We employ antiperiodic boundary conditions in the time direction. In addition, the momenta should have a ratio as close to 0.25 as possible (ideally below 0.3 for nonlocal operators). This is based on empirical arguments, since the Lorentz non-invariant function P 4 appears in O(g 2 a 2 ) contributions in lattice perturbation theory.
2. Implementation of an improvement scheme that eliminates finite-a contamination, such as Eq. (11).
3. Chiral extrapolation using the estimates of Z O obtained from ensembles with different values for the quark masses.

4.
Conversion of the chirally extrapolated Z O to the modified MS scheme (MMS) scheme [41] and evolution to a common renormalization scale. The appropriate expressions for nonlocal operators can be found in Ref. [38].

5.
Elimination of residual dependence on the initial renormalization scale using a linear fit in (aµ 0 ) 2 .

IV. PERTURBATIVE RESULTS
It is instructive to study the finite-a contributions extracted from our perturbative calculation, in order to assess the magnitude of the artifact contamination in the renormalization functions. For demonstration purposes we show the artifacts in the Landau gauge for c SW = 1, and lattice size L = 24. For the bare coupling constant (g 2 0 = 6/β) we choose β = 2.10. These parameters match those of the calculation presented in Ref. [41]. We note that the choice of the optimal tree-level value c SW = 1 is consistent with 1-loop perturbation theory, and not with the value used in non-perturbative calculations. Using c SW = 1 in perturbative calculations is the standard treatment of the clover parameter.
The one-loop lattice artifacts are embodied in the following expression: where The factor of 1 4 is the normalization assuming that the trace is over the Dirac indices. We note that T O (z, p) is a complex function due to the presence of the Wilson line. Numerically, we find that T O is similar for all operators, and has mild c SW dependence. Therefore, we focus on the vector operator, Γ = γ 0 , for c SW = 1.   [33]. Similarly, Fig. 3 corresponds to the imaginary part of T V for z = 1, 3, 6, 9, 12, 15. Given that T O is purely real for z = 0, we include z = 1 instead. For both parts, we observe that the seven classes of momenta in the spatial direction (n x ∈ [0, 6]) contributing to the same z have a different behavior compared to each other that exhibits discontinuities. This holds for any value of z, including z = 0. Such a "wiggling" effect has also been observed in the non-perturbative renormalization of local and nonlocal operators [33,41]. On the contrary, the members of each class of n x (n t ∈ [0, 10]) follow a smooth functional form, a feature also observed in non-perturbative calculations.
There are a number of observations regarding the behavior of the finite-a terms. We find that Re[T V ] is negative for z > 0 for all momentum classes except the ones with spatial components in the class (n t , 0, 0, 0). The artifacts for both Re[T V ](0) and Im[T V ] have negligible values for the class (n t , 0, 0, 0). Another observation is that for n

V. NON-PERTURBATIVE IMPROVED RESULTS
Here we demonstrate the effectiveness of the proposed renormalization prescription using the vertex functions of the vector operator, Γ = γ 0 , which is free of mixing and is used for the calculation of physical matrix elements related to the unpolarized PDF. We checked that the conclusions drawn here are also valid for the axial and tensor operators. We use two ensembles of N f = 2 twisted mass fermions with a clover term, and Iwasaki gluons. The volume is 24 3 ×48, lattice spacing a = 0.0938 fm, and the corresponding values of the pion mass are 235 MeV and 340 MeV. More details can be found in Refs. [33,41]. We find that the level of improvement is the same for both ensembles, and present results for the 235 MeV case. It should be noted that the non-perturbative vertex functions have been obtained using five steps of stout smearing, while the perturbative calculation is performed without a smearing. Therefore, the effect of the subtraction of lattice artifacts can potentially be further improved when the stout dependence is calculated. We refer the reader to Sec. VI for the calculational setup.
The effect of the improvement scheme is best evaluated after the Z O estimates are converted to the MS or MMS schemes, where the residual dependence on the initial RI scale is O(g 4 ). Here, we use the MMS scheme [41], which is employed by the Extended Twisted Mass Collaboration (ETMC). All data have been evolved to a scale of 2 GeV using the conversion factor of Ref. [38].
As mentioned previously, the non-perturbative renormalization functions should ideally be determined on spatially isotropic momenta with P 4 as close to 0.25, which correspond to the fully isotropic case (temporal and spatial). Such choices put very strict limits on the number of RI' scales in a range in which both hadronic contaminations are small and perturbation theory is valid. It is, thus, desirable to expand the range of momenta by relaxing the criterion for P 4. Here, we use momenta that satisfy 0.25 ≤ P 4 ≤ 0.41, and have (ap) 2 up to 7. Fig. 5 shows the raw non-perturbative data for Z V (1), Z V (2), Z V (3), as well as the corresponding improved ones. Similarly, Figs. 6 -7 show Z V for z/a = 4 − 6 and z/a = 7 − 9, respectively. Focusing on the real part of Z V , we find that the effect of the subtraction of lattice artifacts as proposed in Eq. (11) works very well up to z = 7a ∼ 0.65 fm. Beyond that region, the subtraction method does not perform well, particularly for P 4 > 0.29. This is not surprising because such values of z exceed the perturbative region, and the one-loop perturbative results become unreliable. Nevertheless, the method appears to be successful for z < 0.65 fm. In particular, the "wiggling" effect between the different classes of momenta disappears for almost all values of P 4. This is a highly nontrivial feature of the improved estimates. In addition, we observe the formation of good plateaux that allows one to reliably take the limit (aµ 0 ) 2 → 0. The quality of the plateaux improves further if we apply a cut on the data at P 4 = 0.3, as shown in the figures. It is noteworthy that the raw non-perturbative data produce a linear dependence in (aµ 0 ) 2 for P 4 = 0.25, which, however, has a sizeable slope leading to an unreliable estimate at (aµ 0 ) 2 → 0.
We also find interesting features in the imaginary part of Z V . The first observation is that for z ≤ 4a, the subtraction of lattice artifacts results in an estimate which is closer to zero than the purely non-perturbative ones. In fact, the imaginary part of the improved estimates is always closer to zero in the region (aµ 0 ) 2 < 5. This is a desired feature because in dimensional regularization the poles are real to all orders in perturbation theory. The same feature should hold for non-perturbative calculations if a conversion factor to high enough order in g 2 were available. In addition, the improved estimates (Eq. (11)) for the small-z region (z < 0.4 fm) have a smoother dependence on (aµ 0 ) 2 compared to the purely non-perturbative estimates (Eq. (9)). As z increases, discontinuities are observed between different classes of momenta, in particular for (aµ 0 ) 2 3.5. A milder effect is observed for momenta satisfying P 4 < 0.3. It is interesting to observe that the purely non-perturbative data have a smoother (aµ 0 ) 2 dependence in the large-z region compared to the improved ones. This indicates that the real and imaginary parts have different contamination from finite-a contributions. Thus, one can use an improvement scheme in which the subtraction of the artifacts can be applied differently in the real and imaginary parts. Also, including stout smearing in the perturbative calculation as described in Sec. VI, has the potential to further improve the effect of the subtraction.

VI. STOUT IMPROVEMENT
A standard way of improving the behavior of operator matrix elements is implemented by applying stout smearing to link variables [34]. This procedure is typically performed as an iteration of smearing steps; starting from the original link variables U µ (x + a µ/2)), the i th iteration step performs the following replacement [42]: where the definition of Q µ (x) is [the superscript (i), denoting smearing step, is implicit]: The "stout coefficient" ω is a tunable parameter and V µ (x) represents the sum over all staples associated with the link U µ (x).
The implementation of stout smearing on links which are produced in numerical simulations is straightforward; several smearing steps are generally applied in order to optimize improvement in measured quantities. On the other hand, calculations in lattice perturbation theory, in the presence of a smeared action or of smeared operators, are prohibitively complicated even at one loop, given the extraordinary proliferation of terms in the corresponding vertices; as a result, such calculations can usually be performed only when a maximum of two smearing steps have been applied. [35,43].
Despite the sheer complexity of constructing and using vertices with stout links, there are two general classes of vertices which can be expressed in a compact form, even for an arbitrary number of smearing steps: The first class regards vertices in which only one field is a gluon. The second class regards vertices in which two fields are gluons; this class leads to compact expressions only for that part of the vertices which does not vanish when the color indices of the two gluons are contracted (i.e. the part which is free of commutators). Thus, there are several instances of Green's functions (GFs), typically at one loop, which can be thoroughly computed with relative ease for any number of smearing steps. Such instances include: (i) the fermion propagator, (ii) GFs with external fermion-antifermion lines and an insertion of any fermion local bilinear operator of arbitrary twist, (iii) GFs with external fermion-antifermion lines and an insertion of a fermion nonlocal bilinear. The latter instance is of course the focus of the present work; however, the explanation which follows regards generic one-and two-gluon vertices.
Let us first derive the terms up to O(g) (one-gluon) in a link after one smearing step. In terms of the Fourier transformed gluon fieldÃ (0) µ (p), one obtains: [For conciseness, an integral π −π d 4 p e i p·x /(2π) 4 , accompanying each gluon fieldÃ ρ (p), is left implicit.
x ≡ x/a is a four-vector with integer components.] From the above expression it follows that the longitudinal part of the gluon field remains intact, while the transverse part is multiplied by a factor (1−ωp 2 ). This feature propagates to successive smearing steps, with potentially different stout coefficients (ω, ω ): Thus, a succession of N steps with the same coefficient ω leads to: Terms with two or more gluons in the exponent of smeared links (i.e., the O(g 2 ) terms inÃ (N ) µ (p)) are considerably more complicated. However, given that U (N ) µ (x) is special unitary by construction, two-gluon terms in its exponent will necessarily involve a commutator; consequently, such terms will give a vanishing contribution in Feynman diagrams where the color indices of the two gluons end up being identified. In such diagrams, the two-gluon expression arising from a smeared link is simply: Let us apply the above to the zero-, one-and two-gluon vertices of the nonlocal operator O (N ) Γ , defined as in Eq. (5), with N -fold smeared links. Setting z = a n, n ∈ N, we obtain: where, omitting commutators : The appearance of the stout parameter ω exclusively in the combination (1 − ωp 2 ) N provides a point of reference for the numerical values of ω to be employed: Given that the maximum value ofp 2 is 16, and its average value is 8, one may expect that values of ω in the range 1/16 ω 1/8 tend to eliminate the transverse part of the gluon fields after several smearing steps. This, in particular, implies a vanishing result for diagrams d2, d3, d4 (see Fig. 1) in the Landau gauge, leaving only diagram d1; since the latter coincides with its continuum counterpart, it follows that there would be no lattice corrections to this order. Conversely, values of ω beyond ω 1/8 risk a disproportionate increase of lattice artifacts in the transverse gluon field. Non-perturbatively, the value of ω can be tuned using the criterion that the plaquette reaches maximum value for a given number of smearing steps, and typical values for ω are around 0.1.
As a demonstration, we show numerical results for the linear divergence where the coefficient e Nω depends on the gluon action, the steps of stout smearing, N ω , and the value of ω. We separate the coefficient of the linear divergence in the absence of stout smearing, e 0 . Here we use Iwasaki gluons and find e 0 = −12.98. The numerical values of e Nω are shown in Table I for a few choices of ω and N ω ∈ [1,10]. Given the gauge invariance of the linear divergence, the Landau-gauge statements of the previous paragraph, regarding dependence on ω and N ω are directly applicable here. The values ω = 0.1, 0.125, 0.15 have been implemented in numerical calculation of matrix elements of nonlocal operators [41,44,45]. The values ω < 0.1 are included for pedagogical reasons. The data can be compared in different ways leading to various conclusions. Firstly, for all cases, the increase of stout steps leads to an increase of the e Nω value (see, also, left panel of Fig. 8). Because of the different sign of e Nω and e 0 , this leads to a suppression of the linear divergence contribution in the matrix elements. This corroborates the use of multiple smearing steps in physical matrix elements to suppress the linear divergence in the renormalization functions. Secondly, very small values of ω lead to small additions to e 0 , and, thus, the stout smearing does not have any practical benefit. Thirdly, e Nω is not too sensitive to the value of ω when chosen in the range [0.1 − 0.15], which are the values employed in non-perturbative calculations of physical matrix elements. Lastly, the value of e Nω demonstrates convergence as the number of steps increases. This effect can also be seen in the right panel of Fig. 8 where we plot the ratio N ω /N ω−1 . The ratio approaches unity, typically after five steps of smearing for ω = 0.05 − 0.15. While the setup and the conclusions described above are valid in one-loop perturbation theory, they are also indicative for the non-perturbative behaviour.

VII. SUMMARY
In this work, we propose an improvement to RI-type prescriptions for the renormalization of fermion bilinear nonlocal operators containing a straight Wilson line. The method is inspired by its successful implementation to local operators [29][30][31][32][33], and can be generalized to any operator, for instance, nonlocal operators with staple-shaped Wilson line, or gluon operators. The improved renormalization scheme is applied to non-perturbative vertex functions and uses results of the operators' Green's functions from a perturbative calculation in lattice QCD to all orders in the lattice spacing, O((g 2 ) n a ∞ ). The approach is applicable to any order in perturbation theory, but we find that the one-loop level (n = 1) is sufficient to reduce the bulk of finite-a effects. While the method utilizes perturbation theory, the renormalization functions maintain their non-perturbative nature: the subtraction of finite-a effects as outlined in Eq. (11) can be interpreted as an additional finite renormalization containing only higher orders in a.
In our proof-of-concept calculation we implement the improved scheme to two N f = 2 ensembles of twisted mass fermions with a clover term, and Iwasaki gluons. The ensemble has β = 2.10 (a = 0.0938 fm) and lattice size 24 3 × 48. Here, we present results for the ensemble with a pion mass of 235 MeV, which was used in Ref. [41] for the non-perturbative renormalization of PDFs. We note that our existing perturbative calculation is easily adaptable to any gluon action (Plaquette, tree-level Symanzik, Iwasaki and tadpole-improved Lüscher-Weisz (TILW)), and Wilson/clover fermions. Also, the Dirac structure of the operator, Γ, and the parameters g 2 , c SW , the gauge fixing parameter α, (ap) 2 are chosen at the last stage of the calculation.
We find that the proposed renormalization scheme (Eq. (11)) is advantageous compared to the purely nonperturbative one (Eq. (9)), as it has better control of unwanted finite-a effects. As a consequence, one can include in non-perturbative calculations a wider range of renormalization scales, (aµ 0 ) 2 , as well as P 4 values beyond the most democratic ones (P 4 = 0.25), a measure of the Lorentz non-invariant finite-a contributions. For example, in the calculations of the renormalization of nonlocal operators done so far [41,44,[46][47][48][49][50][51][52][53] only about 10 values of the scale were used covering a small (aµ 0 ) 2 interval , typically [2 − 4]; These momenta are chosen based on the P 4 < 0.28 criterion. Here, we are able to include in our non-perturbative analysis 38 values of (aµ 0 ) 2 that cover a range of about [1 − 7]. The improvement of the renormalization functions works well even at P 4 around 0.4 for z ≤ 6, which is a significant advantage. However, caution is needed because of the limitations of perturbation theory that are manifested at high values of z. In this analysis we find improvement for z up to 0.65 fm. Another conclusion from this calculation is that the subtraction of O(g 2 a ∞ ) terms has different effectiveness in the real and imaginary parts of the renormalization function with respect to z. In particular, the real part is still improved for z > 0.65 fm, while the imaginary part worsens. With this in mind, one can employ an alternative scheme to subtract the artifacts, which has different z max for the real and imaginary parts.
As mentioned in the main text, the subtraction of the finite-a terms is done in the absence of stout smearing, which is computationally more expensive due to a fast increase of additional terms in the Green's functions. In this work we lay the foundation for a calculation that contains stout smearing. Our proposal gives, for the first time, a prescription to include an arbitrary number of stout steps with the same computational cost as a single step. This is applicable to cases involving Feynman diagrams with up to two gluons that will be contracted to each other (see, e.g., d4 in Fig. 1), like the cases presented here. This is very powerful, as the non-perturbative calculations are performed with 5-20 stout steps, which is prohibitive in lattice perturbation theory without the approach of Sec. VI.
The applicability of the method can encompass a wide range of operators for both quark and gluon external fields, local and nonlocal cases. For gauge-invariant nonlocal operators, other choices of Wilson lines can be accommodated. Perturbative and non-perturbative results exist in the literature, such as nonlocal operators with staple-shaped Wilson lines [54][55][56] and could be an interesting extension of this work. Another direction that can be immediately pursued is the incorporation of the improved RI-scheme within the hybrid scheme [26], in which an RI-type prescription can be employed in the small-z region.