Hyperasymptotic approximation to the plaquette and determination of the gluon condensate

We give the hyperasymptotic expansion of the plaquette with a precision that includes the terminant associated to the leading renormalon. Subleading effects are also considered. The perturbative series is regulated using the principal value prescription for its Borel integral. We use this analysis to give a determination of the gluon condensate in SU(3) pure gluodynamics that is independent of the scale and renormalization scheme used for the coupling constant: $\langle G^2 \rangle_{\rm PV} (n_f=0)=3.15(18)\; r_0^{-4}$.


Introduction
The expectation value of the plaquette calculated in Monte Carlo (MC) simulations in lattice regularization with the standard Wilson gauge action [1] reads where Λ E is a Euclidean spacetime lattice and P x,µν = 1 − 1 6 Tr U x,µν + U † x,µν .
(1.2) For details on the notation see Ref. [2]. This quantity can be computed using the operator product expansion where a denotes the lattice spacing, and is the so-called non-perturbative gluon condensate [3], which, under some conditions, it is expected to scale like Λ 4 where (in an arbitrary scheme) (1. 6) and so on.
The Wilson coefficient multiplying the gluon condensate is proportional to the inverse of the beta function of α [4,5]. This fixes the Wilson coefficient exactly: Note that C G (α) is scheme-dependent not only through α, but also explicitly, due to its dependence on the higher β-function coefficients: β 2 , etc.. The c k depend on the β i with i ≤ k + 1 via Eq. (1.7). For j ≤ 3 the coefficients β j are known in the Wilson action lattice scheme. β latt 2 has been computed diagrammatically [6][7][8]. The value for β latt 3 that we use [9] is an update of [10], and was obtained by calculating the normalization of the leading renormalon of the pole mass, and then assuming the corresponding MS-scheme expansion to follow its asymptotic behaviour from orders α 4 s onwards. Similar estimates, β latt 3 ≈ −1.37 × 10 6 up to β latt 3 ≈ −1.55 × 10 6 , were found in Ref. [11] using a very different method. For convenience, we also write the expansion coefficients c k defined in Eq. (1.7) in terms of the constants that appear in Eq. (1.5): (1.8) The perturbative sum and the leading nonperturbative correction in Eq. (1.3) are illdefined. The reason is that the perturbative series is divergent due to renormalons [12] (for a review see [13]) and other, subleading, instabilities. This makes any determination of G 2 ambiguous, unless we define how to truncate or how to approximate the perturbative series. Any reasonable definition consistent with G 2 ∼ Λ 4 can only be given if the asymptotic behaviour of the perturbative series is under control. This has only been achieved recently [2], where the perturbative expansion of the plaquette was computed up to O(α 35 ). The observed asymptotic behaviour was in full compliance with renormalon expectations, with successive contributions starting to diverge for orders around α 27 -α 30 within the range of couplings α typically employed in present-day lattice simulations.
Extracting the gluon condensate from the average plaquette was pioneered in Refs. [14][15][16][17], and many attempts followed during the next decades, see, e.g., Refs. [18][19][20][21][22][23][24][25][26][27]. Nevertheless, they suffered from insufficiently high perturbative orders and, in some cases, also finite volume effects. The failure to make a controlled contact to the asymptotic regime prevented a reliable lattice determination of G 2 , where one could quantitatively assess the error associated to these determinations. This problem was first solved in [28]. In such paper, for the first time, the perturbative sum was computed with superasymptotic accuracy for the case of 4 dimensional SU(3) gluodynamics. This allowed to obtain a reliable determination of G 2 that scaled as Λ 4 . One issue raised was to determine to which extent such a result was independent of the scheme used for the coupling constant. The answer to this question can be given within the general framework of hyperasymptotic expansions of renormalizable quantum field theories, as developed in [29][30][31][32] (see [33,34] for the original works in the context of ordinary differential equations). This answer was indeed given in [29], where it was concluded that the error of using the superasymptotic approximation to the perturbative sum was of O( α(1/a)Z P Λ 4 ), where Z P is the normalization of the leading renormalon. This error then sets the parametric precision of the determination of the gluon condensate using the superasymptotic approximation. Note that the scheme dependence of Z P and Λ 4 cancels each other. Therefore, the only remaining/leading scheme/scale dependence of the error is due to the α(1/a) prefactor. It is the purpose of this paper to revisit such analysis and to encode such a result in the more general framework of the hyperasymptotic expansions. In particular, we will reach a better (or more robust) accuracy by incorporating the leading terminant. We will also discuss subleading effects. We confirm that the result we obtain is independent of the scheme/scale used for the renormalization of coupling constant (up to terms that are higher order than the accuracy reached by the hyperasymptotic approximation).
In order to carry out the above program, the first step is to regularize the perturbative sum, which we do using the Principal Value (PV) prescription. Only after regularizating the perturbative sum, the definition of the gluon condensate is unambiguous and the operator product expansion of the plaquette reads This expression is, in practice, formal, as the exact expression of S PV is not known. This would require the exact knowledge of the Borel transform of the perturbative sum. Nevertheless, it is possible to obtain an approximate expression of it with a known parametric control of the error using its hyperasymptotic expansion. The accuracy of this expansion is limited from the information we get from perturbation theory. For the case at hand, we have 10) where N is the maximal order in perturbation theory that is included in the perturbative expansion. Within the hyperasymptotic counting, approximating S PV by S P ≡ N P n=0 p n α n+1 , the first term in Eq. (1.10), corresponds to the superasymptotic approximation, which we label as (0, N P ). Adding Ω G 2 to the superasymptotic approximation corresponds to (4,0) precision in the hyperasymptotic approximation and adding the last term corresponds to (4, N ) precision 1 .
In Eq. (1.10), we take 11) as the order at which we truncate the perturbative expansion to reach the superasymptotic 1 The labeling (D,N) in general is defined in Refs. [30,31].
approximation. By default, we will take the smallest positive value of c that yields an integer value for N P , but we also explore the dependence of the result on c. Note that the value of N P that we use here is slightly different from the value used in [28] to truncate the perturbative expansion with superasymptotic accuracy. In that reference, such number was named n 0 and was determined numerically. We will ellaborate on this difference later. The asymptotic expression of the coefficients of the perturbative expansion was worked out in [2]. We repeat it here for convenience Note that the parameters b 1 and b 2 : that describe the leading pre-asymptotic corrections depend on the expansion coefficients c 0 and c 1 , defined in Eq. (1.7), of the Wilson coefficient of the gluon condensate. Ω G 2 is the terminant associated to the leading renormalon of the plaquette. It can be easily obtained from the general formulas given in [29]. It reads where we take ∆Ω(db) from Eq. (33) of [29] taking γ = 0, and Ω G 2 can also be written in the following way 16) or where , where η c ≡ −4b + 8π β 0 c − 1. The value of Z P was determined approximately (for n f = 0) in [2]: This is the value we will use in this paper. Actually, its error will give the major source of uncertainty in the determination of Ω G 2 , of the order of 40%. The other source of error is due to the fact that only approximate expressions are available for Ω G 2 (see Eq. (1.15), Eq. (1.16), and Eq. (1.17)), as we do not know the complete set of coefficients of the beta function in the lattice scheme. Nevertheless, we can study the convergence pattern of the weak-coupling expansion. We show the results in Table 1 for a representative set of values of α in the interval that we will use later. The first observation is that we observe a very good convergent pattern of the weak-coupling expansion of the terminant using Eq. (1.17) or Eq. (1.15), consecutive terms quickly become smaller. The latter is obtained using the exact numerical determinations of ∆Ω. The second observation is that the strict weak-coupling expansion used in Eq. (1.17) saturates perfectly the exact numerical determination of Eq. (1.15) for analogous precision. On the other hand, if we want to use Eq. (1.16) the convergence is not good. We have to go to β-values (β ≡ 3/(2πα)) rather larger than 6 to get decent accuracy. What lies behind is the fact that Λ latt is not well approximated by its weak coupling expansion at low orders (see the discussion in [28] and [30]). This bad convergence of the weak coupling expansion of Λ latt has to be compensated by a bad convergence of the O(α) corrections of the weak coupling expansion in Eq. (1.16). A similar behavior, but less severe, was observed in [30]. In that paper, the analogous to the Wilson coefficient C G was 1. Nevertheless, the main point is that now the power of Λ is four, whereas in [30], the power of Λ was one. This makes that the relatively bad convergent behavior observed in [30] gets amplified by a factor of four here. 2 Therefore, in the following, we will always use Eq. (1.17) as our approximated expression for Ω G 2 , as it produces a nicely convergent series with a controlled scheme dependence, as the weak coupling expansion is organized in terms of a single parameter: α. The error associated to truncating the expansion in Eq. (1.17) is estimated by observing the convergent pattern of the LO, NLO and NNLO results in Table 1. From LO to NLO, in the worst cases, the differences are close but below 50%, and from NLO to NNLO the differences are below 10%. One could then expect the NNNLO contribution to be at the level of few percent, which can be neglected all together in comparison with the ∼ 40% error associated to Z P .
2 Determination of the gluon condensate Following the notation of [28], we determine the gluon condensate from the following equation: If S PV and P MC were known exactly, this equality is expected to hold up to corrections of O(a 2 Λ 2 ). Nevertheless, neither S PV nor P MC are known exactly. On top of that, we have to account for the fact that C −1 G and the relation between a and β are also known in an approximated way. We now discuss how we determine them and their associated individual errors.
We take the MC data from [36]. Similarly to what was done in [28], in this analysis we restrict ourselves to the more precise N = 32 data and, to keep finite volume effects under control, to β ≤ 6.65. We also limit ourselves to β ≥ 5.8 to avoid large O(a 2 ) corrections. At very large β-values, obtaining meaningful results becomes challenging numerically: the individual errors both of P MC (α) and of S PV (α) somewhat decrease with increasing β. However, there is a very strong cancellation between these two terms, in particular at large β-values, since this difference decreases with a −4 ∼ Λ 4 latt exp(16π 2 β/33) on dimensional grounds, while P MC depends only logarithmically on a. We illustrate this cancellation in Fig. 1. well approximated by its weak-coupling expansion. Consequently, in this scheme, Eq. (1.15), Eq. (1.17), and Eq. (1.16) converge well. It is also interesting to study these equations with n f = 0, in view of future full QCD analyses, which will involve the incorporation of active massless fermions (for preliminary studies see [35]). As expected, one observes a better convergence of the weak-coupling expansion taking n f = 3 than n f = 0.  Figure 1.
The second line is basically indistinguishable with respect to zero with the scale resolution of this plot. The statistical errors are smaller than the size of the points. Equation (1.5) is not accurate enough in the lattice scheme for the β-values used in this paper. Instead, we employ the phenomenological parametrization of Ref. [37] obtained by interpolating non-perturbative lattice simulation results. 3 Equation (2.2) was reported to be valid within an accuracy varying from 0.5% up to 1% in the range [37] 5.7 ≤ β ≤ 6.92, which includes the range β ∈ [5.8, 6.65] we use in this paper. This range corresponds to (a/r 0 ) 4 ∈ [3.1 × 10 −5 , 5.5 × 10 −3 ], and covers more than two orders of magnitude, i.e. in units of energy, we use lattice data in the region 1/a ∼ (3.66 r −1 0 ÷ 13.42 r −1 0 ). For the inverse Wilson coefficient the corrections to C G = 1 are small. However, the O(α 2 ) and O(α 3 ) terms are of similar sizes. We will account for this uncertainty in our error budget.
We now turn to S PV (α). As we have mentioned above, we compute it using the hyperasymptotic expansion. This introduces a parametric error according to the order we truncate this expansion. On top of that, the coefficients p n , obtained in Ref. [2], are not known exactly. They carry statistical errors, and successive orders are correlated. Using the covariance matrix, also obtained in Ref. [2], the statistical error of S P (α) can be calculated. In that reference, coefficients p n (N ) were first computed on finite volumes of N 4 sites and subsequently extrapolated to their infinite volume limits p n . This extrapolation is subject to parametric uncertainties that need to be estimated. We follow Ref. [2] and add the differences between determinations using N ≥ ν points for ν = 9 (the central values) and ν = 7 as systematic errors to our statistical errors. This is the same error analysis as the one used in [28]. We emphasize though, that the order we truncate the perturbative series, N P , is different from the one used in [28] (which was named n 0 in this reference). The difference between both determinations gives an estimate of the parametric error of the determination of S PV (α) by using the superasymptotic approximation S P . The magnitude of Ω G 2 gives an alternative estimate of the error associated to the truncation of the hyperasymptotic approximation. It is also interesting to see the magnitude of changing N P by one unit by fine tunning c from the smallest positive value that yields an integer value of N P to the smallest (in modulus) negative value that yields an integer value of N P . Typically this yields slightly smaller errors. We illustrate this discussion in  [30,31]). If we increase the accuracy of the hyperasymptotic expansion by adding the terminant Ω G 2 to the superasymptotic approximation, the parametric error decreases, and the accuracy reached is (4,0) (note that the statistical error does not change). With this accuracy, the parametric error is ∼ O(e −4 2π β 0 α(1/a) (1+log(3/2)) ) ∼ O((aΛ) 4(1+log(3/2)) ) (see the discussion in [30,31]). Note that 4 log(3/2) 1.6 < 2. Therefore, these effects are parametrically more important than the next nonperturbative power corrections. Compared with the typical size of the terminant Ω G 2 , these effects are suppressed by a factor of order ∼ O((aΛ) 4 log(3/2)) ). In the energy range we do the fits, this yields suppression factors in the range ((aΛ MS ) 4 log(3/2) ) ∈ (0.007, 0.05), where we have taken Λ = Λ MS to be more conservative. This discussion can be affected by powers of α. It is expected that there is an extra suppression factor of α 3/2 (as √ α is already included in the terminants the real suppression factor would be of order α). Depending on the scheme, the size of this extra factor is different. In any case, they go in the direction to make the estimate of the error smaller. We will not dwell further in this discussion of the parametric error of the (4,0) hyperasymptotic accuracy, because we only approximately know Ω G 2 and its error will hide the signal of these O((aΛ) 4(1+log(3/2)) ) effects. For Ω G 2 we use the analytic expression in Eq. (1.17) truncated at O(α 2 ). The error of this expression comes from Z P , and from the truncation of the weak coupling expansion of the terminant. The largest source of error comes from Z P . Due to its size, this error overwhelms the parametric error associated to Ref. [28] (0,N P ) We also show the superasymptotic approximation obtained in [28] truncating at the minimal term determined numerically. The horizontal green band and its central value are our final prediction, and the associated error, for the gluon condensate displayed in Eq. (2.5).
higher-order terms in the hyperasymptotic expansion. Irrespective of the discussion of the error of the (4,0) accuracy, it is nice to see that adding the terminant to the superasymptotic expression makes the jumps that we had with the superasymptotic approximation disappear. Adding the terminant also makes the resulting curve flatter. The dependence in N P (or in other words c) gets much milder too. We illustrate all this in Fig. 2.
In principle, we know perturbation theory to orders high enough to include the last term written in Eq. (1.10) and reach (4, N ) accuracy. Nevertheless, we find that the errors of p n for large n hide the signal. We show in Fig. 3 how the statistical errors grow as we increase N . On the other hand it is rewarding to see that the dependence in c basically vanishes. We elaborate more on the error analysis of the (4, N ) hyperasymptotic approximation in the next section.

