High-loop perturbative renormalization constants for Lattice QCD (III): three-loop quark currents for Iwasaki gauge action and n_f=4 Wilson fermions

This is the third of a series of papers on three-loop computation of renormalization constants for Lattice QCD. Our main point of interest are results for the regularization defined by Iwasaki gauge action and n_f=4 Wilson fermions. Our results for quark bilinears renormalized according to the RI'-MOM scheme can be compared to non-perturbative results. The latter are available for Twisted Mass QCD: being defined in the chiral limit, renormalization constants must be the same. We also address more general problems. In particular, we discuss a few methodological issues connected to summing the perturbative series such as the effectiveness of Boosted Perturbation Theory and the disentanglement of irrelevant and finite volume contributions. Discussing these issues we consider ont only the new results of this paper, but also those for the regularization defined by tree-level Symanzik improved gauge action and n_f=2 Wilson fermions, which we presented in a recent paper of ours. We finally comment to which extent the techniques we put at work in the NSPT context can provide a fresher look into the lattice version of the RI'-MOM scheme.


Introduction
Numerical Stochastic Perturbation Theory (NSPT [1,2]) can be a powerful tool to address perturbative computations in Lattice QCD up to an order which would be impossible to attain with standard, diagrammatic approaches. A few years ago the Parma group applied NSPT to get three-(and even four-) loop Renormalization Constants (RCs) of finite quark bilinears in the scheme defined by Wilson gauge action and Wilson fermions [3]. Very recently, [4] provided in turn both finite and logarithmically divergent three-loop RCs for currents in the scheme defined by tree-level Symanzik improved gauge action and two flavors of Wilson quarks. The inclusion of divergent RCs was made possible by the method first introduced in [5,6]: when an anomalous dimension is in place, finite size effects can be important in NSPT computations and they have to be carefully taken into account. The main result of the current paper is the computation of quark currents RCs in the regularization defined by Iwasaki gauge action and four flavor of Wilson quarks (quenched computations will be reported as well, to enable a comparison). Preliminary results were quoted in [7]. For a complete discussion of our methodology the reader should refer to [4], which has been largely devoted to discuss in some detail the NSPT approach to the computation of renormalization constants (with a main emphasis on the control over finite lattice spacing and finite volume effects).
Both in the case of n f = 2 tree-level Symanzik gauge action and in the case of n f = 4 Iwasaki gauge action, our results can be compared with analogous non-perturbative computations for Twisted Mass fermions [8,9] (the renormalization scheme is massless and thus RCs are the same). In order to do that, perturbative series have to be summed. An important goal of this paper is a discussion of the issues that are related to summing PT series for Lattice QCD.
The overall structure of this paper is as follows: • Section 2 presents an overview of our methodology. It is mainly intended to allow the reader to go through the paper without having to refer to other sources.
• In Section 3 our results for Z S , Z P , Z V , Z A for Iwasaki gauge action and n f = 0, 4 Wilson fermions are presented.
• In Section 4 we address the issue of summing the series, and in particular we deal with the explicit disentanglement of irrelevant (finite a) contributions, which is possible once also finite volume effects have been corrected for. We take into account not only the results of the current paper, but also those of [4] (i.e., we compare the two regularizations).
• Section 5 contains a discussion of different ways of summing the series. The effectiveness of Boosted Perturbation Theory is discussed; it turns out that this is relevant in particular for Symanzik action.
• In Section 6 we briefly discuss to which extent our approach can provide a contribution for an overall better understanding of the lattice version of the RI'-MOM scheme.

Three loop renormalization constants in NSPT
In this section we provide a brief account of our computational strategy. This is basically a summary of the discussion of [4], to which the interested reader is referred for an in-depth description of our method.

