Discretization effects on renormalized gauge-field Green's functions, scale setting and gluon mass

The $SU(3)$ gauge-field propagators computed from the lattice have been exhaustively used in the investigation of the low-momentum dynamics of QCD, in a judicious interplay with results from other nonperturbative approaches, and for the extraction of fundamental parameters of QCD like $\Lambda_{\overline{\rm MS}}$ as well. The impact of the discretization artifacts and their role in the extrapolation of the results to the continuum limit have not been fully understood so far. We report here about a very careful analysis of the physical scaling violation of the Landau-gauge propagators renormalized in MOM scheme and the Taylor coupling, steering us towards an insightful understanding of the effects from discretization artifacts which makes therefore possible a reliable continuum-limit extrapolation.

The key importance of these outcomes makes particularly relevant a careful examination of the impact of regularization artifacts on the lattice QCD Green's functions in Landau gauge, as the one performed in Ref. [26]. The role played by the discretization artifacts, crucial for the extrapolation of results to the continuum limit and for the extraction of QCD parameters therefrom, was not fully understood from the analysis therein performed. Two following publications, a "Comment-to" [68] and its corresponding "Reply" [69], elaborated further on this issue, mainly in connection with the problem of the lattice scale setting, but failed to settle properly the question about continuum extrapolation.
The object of this paper is, precisely, the thorough examination of the scaling violations for Landau-gauge gluon and ghost propagators, after MOM renormalization, and for the running coupling in Taylor-scheme that can be computed from them. Indeed, an ace of our analysis comes from the deep connection between the bare propagators and the renormalized coupling: the latter results when the former are appropriately combined and multiplied by the bare gauge coupling, which is a parameter of the discretized action fixed by the lattice setup [63,70]. Therefore, a consistent analysis of the three quantities is very demanding and we capitalize on the successful description of results obtained from several large-volume lattices, with different set-up's, made only possible by the use of the quenched approximation. We can thus get an insightful understanding of the regularization cut-off effects for gluon and ghost Green's functions, paving the way towards a reliable and very precise extrapolation to the continuum limit.

II. DISCRETIZATION ARTIFACTS ON RENORMALIZED GREEN'S FUNCTIONS
Aiming first at arguing on general grounds, let Γ(k 2 , a) be the bare dressing of the QCD two-point Green's function for either the gluon or the ghost gauge fields; k being the propagated momentum and a standing for a regularization cut-off which drops out by approaching zero (e.g., in lattice regularization, a is the length-dimension lattice spacing; while in dimensional regularization, it corresponds to the dimensionless parameter ε = (4 − d)/2).

