Three-gluon running coupling from lattice QCD at $N_f=2+1+1$: a consistency check of the OPE approach

We present a lattice calculation of the renormalized running coupling constant in symmetric (MOM) and asymmetric ($\widetilde{\rm MOM}$) momentum substraction schemes including $u$, $d$, $s$ and $c$ quarks in the sea. An Operator Product Expansion dominated by the dimension-two $\langle A^2\rangle$ condensate is used to fit the running of the coupling. We argue that the agreement in the predicted $\langle A^2\rangle$ condensate for both schemes is a strong support for the validity of the OPE approach and the effect of this non-gauge invariant condensate over the running of the strong coupling.

cook out as many different renormalization schemes as there are possible kinematical configurations. As the signal for a three-point correlation function, suffering from stronger statistical fluctuations, is much harder to be extracted from lattice simulations than the one for two-point functions, the precision so attained is not comparable with the one achieved when using the T-scheme coupling. The interest of computing α s in different schemes is therefore not to obtain a precise value of Λ QCD but, rather, to test the OPE framework and to gain thus some insight into the nature of the nonperturbative corrections.
In this paper we present the lattice evaluation of the MOM QCD coupling defined through the three gluon vertex for two different kinematical configurations: the symmetric (three equal momenta) and asymmetric (one vanishing momentum) ones. The high-statistics ensemble (800 configurations) of lattice gauge fields we exploit takes into account the dynamical generation of up, down, strange and charm quarks (N f = 2 + 1 + 1). This leaves us with two main "aces" for our game: (i) the Wilson coefficients for the leading contribution in the OPE of the two couplings, as will be seen, differ very much from each other; and (ii) the perturbative running is very reliably known as the same N f = 2 + 1 + 1 lattice configurations provides, via the T-scheme coupling determination, with an accurate estimate of Λ QCD [32], compatible with PDG world average [1], that can be used here. We put ourselves in a near unbeatable position to check the OPE framework, as the nonperturbative contributions supplementing the perturbative running to account for the lattice data of both couplings can be properly isolated and compare to each other. One can see then if they differ as much as OPE predicts.
The structure of the paper is as follows: in section II the renormalization schemes and lattice setup used are described. In section III a reminder of the OPE results has been included. Finally in section IV our main results are presented and we concluded in section V.