Fit
We now implement the discussion of the error of the previous section to Eq. (2.1) for the different truncations of the hyperasymptotic expansion, and for the associated determination of the gluon condensate. The statistical errors of the fits are those of the MC determination of P MC and the statistical error of S PV . The latter is generated by the statistical errors of the coefficients p n . As successive orders are correlated, we use the covariance matrix and, by propagation of the error, compute the statistical error of S PV . This is the same method followed in [28] for the superasymptotic approximation. We then combine the statistical error of P MC and of S PV in quadrature, which is then used to generate the fits. We show the size of these two different statistical errors in Fig. 2.
We now turn to systematic uncertainties. One is the infinite volume extrapolation of the coefficients p n discussed before. Another source of systematic errors is the possible existence of O(a 2 Λ 2 ) corrections to the fit. Looking to Fig. 2, within statistical errors, there is no clear signal of the O(a 2 Λ 2 ) in the whole energy range used. Therefore, our default fits will be in the range β ∈ [5.8, 6.65] and use the difference with fits in the range β ∈ [6, 6.65], as an estimate of these effects. Another source of systematic uncertainties is the incomplete knowledge of C G . We consider the difference between truncating C G to O(α 2 ) and to O(α 3 ) as an estimate of this error. Next, there is a scale error of about 2.5%, translating a into units of r 0 . The other systematic uncertainties are specific to each truncation of the hyperasymptotic approximation used, which we next address. Therefore, we then move to discuss the final error of the different orders in the hyperasymptotic approximation.
*) Precision (0, N ) It is not the purpose of this paper to make a detailed discussion of (how to estimate) the error of the perturbative expansion with precision (0, N ) with N N P . Nevertheless, we can not avoid mentioning that there are scenarios where standard ways to estimate the error of the truncation of the perturbative series can underestimate the error. One such standard methods is to take the magnitude of the last term computed (or this quantity multiplied by α). Here, working in the lattice scheme, we are indeed in such a situation. The magnitude of the next order of the perturbative series is much smaller than the real error of the computation (the difference between the exact result and the truncated sum to order N ), as we move away from the first few orders. We illustrate this in Fig. 4. The reason is that each new term of the perturbative series is only marginally smaller than the previous one. This has an important additive effect before reaching the asymptotic regime. Other renormalization schemes of α (closer to the MS scheme) are expected to work better in this respect (with a smaller ratio between consecutive orders before reaching the asymptotic regime).
*) Precision (0, N P ) We now want to determine the error of the superasymptotic approximation, which we quantitatively discuss. We first give the number obtained from the fit, as well as the The error displayed for the red points is the complete error (statistical plus systematic combined in quadrature) of the p N coefficient obtained in [2] times α N +1 . The black diamond stands for the numerical minimal value of p N α N +1 . The black triangle is p N P α N P +1 using the smallest positive c that makes N P to be integer in Eq. (1.11). Note that the plus/minus error does not display symmetrically in the plot because of the logarithmic scale. errors: The first error is the statistical error of the fit. The following errors are systematic. The second error is the error associated to different infinite volume extrapolations of the coefficients p n . Up to this point, the error runs parallel to the error analysis made in [28]. Nevertheless, unlike in this reference, we do the fit in the range β ∈ [5.8, 6.65]. If we do the fit in the range β ∈ [6, 6.65], as it was done in that reference, the result is -0.04 smaller, a small shift. This is the third error in Eq. (2.4). For both ranges the reduced χ 2 are similar: 0.44 and 0.42 for the range β ∈ [5.8, 6.65] and the range β ∈ [6, 6.65], respectively. On the other hand truncating the partial sum at the numerical minimal term yields, G 2 = 3.18 r −4 0 with χ 2 red = 0.69 for the range β ∈ [6, 6.65], and G 2 = 3.05 r −4 0 with χ 2 red = 1.28 for the range β ∈ [5.8, 6.65]. Looking to the points in Fig. 2, we also observe that the remaining a dependence can not be assigned to a specific slope that can be interpreted as an O(a 2 Λ 2 ) effect, since the sign would flip if taking the sum truncated at the minimal term determined numerically or using Eq. (1.11) (though the latter seems to yield a flatter curve). Therefore, with the superasymptotic precision, we cannot isolate O(a 2 Λ 2 ) effects. The fourth error is the difference of the fit truncating C G to O(α 2 ) or to O(α 3 ). The change is significant. This seems to be due to the low convergence of the weak-coupling expansion in the lattice scheme. We have checked that there is convergence (albeit slow) by including higher-order terms of the beta-function using the estimates obtained in [30]. Following [28], we assign a 2.5% error for the conversion from a to r 0 units. This is the fifth error in Eq. (2.4). The last error is the estimate of the higher-order terms in the hyperasymptotic expansion not included in the superasymptotic approximation. This error has been discussed before. This last error is taken as the difference between the fits including or not the leading terminant. This basically gives the same error than considering the difference of doing superasymptotic fits truncating the perturbative sum at the numerical minimal term or using Eq. (1.11). Other possible ways to estimate the error (like taking c to be negative such that N P changes by one unit) give smaller errors. This error is by far the major source of uncertainty of the superasymptotic approximation. In the last equality in Eq. (2.4) we have combined all these errors in quadrature.