arXiv:1809.05776v1 [hep-ph] 15 Sep 2018
The renormalized dressing function will be then obtained by applying first a well-defined subtraction procedure, implying a particular prescription, and removing next the cut-off. Namely, a) is the renormalization constant, defined for a given scheme and at the subtraction point k 2 = ζ 2 . If MOM prescription is considered, any renormalized two-point Green's function is required to take its treelevel value at the subtraction point, i.e. to amount to unity, such that Z Γ (ζ 2 , a) ≡ Γ(ζ 2 , a) and Thus, it becomes manifest that, without taking the explicit limit making the cut-off to drop out, the subtraction procedure cannot generally prevent from the remaining of a residual dependence on the cut-off. Specially in lattice computations, where simulations are carried out in lattices for which the cut-off is fixed by their discretization spacing, one should either try an extrapolation to the continuum limit by the examination of several appropriate simulations, or somehow care about such a residual cut-off-dependence. Let us specialize, for illustrative purposes, to the quenched gluon-propagator dressing function, denoted as D, which reads at the perturbative oneloop level and in lattice regularization as follows where γ 0 = 13/2 corresponds to the one-loop anomalous dimension and g 0 is the bare gauge coupling, related to the lattice spacing. Apart from the O(a 2 k 2 )-terms, respecting the O(4) symmetry, there also appear other more complicate lattice artifacts, owing to its breaking into the symmetry under the action of the isometry group H(4) of hypercubic lattices on the discrete momentum ak µ ≡ 2π/N n µ , with n µ integers and N the lattice volume in units of the lattice spacing, a. These higher-order artifacts, indicated explicitly in Eq. (2.3), can be generally expressed in terms of the invariants a n k [n] = a n 4 µ=1 k n µ for n = 4, 6, 8. In particular, a leading H(4) correction, at the O(a 2 )-order, is proportional to a 2 k [4] /k 2 [60] and is already present in the tree-level gluon propagator computed in a lattice. As will be pointed below, the O(4)-breaking artifacts can be efficaciously cured by applying the so-called H(4)extrapolation [60,61,71] and we will then focus here on dealing with the leading artifacts preserving the O(4)symmetry.
Therefore, after renormalization but before considering the continuum limit, according to Eq. (2.2) with Γ ≡ D, one would be left with where we have also assumed that H(4)-extrapolation has been applied and the O(4)-invariant corrections have been removed. There, c D is a constant and g R is the renormalized gauge coupling. Eq. (2.4) makes explicit a property of the residual Green's functions' cut-offdependence, implicit from Eq. (2.2), namely its vanishing at k 2 = ζ 2 for so to enable both conditions simultaneously: Γ L (k 2 , ζ 2 ; a) = Γ R (k 2 , ζ 2 ) = 1. Beyond perturbation theory, a nonperturbative resummation for the expansion in the coupling at a fixed cut-off would lead to a more general expression that, extended by analogy also to the ghost case, would read (2.5) where Γ ≡ D (gluon) or F (ghost), and C Γ is a squaredmass dimension function vanishing at k 2 = ζ 2 . The nonperturbative emergence in the expansion of Eq. (2.5) of mass scales, as Λ QCD or the gluon mass, enables that momentum and lattice spacing decouple and we can thus extend the one-loop result of Eq. (2.4) and conjecture that is an effective description of the leading a 2 -contribution for the residual cut-off-dependence, where m g stands for the gluon mass which has been strongly argued to emerge as a result of the so-called Schwinger mechanism to saturate the gluon propagator at vanishing momentum and cure the running coupling from the Landau pole [39][40][41][42][43].
A particular fruitful MOM-like renormalization scheme for the strong running coupling results from the threepoint ghost-ghost-gluon Green's function, defined in Landau gauge and at a subtraction point for which the incoming ghost momentum vanishes [63,70]. Namely, the so called Taylor-scheme running coupling, which has been recently shown to be intimately related to the quarkgap-equation interaction kernel in Dyson-Schwinger approach [72] and to a process-independent effective charge built in analogy to the Low-Gell-Mann QED charge [47]. The particularities of the Taylor-scheme kinematics and the Landau gauge makes the coupling to rely only on the ghost and gluon two-point Green's functions [62,63]. It specifically reads Therefore, the conjecture expressed by Eq. (2.6), translated to Eq. (2.8) allows for a simple separation of k-and ζ-dependence so that one is left with (2.9) where c α = 2c F + c D and d α = 2d F + d D , Λ being a mass-dimension parameter which can be derived from a nonperturbative subleading O(a 2 )-contribution in the bare Green's functions that cancels after MOM renormalization in Eq. (2.2), thus not spoiling the condition at k 2 = ζ 2 , but does not cancel in the running coupling definition, Eq. (2.7).

III. LATTICE DATA: THE ANALYSIS
In what follows, we will examine the validity of Eq. (2.9), and the underlying conjecture about the residual dependence on the cut-off expressed by Eqs. (2.5)-(2.6). For so to accomplish, we will directly obtain Γ L (k 2 , µ 2 ; a), and so α L T (k 2 ; a), from several different lattice simulations with different set-up's parameters and evaluate then whether these quantities relate to their continuum counterparts, Γ R (k 2 , ζ 2 ) and α T (k 2 ), as equations suggest.