3-loop RI'-MOM lattice computations
The lattice is a suitable regulator for the RI'-MOM renormalization scheme [10]. The definition of the latter for quark currents starts from the computation of Green functions on external quark states at fixed momentum p G Γ (p) = dx p| ψ(x)Γψ(x) |p .
The G Γ (p) are then amputated to get vertex functions (S(p) is the quark propagator) By projecting on tree-level structure one gets the quantities O Γ which enter the definition of the currents renormalization constants By choosing different Γ one obtains the different currents, e.g. the scalar (identity), pseudoscalar (γ 5 ), vector (γ µ ), axial (γ 5 γ µ ). The master formula (Eq. (1)) is defined in terms of the quark field renormalization constants which in turn reads We adhere to the standard recipe of getting a mass-independent scheme by defining everything at zero quark mass.
A main point in our strategy is to get the (divergent) logarithmic contribution to the renormalization constants from continuum computations: NSPT is only in charge of reconstructing the finite parts. The typical renormalization constant we want to compute (in the continuum limit) reads where the lattice cutoff (a) is in place and the expansion is in the renormalized coupling. Divergencies can show up as powers of l = log(µa) 2 . By differentiating Eq. (3) with respect to l one obtains the expression for the anomalous dimension γ = 1 2 d dl log Z whose expansion can be read from continuum computations [11] This is a scheme dependent, finite quantity, with no dependence on the regulator left. By imposing that the expression we get by differentiating Eq. (3) matches Eq. (4) we can obtain the expressions of all the d (i>0) n (n ≤ 3), which are thus expressed in terms of the γ m≤n , the d (0) m≤n and the coefficients of the β-function; the latter come into place since part of the dependence on µ in Eq. (4) is via the coupling α(µ).
In the above discussion there was no reference to a (covariant) gauge parameter λ. This is legitimate, since we compute in Landau gauge, i.e. λ = 0. In a generic (covariant) gauge, one has a dependence on λ entering Eq. (3). Moreover, the gauge parameter anomalous dimension comes into place in linking Eq. (3) to Eq. (4). Since the non trivial dependence on the gauge parameter anomalous dimension is itself proportional to λ, all this is immaterial in Landau gauge: if one keeps track of all the λ-dependence and then puts λ = 0 one gets the same result which is got by ignoring λ from the very beginning.
In our computations the Zs are expressed as expansion in the bare lattice coupling α 0 Eq. (5) is obtained from Eq. (3) by plugging into the latter the matching of the renormalized coupling to the lattice bare one.

2-loop matching of α IWA to continuum
The matching of α IWA to a continuum coupling is only known to one-loop [12]. Since we need a two-loop matching to get Eq. (5), we had to compute it. This was done by first matching α IWA to an intermediate scheme, which was chosen to be a potential scheme. The matching of the latter to MS is known [13] and the results for the anomalous dimension we can read from [11] are obtained as expansions in α M S . Thus, computing the matching of potential coupling α V to α IWA is a possible solution. Here and in the following our notation only enlightens the dependence of the scheme on the gluonic action: the dependence on Wilson fermions has to be assumed as well, when we refer to the four flavors case. The strategy of the computation is that of [14,15]. The interested reader can find more information on technical details both in [4] and in [16].
We started from the NSPT computation of Wilson loops W (R, T ) from which we got Creutz ratios A potential for static sources at distance r = Ra can now be defined and a coupling out of it according to where one can see that in a lattice regularization a residual mass δm comes on top of the coupling. Here we need to rely on an approximation, since we can not compute the limit in Eq. (6). As a consequence of the same observation, our results are not in the continuum limit. Despite this, we could obtain a decent estimate of the matching that in perturbation theory reads (r = Ra) where the expansion coefficients are a function of scale parameters Λ and coefficients of the Reconstructing one-loop result was a check that the procedure is viable, and at two-loop we could finally obtain where the new piece of information is contained in the quantity X: in the following we will refer to the latter.