*) Precision (4, 0)
We now add the terminant to the superasymptotic approximation. We obtain The error analysis follows to a large extent the error analysis of the superasymptotic approximation. The first error is the statistical error of the fit, with a smaller than one reduced χ 2 : χ 2 red = 0.43. The following errors are systematic. The second error is the error associated to different infinite volume extrapolations of the coefficients p n . We emphasize again that we do the fits over the whole range β ∈ [5.8, 6.65]. If we do the fit in the range β ∈ [6, 6.65], the result is +0.09 larger with a rather small χ 2 red = 0.019. This is the third error in Eq. (2.4). Having a look to the points in Fig. 2, the remaining a dependence is very small but may point to a small negative slope. If anything, this effect is only visible for the largest distances. At short distances, the a dependence is completely hidden by the errors, which reflects in this very small χ 2 red , but even at the largest distances, the errors hide any meaningful signal of these effects. Note that this possible remaining a dependence can be associated to higher-order terms of the hyperasymptotic expansion of S PV , which would then scale as O((aΛ) 4(1+log(3/2)) ) rather than to genuine nonperturbative corrections that would scale as O(a 6 Λ 6 ). In this respect, it is worth noting that this small slope somewhat tends to disappear as we work with precision (4, N ), albeit with a huge error (see Fig. 3). The fourth error is the difference of the fit truncating to O(α 2 ) or to O(α 3 ) the perturbative expansion of C G . The fifth error is the one associated to the conversion from a to r 0 units. The last error is the error associated to Z P , the normalization of the leading renormalon. The error of this quantity is heavily correlated to the knowledge of the coefficient b 2 (see the discussion in [2]). Therefore, to estimate this error, we correlate the change of Z P to setting b 2 = 0. Comparatively to this error, the subleading terms of the weak coupling expansion in Eq. (1.17) produce a smaller change and can be neglected. We now discuss the error associated to the truncation of the hyperasymptotic approximation. As discussed before, the leading contributions to the hyperasymptotic expansion of S PV that are not included in the (4,0) precision are expected to scale as O((aΛ) 4(1+log(3/2)) ) and to be suppressed by a factor O((aΛ) 4 log(3/2) ) (times α) with respect to the typical size of Ω G 2 . This produces corrections at the level of the one/two MeV level. Therefore, unlike in the case of the superasymptotic approximation, the errors of the hyperasymptotic approximation (4,0) are small and can be considered to be included in the other errors, in particular in those associated to Z P and the different ranges we do the fits, as they measure our incomplete knowledge of the perturbative expansion (independently of the truncation of the perturbative expansion in C G ). Finally, we combine all the errors in quadrature producing the last equality in Eq. (2.5). This is our most precise prediction for G 2 PV , which we display in Fig. 2.
The central value we obtain does not change much with respect to the central value obtained [28]. Nevertheless, this is to some extent by accident, as the fit is made over different energy intervals. On the other hand the superasymptotic approximation truncated at the numerical minimal term appears to approach better the central value 4 . We saw that also for the self-energy of the static quark [30]. Nevertheless, the error is larger because the points are more scattered around, and because of the intrinsic inaccuracy of the superasymptotic approximation. In our case, the total error is basically shrunk by a factor 1/2. Note that the statistical error and the error associated to the infinite volume extrapolation of the coefficients are smaller now. The improvement in the quality of the fit can also be observed by the flatter curve we have now.
*) Precision (4, N ) We may try to increase the accuracy reached with the (4,0) hyperasymptotic approximation by adding the last term of Eq. (1.10). Nevertheless, the errors quickly grow and get out of hand. This is mainly due to the error of the coefficients of the perturbative expansion. We have repeated the same error analysis than in the previous item for N = 31, 32, 33, 34. We show the obtained central values and errors in Fig. 5. We see how the errors quickly grow. Actually, the most important source to the error comes from the infinite volume extrapolation of the perturbative coefficients p n .