A. Set-up's and scale setting
We have produced the SU (3) lattice gauge field configurations U µ (x) from the Monte Carlo sampling using the standard Wilson gauge action, where β ≡ 6/g 2 0 (a); and next gauge fixed them to the minimal Landau gauge as explained, for instance, in Ref. [25]. The set-up's parameters can be found in Tab. I. Then, the gauge field is defined as 2) withμ indicating the unit lattice vector in the µ direction. The two-point gluon Green function is then computed in momentum space through the following Monte-Carlo average where λ a stands for the Gell-Mann matrices and the trace is evaluated in color space. Then, the gluon dressing function results from taking the appropriate trace of the propagator, On the other hand, the Landau gauge ghost propagator result from the Monte-Carlo averages of the inverse of the Faddeev-Popov operator, M . Namely, Thus, Eq. (3.5) and Eq. (3.6) define the bare gluon and ghost dressing function which are obtained as the appropriated projection of the lattice propagators. Statistical errors have been derived by applying the Jackknife procedure. More details of the computations can be found in Refs. [25,73]. Before applying the MOM prescription to get the renormalized dressing functions, as above stated, the H(4)-extrapolation is applied to deal with the O(4)breaking artifacts. The dressing functions, being scalar form factors of two-point Green's functions, computed in lattice QCD, are not invariants under O(4) but under H(4) transformations. Therefore, the prescribed recipe implies the average of results obtained for momenta corresponding to the same H(4) orbit (all the lattice four-momenta with the same k [n] invariants) and, next, an extrapolation towards the continuum limit by the subtraction of the O(4)-breaking contributions, fitted as smooth functions of k [n] for all orbits sharing the same k 2 . More details for the H(4)-extrapolation procedure can be found in Refs. [60,61,71]. On top of this, for all the set-up's, we will apply an upper cut in lattice momenta: ka(β) ≤ π/2. a(β)/a(β0) β = 5.6 β = 5.7 β = 5.9 β = 6.02 β = 6.202 [74]  Finally, the scale setting is a key issue for the combined analysis of data resulting from different set-up's with several β's. Specially, such a combined analysis relies on an accurate relative lattice calibration; i.e., the knowledge of the ratios of lattice spacings for any pair of β's. Indeed, a thorough discussion about how deviations in the lattice calibration might impact on the scaling of renormalized propagators has been the object of a Comment [68] to Ref. [26] and its corresponding Reply [69]. With this in mind, we can follow Ref. [75] and write where the coefficients a j , obtained by a fit of Eq. (3.7) to lattice spacings obtained from applying the Sommer's parameter method for 5.7 ≤ β ≤ 6.92 [74], appear gathered in Tab. II. Nevertheless, cut-off effects borrowed by the scale-setting procedure, the determination of the heavy quark potential at intermediate distances in the case of the Sommer's parameter, can induce significative deviations at low β for the lattice spacings obtained with two different procedures. The effect of these deviations will be anyhow cured by the extrapolation to the continuum limit when is properly made. However, as suggested in Ref. [68], the scaling of a renormalized Green's function obtained for different β's, when it exists, provides with a strong criterion guaranteeing the negligible impact on them of cut-effects from the lattice scale setting. This would be an ace for our analysis and we can thus, as done in Ref. [68], assume that O(4)-breaking artifacts dominate the cut-off deviations for the gluon propagator so that, after applying H(4)-extrapolation, the scaling of the gluon dressing function can be imposed as the condition in order to fix the ratios of lattice spacings. As it will be seen below, this assumption underlying the validity of the relative scale setting, i.e. the scaling of the gluon dressing function after renormalization and H(4)extrapolation, can be explicitly confirmed a posteriori.
In so doing, we obtain the results of Tab. II for the ratios a(β)/a(β 0 ) and fit Eq. (3.7) to them, thus obtaining a refined set of coefficients a j , also collected in Tab. II, reliable only from 5.6 ≤ β ≤ 6.2. We have then displayed the results from Eq. (3.7) with the two sets of coefficients a j , from [74] and the one herein obtained, in Fig. 1 and show that both agree pretty well when β, β 0 ≥ 5.9 and that deviations appear only if the ratios involve lower values of the gauge-coupling parameter, distancing a simulation set-up from the continuum limit; the lower is β the larger the deviation. Still, if β 0 =6.2, the deviation amounts only a 3.6 % for ratios with β=5.7, in the lower border of the validity range of Ref. [74]'s results, and increases up to a 5.3 % with β=5.6, outside this range.
In what follows, we will apply the ratios of lattice spacings from Tab. II obtained by requiring the scaling of the gluon dressing functions and will thus express all momentum-and mass-dimension quantities in units of 1/a(5.8). A conversion to usual physical units can be done by implementing an absolute calibration for one of the lattices at a particular β and applying next the ratios from Tab. II. When needed, we will take r 0 Λ MS =0.586 and ln (a(β)/r 0 ) from Eq.(2.6) in Ref. [74], Λ MS =0.224 GeV from Ref. [63] and the ratio a(β)/a(5.8) from Tab.II to obtain 1/a(5.8) = 1.372 GeV and make the conversion to physical units (this particular result is obtained for β = 6.2 but, as it is apparent from Fig. 1, results derived from any other β between 6.0 and 6.2 only differ by less than a 0.5 %). Now, we are in good position of placing Eq. (2.9) for the lattice Taylor coupling under scrutiny.