3-loop critical mass
Staying at zero quark mass in our three-loop NSPT computation requires the knowledge of the Wilson fermion critical mass at two-loop, which is known from the literature [17].
From now on, we switch to β −1 as the expansion parameter for our results. Also, we introduce a hat notation to denote dimensionless quantities, e.g.p = pa (if needed, explicit factors of a will be later singled out).
The critical mass is computed from the inverse quark propagator (m W (p) = O(p 2 ) is the irrelevant mass term generated at tree level) More precisely, in the self-energyΣ(p,m cr , β −1 ) we single out the components along the the (Dirac space) identity, the one along the gamma matrices and the irrelevant one along the remaining elements of the Dirac basiŝ Σ(p,m cr , β −1 ) =Σ c (p,m cr , β −1 ) +Σ γ (p,m cr , β −1 ) +Σ other (p,m cr , β −1 ), The critical mass can be read fromΣ c 1 at zero momentum The known one-and two-loop values of the critical mass were inserted (as counterterms): this is enough to have massless quarks at the order we are interested in (three-loop). The novel three-loop result for the critical mass is not relevant for the computations at hand: it is simply a byproduct.
Computations were performed on different lattice sizes: 32 4 , 24 4 , 20 4 , 16 4 , 12 4 . Left panel of Fig. 1 shows the three-loop computation ofΣ c at different values of momentum on a 32 4 lattice, in the n f = 4 case. We are interested in the zero momentum value, which can be got by fitting our observable as an expansion in hypercubic invariants, e.g.
This is a general feature of all our computations. On each lattice size we got a different value and an infinite volume result could then be obtained by extrapolation. Right panel of Fig. 1 displays results plotted as a function of N −2 , which is the power that best fits our data. N = L/a is the only significant quantity in the NSPT context (there is no value one can attach to the lattice spacing a). The infinite volume extrapolation was first fitted by keeping only the single power −2. We then checked that the central values of fits performed adding other powers were consistent with that result, within the error of the latter. This procedure, thought limited by the number of available sizes, proved to be accurate enough, as we could check at one-and two-loop level, for which the expected zero value of the critical mass was obtained in the infinite volume limit. Our final results for the three-loop contribution to the critical mass arem