II. RENORMALIZATION SCHEMES.
The starting point for this calculation shall be the gauge configurations produced by the European Twisted Mass collaboration (ETMC) for 2 + 1 + 1 dynamical quark flavors that provide a realistic description of the QCD dynamics including heavy flavours. These gauge field configurations, after fixing Landau gauge, allow to compute the renormalized running coupling in momentum substraction schemes. In particular we will focus on the coupling defined from three-gluon vertices.
The three-gluon vertex ( Fig. 1) can be computed from the lattice for any momenta p 1 , p 2 and p 3 satisfying p 1 + p 2 + p 3 = 0. In particular, we will concentrate on the symmetric three gluon vertex (p 2 1 = p 2 2 = p 2 3 ) and the asymmetric one (p 3 = 0 and therefore p 2 = −p 1 ). The renormalized coupling can be straightrowardly defined from gluon propagators and vertices (a detailed description of the procedure can be found in [3]). The renormalized coupling is defined by: where µ 2 is the renormalization scale, to be fixed for each renormalization scheme, G 2 (p 2 ) is the bare gluon propagator extracted from the lattice: is the gluon field renormalization constant, and G (3) (p 2 1 , p 2 2 , p 2 3 ) is the scalar function extracted from three gluon vertex This scalar function is defined as the coefficient of the tree level tensor and is obtained after projecting the vertex onto the adequate tensor as described in [3].
This procedure allows to compute the running coupling α(µ 2 ) = in momentum substraction schemes both from the symmetric three gluon vertex (MOM) and from the asymmetric one ( MOM).
The symmetric vertex requires the three momenta p 1 , p 2 and p 3 to satisfy the constrain p 1 + p 2 + p 3 = 0 simultaneously with p 2 1 = p 2 2 = p 2 3 which is rather rare in the lattice. It means that there are rather few momenta where the vertex can be evaluated. For the asymmetric one, the constrain is less restrictive and the vertex can be evaluated at any lattice momenta p.
As mentioned above, the two and three-point gluon Green functions will be computed from lattice gauge field configurations simulated at N f =2+1+1 by the ETM collaboration [36,37]. The details of the computation can be of the background field method and the pinching technique [34,35]. found in [14] and references therein. We have exploited here a set of 800 configurations for a lattice volume 48 3 x96 and β = 2.10, where the lattice parameters κ c = 01563570, µ l = 0.002, µ σ = 0.120, µ c = 0.385, are fixed so that light quark mass is set to ∼ 20MeV and the strange and charm are set to ∼ 95MeV and 1.5GeV respectively (see [14] and references therfein for the details of the simulations). According to Ref. [32], the lattice spacing corresponding to this set-up is a(2.10) = 0.0583 (11)fm.
Contrarily to the continuum ones, the lattice scalar functions do not depend only on the momentum squared p 2 , due to the lattice discretization. Indeed, the lattice discretization breaks the O(4)-symmetry introducing an spurious dependence on the invariants of the group H(4). These O(4)-breaking lattice artefacts can be efficiently removed by using the so-called H(4)-extrapolation procedure [38][39][40]. This method works efficiently for gluon propagator and asymmetric vertex, where there is a high number of momenta at which the Green function can be evaluated. The correction of lattice artifacts for the symmetric vertex is not so efficient due to the lower number of lattice momenta and the fact that it depends on three momenta instead of one. This introduces a limitation for the larger momenta that can be used in the schema MOM that, in practice, implies a smaller fitting window than in MOM scheme. As will be seen in next section, we take momenta below a(2.10)p ≃ 1.6 (around 5.5 GeV, in physical units) for the fit in the MOM case, while the fitting window is restricted only to momenta below a(2.10)p ≃ 1.3 (around 4.5 GeV) for the fit in MOM, to avoid the noise induced by the non-properly-cured lattice artefacts.

III. OPE NONPERTURBATIVE PREDICTIONS.
The running of the strong coupling constant with momentum, obtained from QCD perturbation theory corrected by a nonperturbative leading OPE power contribution, can rather generally read [12,14] where the subindex R specifies any particular renormalization scheme and α pert gives the running behaviour perturbatively obtained from the integration of the QCD beta function at that R scheme, with h R = α R (µ 2 )/(4π) and β 0 = 11 − 2/3N f , β 1 = 102 − 38/3N f , being scheme-independent coefficients. The result for α pert from the integration of Eq. (4) and its conventional perturbative inversion, in terms of momenta and the QCD scale Λ R , can be found in [1]. Within the bracket, c R is given by the tree-level Wilson coefficient contribution, γ A 2 0 is the first coefficient for the local operator A 2 anomalous dimension, determining the Wilson-coefficient leadinglogarithm contribution, which is found to be scheme-independent; and R(α, α 0 ) encodes the higher-order logarithmic corrections for the leading Wilson coefficient [16]. In refs. [12,15,32], Eq. (3) particularized to the MOM T-scheme accounted very accurately for the running of the corresponding lattice data with momenta. This allowed for a precise determination of Λ T , and hence Λ MS pretty in agreement with the PDG [1] "world average".
Hereupon, we will mainly concentrate on the running coupling renormalized in both MOM and MOM schemes, for which the three-loop beta coefficients are known 2 . For N f = 4, one is left with and with the following ratios that can be exactly obtained from the first non-trivial coefficient for the expansion of MOM and MOM couplings in terms of the MS one. The tree-level Wilson coefficients for the strong coupling in both schemes have been also studied [7,8] and their computation gives It is worthwhile to recall that, in the MOM case, as a consequence of the soft gluon field in the three-gluon Green function defining the vertex, Eq. (3) only results after the factorization of a leading higher-dimension condensate 3 in the OPE, induced by a vacuum insertion approximation [8] (the same vacuum insertion approximation have been proven to work for the OPE expansion of the ghost-ghost-gluon Green function [41]). The function R(α, α 0 ) is related to the A 2 anomalous dimension but also depends on the scheme we used to define the coupling [16]. For the MOM T-scheme [11], as the coupling can be directly related to gluon and ghost propagators involving no three-point Green function, it has been computed at the O(α 4 )-order [12,16]. Its computation is nevertheless cumbersome when dealing with three-gluon Green function is needed. Thus, in the following, we will take R(α, α 0 ) = 1 and will work at the leading-logarithm approximation.