B. The Taylor coupling
In order to need no further assumption concerning the nonperturbative running for the continuum α T (k 2 ), we proceed by analyzing ratios of α L T (k 2 ; a) estimated from different lattice set-up's. Indeed, according to Eq. (2.9), these ratios differ from 1 by where the overlined quantities denote that they are expressed in units of 1/a(β 0 ); i.e., # ≡ # a(β 0 ). Therefore, the deviations from 1 for the ratios of α L T (k 2 ; a) appear to be described by an expression, Eq. (3.8)'s second line, not depending on the lattice parameter's set-up, weighted by the factor a 2 (β)/a 2 (β 0 ) − 1. Specially, according to our conjecture in Eq. (2.5), the parameters m g or Λ should presumably result from a nonperturbative resummmation up to all orders in g R , this is the physical interaction and they will thus be universal (except for their possible borrowing of higher-order o(a 2 )corrections by practical fitting); while c α and d α rely on the discretization and will generally depend on the details of the lattice action, the gauge-field definition or the gauge-fixing. However, once these details are fixed and remain the same for all the simulations, as we did, the results obtained for different choices of the bare gauge coupling, g 2 0 (a) = 6/β (see Tab. I), can only differ by the effect of the factor a 2 (β)/a 2 (β 0 ) − 1. This is a main feature that is made strikingly apparent by Fig. 2. To produce it, we first compute α L T (k 2 ; a) after Eq. (2.7) for all the set-up's indicated in Tab. I. In particular, we did it for β = 5.6, 5.7, 5.8, 6.02 and the larger volume of 5.9, and evaluate next Eq. (3.8)'s l.h.s. with β 0 = 5.8, over the momentum intervals where the data for β and β 0 overlap. In order to compute the ratios, the lattice running couplings for β 0 have been estimated at the same momenta as those for each β by performing an interpolation with a Legendre polynomial and propagating errors. Thus, the data displayed in Fig. 2 have been obtained directly from the bare gluon and ghost Green's functions without any renormalization other than the multiplication by 6/β, which introduces no additional freedom for data rescaling. As it can be clearly seen, (i) Eq. (3.8) fits well the data (obtained for five different values of β), explaining satisfactorily their structure and dispersion in terms of β, only controlled by the factor a 2 (β)/a 2 (β 0 ) − 1; (ii) data clearly deviates from a linear behavior on k 2 a 2 (β 0 ), strongly supporting the introduction of the Eq. (2.6)'s nonperturbative logarithmic term, and (iii) are consistent with the emergence of a mass scale preventing from the zero-momentum logarithmic divergence.
For the fit, we have taken m g from Ref. [46] and are so left with the mass-dimension Λ and the dimensionless c α mg Λ −cα −dα 0.331 0.443 0.013 0.237  Fig. 2, obtained for five different ensemble of lattice data with β = 5.6, 5.7, 5.8, 5.9 and 6.02, with mg = 0.455 GeV taken from Ref. [46] and expressed in units of 1/a(5.8) = 1.372 GeV (mg). and d α as free parameters that take the best-fit values gathered in Tab. III.
It is worthwhile to highlight a striking feature shown by Fig. 2 (specially in the lower panel): the data for all the different β happen to cross zero fairly at the same momentum, thus implying that α L T (ζ 2 0 ; a(β)) = α T (ζ 2 0 , a(β 0 )), where ζ 0 is the same physical momentum for all β. In other words, the O(a 2 )-corrections for α L T (k 2 ; a) in Eq. (2.9) become quenched at the same nonzero physical momentum, for which the physical running coupling α T (k 2 ) is recovered, irrespective of the values of the gauge coupling and lattice spacing that are used for the lattice set-up. This momentum can be estimated as from Eq. (3.8) and Tab. III. Actually, Fig. 2 tells us that the Taylor coupling directly computed from the bare gluon and ghost Green's functions obtained from four different lattice set-up's (β ranging from 5.6 to 5.9) only coincides with each other very much near k = ζ 0 . A fifth simulation at β=6.02 is performed in a lattice volume so small in physical units that it cannot reach the zero-crossing low-momentum region without being significantly affected by volume artifacts (as will be discussed below). Its data obtained for k > 0.6 appear however to behave well according to Eq. (3.8) with the parameters displayed in Tab. III.