Fitting irrelevant and finite volume effects
In Eq. (1) currents renormalization constants are defined in terms of the quark field renormalization constant. The latter can be computed -see Eq.
(2) -from the quark self-energy, more precisely from its component along the gamma matrices, which at any finite value of the lattice spacing readŝ Notice the tower of irrelevant contributions which go on top of the one expected in the continuum limit. All these contributions are contained in the definition of Σ γ (p, ν). In the latter one recognizes a dependence on the direction ν, which comes via the dependence on the length |p ν |. In the continuum limit this dependence drops out and, once one subtracts the logarithmic contribution discussed in subsection 2.1, the value of the finite part of Z q (µ = p) is given by lim a→0Σ . This observation on the dependence on ν of Σ γ (p, ν) becoming immaterial in the continuum limit has to be be born in mind also in the following, e.g. in Eq. (15).
The second ingredient in Eq. (1) is given by the quantities O Γ . These have their lattice counterpartsÔ For the vector and axial currents we eliminate dependences on directions like the one we have just discussed in the case of Σ γ (p, ν), e.g.
The reason for getting rid of this dependence in this case while retaining it in the case of Σ γ (p, ν) will be clear in a moment.
Our NSPT computations are performed on different, finite lattice sizes N = L/a (our lattices are always isotropic); therefore one must expect finite size corrections. On dimensional grounds, we can expect a dependence on pL. Since we want to compute the currents renormalization constants in both the continuum and the infinite volume limit, we need to take two limits. This is done in the following form Eq. (15) is our key formula and deserves a few comments: • We only compute the finite parts. As it was made clear in subsection 2.1, we know all the relevant logarithms entering the quantities we are concerned with. This means in particular that we can subtract their contribution: this is the meaning of the subscript . . . | log subtr in the definition of O Γ (p, pL, ν).
• Since in Eq. (15) we take the limits a → 0 and L → ∞ and we subtract the logs that come from the anomalous dimension γ O Γ , Σ γ (p, pL, ν) reconstructs the contribution of Z q (and in this limit the dependence on ν drops out). Notice that this is true because of the ratio that is taken. In order to determine Z q itself one should look at Σ γ (p, pL, ν) alone and perform the subtraction of different logs (i.e., those connected to the quark field anomalous dimension).
• On a fixed lattice volume, the a → 0 limit of O Γ (p, pL, ν) can be evaluated by computing the quantity for different momentap and fitting the results in terms of hypercubic invariants, e.g. those listed in Eq. (13). The possible terms are dictated by symmetries of both Σ γ (p, pL, ν) andÔ Γ (p, pL) (a formal power counting fixes how many terms one should retain).
• We want to account for the limits a → 0 and L → ∞ simultaneously. This is done by computing the quantity O Γ (p, pL, ν) on different volumes and performing a combined fit. The combined fit is made possible by defining finite size corrections according to where the main rationale for the last (approximate) equality is that we neglect corrections on top of corrections. Since p µ L = 2πnµ L L = 2πn µ , there is only one finite size correction for each 4-tuple {n µ | µ = 1, 2, 3, 4} and no functional form has to be inferred for the correction.
All in all, a prototypal fitting form of ours reads where in order to make things easy we limited to a very moderate order in a (actually, less than what we use in realistic fits). The relevant contribution to the finite part is c 1 . We recall once again that this is a combined fit to data taken on different lattice sizes, with the same ∆ O Γ (pL) applying to all the data corresponding to the same momentum 4-tuple {n 1 , n 2 , n 3 , n 4 } (resulting in different values of momenta on different lattice sizes). Notice that the inclusion of ∆ O Γ (pL) makes the fit not constrained with respect to an overall shift. This is cured by including in the fit a few measurements (of the order ∼ 1, 2, 3) taken on the largest lattice in the high momentum region: these are assumed free of finite size effects and act as a normalization point.
Finite size effects can be easily spotted in the left panel of Fig. 2, where we plot oneloop O S (p, pL, ν) on both 32 4 (black symbols) and 16 4 (red, filled symbols), in the n f = 4 case. One can see that data are arranged in families (different symbols): this is a direct consequence of the dependence on ν one can see e.g. in Eq. (17). All in all, there is one family for each length |p ν |. Since finite size effects are there, families do not join smoothly across different lattices. On the right panel of Fig. 2 one can inspect how things change taking into account finite size corrections ∆ O S (pL): families do join smoothly. This family mechanism provides a very effective handle to detect finite size effects: this is the reason for retaining the ν dependence in the definition of O Γ (p, pL, ν); to be definite, we retain it in the numerator (i.e., in the contribution connected to Z q ); keeping it also in the denominator would result in a too odd fitting form.