Discussion and conclusions
In this paper, we have given the hyperasymptotic expansion of the plaquette with a precision that includes the terminant associated to the leading renormalon. Subleading effects are also considered. The perturbative series is regulated using the PV prescription for its Borel integral. We use this analysis to give a determination of the gluon condensate in We emphasize that this result is independent of the scale and renormalization scheme used for the coupling constant. Even if the computation was made in the lattice scheme, the result is the same in the MS scheme within the accuracy of the computation. G 2 PV (n f = 0) was computed with superasymptotic approximation in [28]. Here we have improved over this determination, principally by including the terminant associated to the leading renormalon. Adding Ω G 2 elliminates the jumps one has when using the superasymptotic approximation. The result is now much more smooth and flatter, to the point that, within errors, we can not isolate O(a 2 Λ 2 ) effects. We still observe some small bending, which could also be due to higher-order perturbation theory. Overall, we are able to shrink the error by around a factor 2.
In the lattice scheme, the impact of adding Ω G 2 is small compared with the size of the NP gluon condensate regulated using the PV prescription, of order 10%. In the MS scheme, the contribution of the terminant would be larger by a factor α MS /α latt , which could easily enlarge the contribution by a factor 2. Note though that these statements are dependent on the value of c used to fix N P .
As we have seen in this paper, at present, the limiting factor for improving the determination of the gluon condensate in pure gluodynamics is the error of perturbation theory. All systematic sources of error have its origin in the errors of perturbation theory (even what we call statistical errors of Eq. (2.1) are dominated by the statistical errors of the coefficients p n ). More precise values of these perturbative coefficients, and its knowledge to higher orders, would yield a more precise determination of the normalization of the renormalons, Z P , and would allow working with hyperasymptotic accuracy (4, N ). Nowadays, if we try to reach this accuracy, we find that the error of the coefficients are too large to get accurate results. The situation with active light quarks is in an early stage but starts to be promising. The coefficients of the perturbative coefficients have been computed at finite volume in [35] for QCD with two massless fermions. More data at different volumes, and the infinite volume extrapolation of these coefficients, would then allow to give a determination of the gluon condensate in QCD with two massless fermions.
It is interesting to show how our results fit general expectations for superasymptotic and hyperasymptotic expansions. Fig. 6 nicely display, for a four-dimensional gauge theory, the standard behavior expected for perturbative series that are asymptotic to an observable (we take β = 6 for illustrative purposes). We first discuss the blue points. First, as we add more terms to the perturbative series, such a perturbative series gets closer to the MC simulation of the plaquette. Nevertheless, as it is also expected for an asymptotic series, the rate of convergence diminishes till reaching an inflection point. This point is not a n ]α n+1 ) (black squares). For N = N P we draw P MC − ( N P n=0 p n α n+1 + Ω G 2 + π 2 36 C G (α) a 4 G 2 PV ) (red diamonds). For N > N P we draw P MC − ( N P n=0 p n α n+1 + Ω G 2 + π 2 36 C G (α) a 4 G 2 PV + N n=N P +1 [p n − p (as) n ]α n+1 ) (red diamonds). The error in all cases is the statistical error of the sum N n=0 p n α n+1 and P MC combined in quadrature. Note that the plus/minus error does not display symmetrically in the plot because of the logarithmic scale, and also because of the logarithmic scale the error looks different for different points located at the same N . In the small box a zoom of the points for N ≥ 27 are shown in non-logarithmic scale. minimum. Adding extra terms to the perturbative series one approaches to the observable even if the perturbative series is divergent here. The reason one does not reach the minimum is the non-zero value of the gluon condensate using the PV prescription. If we also subtract the gluon condensate plus the perturbative expansion, we are in the same situation than in the large-β 0 models that were studied in [29][30][31], where the nonperturbative contribution is zero by construction. See Figs. 7 in [29,30] for illustration. We then see the minimum (the maximal accuracy reached by the theory and how the series deteriorates if one continues adding extra perturbative terms). If at this points one adds Ω G 2 and the modified perturbative series where the leading renormalon is subtracted, one gets a plateau, and one can determine the gluon condensate. We illustrate this behavior for β = 6 in Fig. 7. If one also subtracts the gluon condensate, as we do in this figure, one also has the jump in the precision achieved, as also observed in the large-β 0 plots (see again Figs. 7 in [29,30] for illustration).
We finally mention that the nonzero value of G 2 PV shows that the PV regularization of the perturbative sum, even if computed exactly, would differ from the Montecarlo simulation of the plaquette by a term of O(a 4 Λ 4 ). This may affect the conjecture that the resummation technique of the perturbative expansion proposed in [38] for the Adler function would not need such nonperturbative corrections. This should be further investigated.