C. Continuum limit and finite volume effects
After the painstaking scrutiny we have made in the previous subsection for a reliable description of discretization lattice artifacts from the Taylor coupling, Eq. (2.9) can be applied to take the continuum limit and thus obtain the nonperturbative physical running of the coupling, α T (k 2 ), from its corresponding lattice estimate, α L T (k 2 ). Still, for so to do, one needs to make sure that finite volume effects are under control and do not appear entangled with discretization artifacts in the lattice estimates of the lattice coupling. As it happens in the analysis of an amplitude in spectroscopy, corrections of the order of exp (−mL) are expected for any lattice correlation function, where L is the physical size of the hypercubic lattice and m the mass of the physical bound states which propagates all over the lattice. However, in a quenched theory, even the dominant contribution is negligible when the lattice size is of a few fm's, as it should come from the lightest glueball state, for which mL ∼ O(10). On the other hand, when computing the gauge-field correlations functions wich take part in Eq. (2.7), a sizeable effect should appear when the associated gauge-field wavelength is, at least, of the same order as the lattice size, i.e. when Lp < ∼ 1. Thus, the larger is the momentum for which the correlation function is evaluated and the shorter the associated gaugefield wavelength, the less impact the volume effects have. In practice, this impact can be estimated in Fig. 3, where we display the results for the continuum α T (k 2 ), obtained with Eq. (2.9) and the parameters of Tab. III, for the lattice estimates from three simulations made at β=5.9 in three different lattice volumes, V =3.48 4 (L/a=30), 5.56 4 (L/a=48) and 7.42 4 fm 4 (L/a=64) and a fourth one at β=6.02 and V =(3.35 fm) 4 (see Tab. I). We can there clearly appreciate that (i) the two simulations made in lattices of L=5.56 fm and L=7.42 fm at β=5.9 exhibit results plainly compatible within the errors for all comparable momenta, while the one for the same β and L=3.48 fm shows a significant volume effect only for k < ∼ 0.6; (ii) the results from the two simulations made in the same smaller physical volume and different β appear also to fully agree over the whole range of momenta.
We can thus conclude that the finite-volume artifacts affecting the Taylor coupling via gauge-field correlation functions are controlled by the lattice physical volume; and that, in quenched QCD, they are negligible in practice when this physical volume is above (6 fm) 4 . The latter agrees with the findings of Ref. [26]. Therefore, the four simulations at β=5.6,5.7,5.8 and 5.9 we dealt with in Subsec. III B, made in lattices of more than 6 fm, can be taken, in very good approximation, as free of finite-volume artifacts. In Fig. 4, it is shown how the results from these simulations look like, before (upper panel) and after (lower) the continuum extrapolation made through Eq. (2.9). For the sake of comparison, the results obtained at β=6.02 and in a volume V =(3.35 fm) 4 lattice appear also displayed and, after extrapolation, all data for k < ∼ 0.6 dropped. The scaling found for five different simulations, with five different values of β, is extremely good, all data lying strikingly on top of each other after applying Eq. (2.9) with the gluon mass borrowed from literature and the three other fitted parameters shown in Tab. III.