IV. RESULTS.
Then, Eq. (3) can be fitted to the lattice data, with g 2 A 2 as a free parameter, for both MOM and MOM couplings, with their perturbative predictions obtained by the integration of the beta function with the coefficients given in Eq. (6) and the ratios of Λ's in Eq. (7). We will take Λ MS = 314 MeV, as an input from ref. [32] 4 where, as above mentioned, the MOM T-scheme coupling is computed from the lattice and confronted with Eq. (3), properly particularized (ref. [32] upgrades the previous results of refs. [14,15]).
Figs. 2 shows the lattice results for the running coupling constant and the best fits with Eq. (3) and g 2 A 2 from Tab. I, in the MOM and MOM schemes. Within the framework of OPE and SVZ sum-rules approach, the gluon condensate needs to take similar values in the OPE for the two couplings. The agreement in the values of the extracted condensates is therefore a strong indication of the validity of this approach, at least for a window of momenta not lying in the deep IR domain. As the nature of the OPE condensates is the object of a recent controversy [43,44], It is worth to point out that our check is validating the sum-rules factorization but does not tell necessarily anything about the nature of the condensates.

MOM
MOM T-scheme g 2 A 2 (GeV 2 ) 5.5 ± 0.8 6.0 ± 1.  R,µ 0 renormalized at µ0 = 10GeV extracted from the best fit of the lattice running couplings αMOM(µ) and α MOM (µ). For comparison, the MOM T-scheme result, estimated from ref. [32] data as explained in the text, is also included. The values of the condensates obtained in this paper are sizably larger than the ones reported in [14,15,32] for the Taylor scheme. It is however important to note that, in the Taylor scheme analysis of those papers, the beta function is expanded at the four-loop level for its integration and that the Wilson coefficient has been computed up to the O(α 4 )-order. Here, in our analysis of MOM and MOM couplings, we have only used a three-loop beta function and Wilson coefficient at the leading order. Indeed, the effect of including higher orders in either the perturbative part of Eq. (3) or the Wilson coefficient is well known to reduce the value of the condensate [11,45]. Alternatively, for the sake of a consistent comparison, we repeated the analysis of Ref. [32] under the same approximation level (the best which can be here coherently attained) applied for the current one: the perturbative coupling expanded only up to three-loops and the Wilson coefficient kept at the leading-logarithm approximation. One then obtains always at the renormalization point µ = 10 GeV. This last result happens to be exactly the same as the one for MOM and to lie in the same ballpark as that for MOM, as can be seen in Tab. I. Indeed, one can also apply both Λ MS from [32] and g 2 A 2 from Eq. (9) to Eq. (3) and, without free parameters to be fitted, account for the lattice data for MOM and MOM couplings with χ 2 /d.o.f. that would be then 1.37 for the former and 0.76 for the latter. It is worthwhile to emphasize that refs. [14,15,32] exploited the same (800) lattice configurations for the gauge fields at a bare coupling, β = 2.1, here analysed, but also configurations at β = 1.90 for three different light-quark twisted masses (500 each) and at β = 1.95 (150). The latter gives us the grounded conviction that lattice artefacts are properly under control in obtaining Λ MS and g 2 A 2 from the Taylor coupling, as done in Eq. (9). Therefore, that Eq. (3) successfully describes the MOM and MOM coupling data for a given momentum window, with the same Λ MS and g 2 A 2 , strongly indicates that lattice artefacts appear to be negligible also for them, after H(4)-extrapolation, within such a window.
That the condensates obtained from both MOM and MOM takes the same value is a demanding result, as the Wilson coefficient is three times larger for the former than for the latter. This implies that deviations from the perturbative behaviour should be very different in both cases, being consistent with the ratio of 3 given by Eq. (8). This can be seen in Fig. 3.a and, otherwise presented, as follows: Eq. (3), in the leading logarithm approximation, and Eq. (8) left us with: which provides with a very demanding consistency check for the OPE and SVZ sum-rules approach, which is totally equivalent to the compatibility of condensates in Tab. I. To perform this check, we need to compute Eq. (10)'s l.h.s. from the lattice data for the MOM and MOM three-gluon coupling and their perturbative predictions obtained again by the integration of the beta function with the coefficients given in Eq. (6), the ratios of Λ's in Eq. (7) and Λ MS = 314 MeV from ref. [32]. However, as the lattice momenta for MOM and MOM differ, and aiming at employing as large a statistics as possible, Eq. (10)'s l.h.s. will be indirectly computed by fitting both numerator and denominator, within as large as possible a momentum domain for each, to as suggested by Eq. (3), with a R as the only free parameter to be fitted. This can be seen in the plots of Fig. 3.b, where it clearly appears that, as the running given by Eq. (11)'s r.h.s. is near the same for both schemes, a MOM is to be rather larger than a MOM . Thus, once both parameters are fitted, they can be applied to compute Eq. (10)'s l.h.s., that compares remarkably well with Eq. (10)'s r.h.s. evaluated through the OPE results given by Eq. (8). It should be noticed that, in Eq. (12), α MOM /α MOM = 1 + O(α) has been applied, as corresponds to our approximation level (see Eq. (10)). We performed the fit in a momentum window p ∈ (2.8, 4.5)GeV in the MOM case and p ∈ (2.8, 5.5)GeV in MOM. The errors for the fitted parameters, a R , have been computed by applying the jackknife procedure, and propagated then into the final result for the ratio.