Results
In Table 1 we give a brief account of our statistics. One can see that computations were performed on five different sizes (32 4 , 20 4 , 24 4 , 16 4 , 12 4 ) for both n f = 0 and n f = 4. The latter are our main point of interest. While with only two values of n f we can not determine the coefficients of the n f dependence 2 , it is interesting to have at least an indication of how sensitive results are to the number of flavors.
Since NSPT requires the (order by order) integration of the Langevin equation, a finite order integration scheme for this stochastic differential equation is needed: the Euler scheme is our choice here. An → 0 extrapolation (a linear one, in the case at hand) is needed to remove the effects of the finite time step . In Table 1 we provide the statistics collected on each different size. Notice that the values of are chosen different for n f = 0 and n f = 4 (see the discussion in [2]).
Configurations were saved on which we can still measure different observables. They were saved with frequencies which were chosen having a rough analysis of the autocorrelations in place. The analysis of the different observables is in charge of dealing with the residual autocorrelation effects.
In Table 2 we report the coefficients of the three-loop expansion of Z S , Z P , Z V and Z A 3 . Results are reported for both n f = 4 and n f = 0. We remind the reader that   12  210  210  208  220  209  209  16  179  185  148  111  119  116  20  84  84  82  71  70  70  24  72  74  79  54  54  52  32  49  47  49  42  46  45 the expansion parameter is β −1 . We stopped at three-loop order given the knowledge of anomalous dimensions which we can get from [11]. For finite quantities there is in principle no limitation (other than practical ones dictated by statistics). Notice that [19] could now open new opportunities for even higher order computations. We quote the analytical oneloop results [12]: the comparison is a first proof of the effectiveness of our method. The reader could notice that we have somehow less systematic deviations from analytic results here than in [4]. This is due to the fact that we took the normalization points for finite size effects at slightly higher values of momentum: see the discussion after Eq. (17). A more stringent confirmation of our results comes from the fitting of irrelevant contributions. The latter is a key ingredient of our approach, since fitting irrelevant contributions compliant to lattice symmetries is our handle on continuum limit. A comparison to results in [18] made us confident in our results: the leading irrelevant terms that we fitted are consistent with the diagrammatic results.
The errors we quote are dominated by the stability of fits with respect to the change of fitting ranges, functional forms, number of lattice sizes simultaneously taken into account. Notice that three-loop results for Z S and Z P have an extra source of error in the indetermination in the coupling matching parameter X (see Eq. (9) and the discussion over there).
All in all, our new results, i.e. two-loop and three-loop contributions to Z S , Z P , Z V and Z A for Iwasaki action seem to be quite moderate, in particular two-loop: taking into Table 2: One-, two-and three-loop coefficients of the renormalization constants for quark bilinears for both n f = 4 and n f = 0. Expansions are in β −1 . One-loop analytical results are reported for comparison. account the typical values of β that are relevant to numerical simulations (β ∼ 2), three-loop contributions are typically larger than two-loop. Coefficients themselves are smaller than the ones found in the case of tree-level Symanzik improved action [4], but this is not per se any significant. First of all, the comparison in magnitude of two-and three-loop coefficients has to be corrected for the different β value regimes one is interested in (there is roughly a factor of 2). What is even more important is the weight of two-and three-loop contributions with respect to the leading one: what really makes the difference in between the two different regularizations is the relative weight of one-loop contributions themselves.
Another general feature that emerge from our computations is that irrelevant corrections from hypercubic invariants which are not O(4) invariant appear to be in general quite significant for this action. All this is of course a numerical accident, but it is a relevant one when it comes to summing the series and assessing irrelevant effects.