D. Gluon and ghost propagators
To obtain the results displayed in the previous subsection for the Taylor running coupling, one essentially needs to deal with bare gauge-field two-point Green's functions and the bare coupling, g 2 (a) = 6/β, properly combined as Eq. (2.7) reads; and the relative lattice calibration described in Subsec. III A, which left us with the ratios of lattice spacings of Tab. II. We thus apply the H(4)-extrapolation to cure from O(4)-breaking artifacts and express all the Green's functions, so obtained from different lattice set-up's, in terms of the momentum expressed in the same non-standard but physical units, 1/a(5.8); and produce then the plots of Fig. 2. They make strikingly apparent that the remaining cut-off corrections behave as Eq. (2.9) dictates. We have employed a relative lattice calibration which agrees well with that in Ref. [74] but introduces marginal deviations at low β of as much as about 3.6 %. This calibration is anyhow grounded on the assumption that O(4)-invariant artifacts have a negligible impact on the gluon propagator, implying thereupon that its dressing function should exhibit a physical scaling after renormalization and H(4)-extrapolation. The latter is a factual statement 1 discussed in Ref. [68] and, precisely, applied therein to refine the lattice scale setting. Here, it merely implies that, after MOM renormalization and H(4)-extrapolation, both gluon and ghost propagators should behave as dictated by Eqs. (2.5)-(2.6) with (3.10) c α and d α given in Tab. III. This is, strikingly again, confirmed by Figs. 5 and 6. In them, and owing to a sampling of gauge-field configurations of O(1000), we displayed data for the renormalized gluon and ghost propagators with statistical errors of the order of, respectively, one and five per mil. Even at this impressive level of statistical accuracy, the ratios of gluon dressing functions obtained in large-volume lattices at four different β's do not differ from 1 within the errors, except for deeply low momenta (see upper panel of Fig. 5). There, for k < 0.2, results at β=5.6 and V 1/4 =11.3 fm deviate by around a 2 % from those at β=5.7,5.8 and 5.9 and V 1/4 7 fm which, on their side, remain compatible with each other. This strongly suggest that this slight deviation is caused by a still-remaining finite-volume effect. This systematic effect is very nearly negligible, made only apparent by the huge statistics here employed, and only happens at very low momenta: namely, there is no impact from it at k 0.29 (highlighted by a red dashed line in Figs. 5 and 6), the momentum for which all the lattice estimates for the Taylor coupling coincide in Fig. 2. The same is shown in the lower panel of Fig. 5, where the gluon propagators from the four simulations appear plotted and lie, very accurately, on top of each other. This is a nonobvious result firmly confirming the starting hypothesis of the relative scale setting and the efficiency of the H4extrapolation.
In Fig. 6, the ratios of ghost propagator dressing functions from the four lattice simulations are clearly proven to behave according to [Upper panel] The ratios of lattice gluon propagator dressing functions after MOM renormalization, DL(k 2 ; a(β))/DL(k 2 ; a(β0)), with β0=5.8 and β=5.6,5.7 and 5.9, are shown to be compatible with 1 for all momenta except when k 2 a 2 (5.8) < ∼ 0.05, where nearly negligible volume artifacts takes place. [Lower panel] The gluon propagator, ∆(k 2 ) = DR(k 2 )/k 2 , for β=5.6,5.7,5.8 and 5.9 exhibit a nearly perfect physical scaling. Here, we plotted in terms of k in units of 1/a(5.8), instead of the squared momentum, to make more apparent the domain of deeply low momenta. In both cases, the red dashed line is placed at ka(5.8) 0.29, to indicate where α L T (k 2 ; a) from all the different simulations crossed.
nearly perfect physical scaling (lower panel). The renormalization point is chosen here to be ζ=0.8, lying thus within the momentum range, and far from its borders, for the four simulations. This is why, precisely, we do not include here the data from the simulation at β=6.02 in the small volume (V 1/4 = 3.35 fm): they are significantly affected by finite-volume artifacts at this renormalization point and, therefore, the MOM renormalization prescription will contaminate with these artifacts the whole momentum range. Its momentum range and β=5.8 simulation's offer anyhow a reliable overlap which makes possible the determination of the ratio a(6.02)/a(5.8) in Tab. II.
Furthermore, as explained in Fig. 1's  FIG. 6. The ratios of lattice ghost propagator dressing functions after MOM renormalization, FL(k 2 ; a(β))/FL(k 2 ; a(β0)), with β0=5.8 and β=5.6,5.7 and 5.9, are shown differ by 1 according to Eq. (3.11) (solid lines). [Lower panel] After extracting the ghost dressing function FR(k 2 ) from simulations at β=5.6,5.7,5.8 and 5.9 with Eqs. (2.5)-(2.6), a nearly perfect physical scaling is exhibited. Here, we plotted in terms of k in units of 1/a(5.8), instead of the squared momentum, to make more apparent the region of deeply low momenta. In both cases, the red dashed line is placed at ka(5.8) 0.29, to indicate where α L T (k 2 ; a) from all the different simulations crossed.
ploited the overlap with β=5.8 and so extracted the ratio r = a(6.202)/a(5.8) in Tab. II. Then, we calculated a(β)/a(6.202) = r −1 a(β)/a(5.8) and plotted the results in Fig. 1. Indeed, as can be seen in the plot, our estimates for β > 5.9 appear to be in nearly perfect agreement with those of Ref. [74]. It should be reminded that no fit is made here to produce the results displayed in Figs. 5 and 6. After the scale setting, we merely apply a MOM renormalization prescription, compute the ratios of dressing functions and, when needed, use Eq. (3.11) with the parameters of Tab. III, obtained for the Taylor coupling, and Eq. (3.10). It is important to highlight that the connection between the lattice Taylor coupling and the dressing functions relies on the field Theory and its renormalization: Eq. (2.7) for α L T (k 2 ; a) involves the bare two-point Green functions and the bare gauge coupling which is directly given by the parameter β in the lattice gauge action. Their dependence on the lattice spacing is singular for each but cancels in Eq. (2.7), remaining only a residual one which vanishes in the continuum limit. Such a residual dependence is also related to the one from the MOM renormalized Green's functions. On the other hand, the way in which the lattice spacing relates to β depends on the lattice action and determines the physical scaling of quantities from lattice simulations; namely, the dressing functions in the continuum limit. How all this takes place and match is highly non-trivial. This is what we have exposed here.