V. SUMMARY AND CONCLUSIONS
We used lattice gauge field configurations, generated with four twisted-mass dynamical quark flavours (two light degenerate and two heavy) within the framework of ETM collaboration, to compute the running of the QCD coupling constant, α s , defined from the symmetric and asymmetric three-gluon vertices. This leads us to two different MOMtype renormalization schemes for the coupling, where their running with momenta, roughly from 3 to 6 GeV, has been described only after supplementing the well-known perturbative prediction with a non-perturbative correction dominated by a non-vanishing dimension-two gluon condensate in Landau gauge, A 2 .
As can be seen in Fig. 3.a, the perturbative estimate for the coupling in MOM scheme is larger than the one in MOM. This is a direct consequence from the way their three-loop β coefficients in Eq. (6) and the ratios in Eq. (7) compare to each other. Contrarily, the OPE nonperturbative leading correction for MOM is predicted to be three times larger than for MOM, irrespectively with the value of A 2 . Being so that altogether, for a sufficiently large condensate, the relative strength between the two nonperturbative couplings will result reversed with respect to the perturbative case. We indeed found the lattice data for the MOM and MOM coupling to follow this reversed pattern. We also measured the ratio between their nonperturbative corrections and found it to be 2.7 (8), strongly supporting the OPE approach to account for the nonperturbative contributions. It is worth to recall that the prediction for this ratio is only relying on the Shifman-Vainshtein-Zakharov technology [46,47] to compute the OPE Wilson coefficients and the universality of the involved condensate.
Finally, the value of the condensate is also found to be consistent with that resulting from the analysis of the running of the coupling in the MOM T-scheme, in ref. [32], after restricting the calculation there to the same order than here. This provides thus with additional support for the very accurate estimate of Λ MS obtained therein.