Summing the series
We now come to the issue of summing the series. In this and in the following Section we will deal not only with the results of this paper, but also with the ones for the regularization defined by tree-level Symanzik improved gauge action and n f = 2 Wilson fermions, i.e. those of [4].
For both cases one can compare perturbative and non-perturbative results, which for Iwasaki can be found in [9]. As already said, [9]  where the first error is the statistical one, while the second is a rough estimate of the truncation effects, which we simply take as the highest order contribution. The latter recipe is the conventional one: one is of course well aware of the its roughness, even if making use of it at three-loop level is more than what is usually done. On the other side, we have already made the point that two-loop contributions are indeed small for the Iwasaki case (and indeed even smaller for the n f = 4 than for the n f = 0 case). Since three-loop contribution is thus relatively important, the net effect is an estimate of truncation errors which is quite large. The other finite renormalization constant is Z A , for which we obtain One can directly inspect a fair agreement with the results of [9]: basically the errors that result from our procedure make the perturbative and non-perturbative results fully consistent. Actually the (smaller) statistical errors would be enough to obtain a substantial agreement with non-perturbative results in the case of the finite renormalization constants. Notice that non-perturbative results in [9] are presented in two variants, referring to the different prescriptions "M1" and "M2" which the authors discuss. Here it suffices to say that the differences are to be ascribed to different treatments of irrelevant effects (and thus they are systematic effects): in general the two methods differ more for the divergent than for the finite renormalization constants. We will have more to say on irrelevant effects later in this Section and then again in Section 6.
We now move to the results we got for the regularization defined by tree-level Symanzik improved gauge action and n f = 2 Wilson fermions (i.e. those in [4]). We stress once again that in this case the two-and three-loop coefficients are larger, but this should be corrected by taking into account the regime of β one is interested in (typical values of β for tree-level Symanzik are roughly double of those for Iwasaki). Moreover, convergence properties of the series are dominated by the relative weights of one-loop and higher orders contributions. The results we obtain summing the series we computed in [4] can be compared to the nonperturbative ones in [8]. In this case we make our comparison at the largest value of β which is discussed in [8] (the reason for this will be clear in a moment). For Z V our results sum to Conventions with errors are the same as before. In this case deviations are manifest, in particular for Z S and Z P . This in the end does not come as a surprise, given the observations we have already made: convergence properties are strongly controlled by the relative weight of one-loop and higher-order contributions. This is the reason for not attempting to sum the series at values of β smaller than the largest one. While there is a tendency to converge for finite quantities, logarithmically divergent constants are fairly away from each other in the perturbative and non-perturbative computations. This clearly motivates the step forward of summing the series in different couplings, which will be addressed in the following Section. Before we move to that issue, we present a first discussion of how we can assess the impact of irrelevant effects once we sum the series.
In Figure 3 we plot the results we obtain when we sum Z V for the Iwasaki (left panel) and the Symanzik (right panel) case (values of the coupling are once again β = 2.10 and β = 4.05 respectively). The figures have to be understood as follows. We have at our disposal the fits for the one-, two-and three-loop quantities O loop order). Notice that this time we do not write down any pL dependence because that has been eliminated by fitting the ∆ O V (pL). The dependence onp is the irrelevant one: this is trivially true for finite quantities, but would be true also for divergent constants since we always compute the finite parts only. The values one is interested in (the finite contributions) are marked by the red, continuous lines: these are just the values for Z V we reported before. The reason why the values marked by the red, continuous lines can be slightly different from the central values we wrote down a few lines above is easy to understand: these results come from real, individual fits, which always return both the finite and the irrelevant parts. Deviations from the central values give an idea of the typical deviation of results coming from one fit or another (in other terms, these deviations are of the order of the statistical errors that we quoted before). On top of the red, continuous lines we plot blue, empty circles. These are just the sum of both relevant and irrelevant contributions. In other terms, these are the values one obtains at different momenta and at finite values of the lattice spacing according to our perturbative computations. Notice that in abscissa we report values of momentum in dimensionless units (in other terms, there is no value for the lattice spacing involved). Figure 3 displays the quantity 1 , ν), i.e. this time we average on directions, which is the common practice. Notice that irrelevant effects come out of a fit that is necessarily an effective one: we have to stop at a given order in the lattice spacing. We stress nevertheless that the fit is performed at fairly large orders (typically a 6 ) and at three-loop level.
One can easily see how different is the impact of violations of (continuum-like) rotational symmetry in the two cases. It is true that one often tries to minimize these effects by a convenient choice of the momenta. One should nevertheless keep in mind the trivial observation that the amount of violation is not decided by the choice of momenta: one should try in any case to fit terms compliant to the lattice symmetries.
for Symanzik action at β = 4.05. The red, full circles mark the contributions coming from terms which are not depending on the length |p ν | (see the text for a few more comments). In other words, in the right panel the blue, empty circles single out the families. None of them smoothly guides the eye to the extrapolated finite result.
In Figure 4 we compare (in the case of Symanzik at β = 4.05) the observable relevant for computing Z V averaged (left panel) or not averaged (right panel) on directions. In the left panel one simply recognizes the same quantity which is plotted in Figure 3. In the right panel we plot 1 , and there is one subtlety which is related to the curves we have already referred to as families. In both panels, the red, full circles are the contributions which are not dependent on the length |p ν |. One can understand what we mean by referring to Eq. (14): red, full circles single out the contribution of terms like theΣ (0) γ (p). One can see that in the right panel the blue, empty circles do not point in a trivial way to the result one is interested in. In other terms, when we keep the families structure, only the red, full circles are the ones smoothly guiding the eye to the correct extrapolated result.