IV. SUMMARY AND CONCLUSIONS
We have carefully examined the physical scaling violations of two-points gluon and ghost Green's functions, when they are obtained from a fixed cut-off simulation in lattice regularization, after MOM renormalization and before extrapolation to the continuum limit. Specially, we performed a combined analysis of the gaugefield propagators and the Taylor coupling (obtained from them) seeking a consistent description of results from many lattice simulations, with different lattice spacings ranging widely from 0.07 to 0.24 fm. It should be highlighted that, when the physical scale is properly set, the H(4)-extrapolation cures efficaciously the gluon propagator from cut-off deviations up to the order O(a 2 ) and, after renormalization, the results thus obtained show a very striking physical scaling. This is not the case either for the ghost propagator or the Taylor coupling, which are affected by sizeable O(4)-breaking artifacts. However, we can accurately deal with these artifacts and get thus an insightful understanding of the impact of discretization cut-off effects on the two-point Green's functions, which makes therefore possible and reliable their very precise continuum extrapolation. This will be of very much help in any future work aiming at extracting QCD parameters from these lattice Green's functions, or just at their comparison when they are obtained from lattice set-up's with different discretization spacings.
It is furthermore remarkable that the violations of the physical scaling herein scrutinized behave, within our approximation order and after H(4)-extrapolation, as the squared lattice spacing times a function of the physical momentum saturated by the gluon mass in the IR limit. The latter appears to suggest that the emergence of the gluon mass become also manifest in the nonperturbative cut-off effects, within a renormalization prescription, before removing them by taking the appropriate limit.