Summing the series in different couplings
The Symanzik case displayed not so brilliant convergence properties. Thus, that is the prototypal situation in which one would like to go for what is usually, generically referred to as Boosted Perturbation Theory [20]. One re-expresses the series as expansions in different couplings, of course looking for better convergence properties that in the case one starts with. Often one deals with this having only a one-loop result available. As discussed in [3], this is at risk of being an empty exercise. At one loop, nothing changes but the value of the coupling itself. So, the effectiveness of the procedure relies on the optimal choice of coupling and scale that are really relevant for the computation at hand. Actually this choice has to be regarded as so good that one-loop captures essentially the complete result. This does not need to hold true and can be strictly speaking only assessed a posteriori. Only having at least a two-loop result available one can inspect how the series actually reshuffle and one can hope to learn something more on the convergence properties. Our question is: can a three-loop computation be reliable enough to gain solid, new pieces of information?
We here compare results obtained as expansions in the couplings which were also used in [3], i.e.
P is the basic 1 × 1 plaquette, for which we do have an expansion in β −1 ; so, we can work out the expansions we are interested in. x 2 and x 3 are quite popular as boosted couplings. In the end, we want to see whether results coming from summing series in different couplings do or do not all approach the same result. The definition of x 1 can be useful with this respect. Convergence properties in the Iwasaki computations are fairly good in the original coupling; we will focus the case of Symanzik action, looking for better convergence. There is an overall ambiguity we have to live with: we do not have non-perturbative simulations in the same setting we are dealing with (Symanzik action and n f = 2 Wilson fermions).
In view of this limitation, we have no non-perturbative value for the different couplings.
We have indeed estimates which come in turn from summing perturbative expansions of the plaquette (actually at higher orders than three-loop): these are the values we plug in. This is admittedly a limitation. Still, if we take into account the order of magnitude of the error one can attach to the value of the coupling, it turns out that this is dominated by the other errors, typically the truncation errors which are still the dominant ones. The latter are estimated as done previously (i.e., as the highest order contribution) and will be the only ones reported in the following.
Let's start from Z V , for which we previously got the result Z V (β = 4.05) = 0.710(2)(28). Switching from x 0 to x 1 we obtain Z V (β = 4.05) = 0.686(21) x 1 -expansion Next step is the expansion in x 2 (notice that the value of the couplings are getting larger and larger as we proceed), from which we get We notice that results are quite close to each other and they both approach the result of [8]. We get even closer when we switch to x 3 , an expansion in which returns Z V (β = 4.05) = 0.661(55) x 3 -expansion While the central value is now literally on top of the non-perturbative result, the error has become pretty large. This is simply the effect of the fact that the series has started oscillating: already at one loop one gets essentially the result 0.66, and then two-and threeloop contributions basically cancel each other. We proceed to Z A , for which in the original coupling we reported Z A (β = 4.05) = 0.788(2)(18)), and for which we now get respectively: Z A (β = 4.05) = 0.773 (12) x 1 -expansion Z A (β = 4.05) = 0.775 (9) x 2 -expansion Z A (β = 4.05) = 0.763(26) x 3 -expansion Once again, in the last case the series has already started oscillating. All in all, it is fair to say that results for finite constants display a tendency to get closer to the non-perturbative ones. Actually, Z V changed more than Z A , which is good, since the latter was deviating more than the former from non-perturbative results.
We proceed to the logarithmically divergent renormalization constants. For the scalar current the result obtained in the β −1 -expansion reads Z S (β = 4.05) = 0.753(4)(30); now we get Comparing to the results in [8], one can still see quite important discrepancies for divergent renormalization constants, which did not hold true in the case of finite constants. It could well be that one simply needs more terms to definitely assess the convergence properties, but there is another issue which could be considered. We have already noticed that "M1" and "M2" results in [8] differ much more in the case of Z S and Z P than in the case of Z A and Z V . One method tries to gain more information from the lower momenta region than the other. To be more precise, one simply subtracts the known leading one-loop a 2 irrelevant effects and look for a plateaux region, while the other tries to fit extra irrelevant effects in the lower momenta region. This region is just the theater of a subtle interplay of UV and IR effects: from one side higher powers of pa are suppressed (and this is good to assess irrelevant contributions), but from another side that is just the region which is prone to suffer from finite size (IR) effects.
All in all, Boosted Perturbation Theory apparently solves the problem of the discrepancies in between perturbative and non-perturbative results for Z V and Z A . This sounds good, also in view of the fact that different boosted couplings basically point to consistent results. Discrepancies are still there for Z S and Z P . While there is of course the possibility that even higher order terms should be included, there is another explanation that could hold true. Given the interplay of IR and UV effects, there is a possibility that non-perturbative computations could suffer from finite volume effects. These effects are not expected to be the same that we get (and correct for) in our NSPT setting, but could be possibly assessed: more on this in the following Section.
As a final comment, we go back to the Iwasaki case, for which basically there was no compelling reason to go for boosted couplings (of course one could nevertheless do it). We stress that the latter observation could be done only in view of the control of the series at three-loop level.

Some general remarks on lattice RI'-MOM
There is something interesting that one can learn from our computations, not only with respect to a comparison of perturbative and non-perturbative results.
First of all, we put forward a method to assess (the possible presence of) finite size effects. One can see that there is in principle no reason why one should not attempt the same in the non-perturbative case. We are actually working on this [21].
Moreover, the high-loop computations which are enabled by NSPT can provide a new handle to correct non-perturbative computations with respect to irrelevant contributions. We have briefly sketched this in [22]. Quite interestingly, another group is working on the same ideas [23]. Basically this amounts to the following simple recipe: • One needs both a high-order NSPT computation and a standard non-perturbative computation.
• First of all, one should try to assess the possible presence of finite size effects in both cases. We stress once again that these do not need at all to be the same. Once assessed, they should be corrected in both computations.
• Once both results are corrected for (possible) finite size effects, one can take the irrelevant effects as estimated via the fitting procedure we described (here and) in [4] and subtract them from the non-perturbative data.
Subtracting irrelevant effects is by now a common practice. There are many approaches to this, requiring different combinations of perturbative computations and fitting of terms compliant to the lattice symmetries: see [24] for a recent contribution. In the end, our proposal is basically yet another variant, whose merits are worth investigating.

Conclusions and prospects
This work is a little landmark at the end of a path that we took a few years ago. The point we wanted to make is that the three-loop computation of Renormalization Constants for Lattice QCD is a realistic goal. There is in principle no sharp constraint on computing finite constants, while for logarithmically divergent ones there is a limit because continuum computations are available at three-loop order in the RI'-MOM scheme. These results make it possible to derive the leading logarithmic contributions one has to account for in the lattice regularization of the same RI'-MOM scheme. As for the finite parts, Numerical Stochastic Perturbation Theory can do the job. All this is under control because we can assess both finite lattice spacing (UV) and finite volume (IR) effects.
As a general conclusion, it is fair to say that the NSPT approach to the computation of Renormalization Constants for Lattice QCD can provide at least two valuable contributions. First of all, it is a completely independent approach with respect to non-perturbative computations, with different systematic effects. From another point of view, NSPT techniques provide a new method to correct non-perturbative computations with respect to irrelevant contributions.
Last but not least, the method we suggested for the correction of finite size effects could be useful in non-perturbative cases as well.