Gradient flow, renormalon ambiguity, and the gluon condensate

We study the renormalon ambiguity in the perturbative expansion of~$\langle E(t)\rangle$, the vacuum expectation value of the `energy density operator' $E(t,x)=\frac{1}{4}G_{\mu\nu}^a(t,x)G_{\mu\nu}^a(t,x)$ defined by the Yang--Mills gradient flow, within the so-called large-$\beta_0$ approximation. The leading uncertainty turns out to be of the order of~$\Lambda^4$, where $\Lambda$ is the dynamical scale. We argue that the perturbative expansion of~$\langle E(t)\rangle$ is identified with the coefficient of the identity operator, $c_{\mathbbm{1}}(t)$, in the small flow time expansion of~$\langle E(t)\rangle$. Hence, the higher order terms in the expansion, such as the gluon condensate, are ambiguous unless one properly separates the renormalon ambiguity from~$c_{\mathbbm{1}}(t)$. Based on this reasoning, we give a crude estimate of the gluon condensate defined in renormalon-free and observable-independent ways, after we analytically remove the renormalon ambiguity from~$c_{\mathbbm{1}}(t)$. For this estimate, we use the lattice data of~$\langle E(t)\rangle$ in the quenched QCD.


Introduction
Each coefficient in the perturbative expansion in quantum field theory typically grows factorially as a function of the order of perturbation theory and thus the perturbation series is an asymptotic series at best [1].One of the origins of this growth is a factorial increase of the number of Feynman diagrams with respect to the order.In renormalized perturbation theory, there is another origin: there exists a class of Feynman diagrams whose amplitude grows factorially [1,2].This kind of factorial behavior produces the so-called renormalon ambiguity in perturbation theory of order e −4πu/(β0α) ∼ (Λ 2 /µ 2 ) u , where β 0 is the one-loop coefficient of the beta function, α the renormalized coupling constant and u parametrizes the "strength" of the renormalon; Λ is the renormalization group invariant mass scale and µ a typical energy scale in the problem under consideration.We assume that β 0 > 0. Then this infrared (IR) renormalon exhibits a difficulty in extracting the low-energy physics at µ ∼ Λ from the perturbation series [3].
In this paper, we consider the ambiguity induced by IR renormalon in the perturbative expansion of E(t) , the vacuum expectation value of the "energy density operator" E(t, x) ≡ 1 4 G a µν (t, x)G a µν (t, x) defined by the Yang-Mills gradient flow [4,5].The Yang-Mills gradient flow is a one-parameter evolution of the gauge field A µ (x) according to the flow equation, 1 where t ≥ 0 is called the flow time and the initial value of the evolution is given by the original gauge field A µ (x); G µν (t, x) is the field strength of the flowed gauge field B µ (t, x), and the covariant derivative is also defined with respect to B µ (t, x), Using the field strength (1.2), the gauge invariant "energy-density operator" is constructed as G a µν (t, x)G a µν (t, x). (1.4) As the renormalizability theorem [6] implies, its vacuum expectation value is a renormalized finite quantity although it is a certain combination of the bare gauge fields through the flow equation. 2This physical observable is thus quite useful for the scale setting and the non-perturbative definition of the gauge coupling in the context of lattice gauge theory; see a review [8] and a recent paper [9] and references cited therein.
In this paper, we study the IR renormalon ambiguity in the perturbative expansion of E(t) (1.5) using the so-called large-β 0 approximation [10,11].The size of the ambiguity turns out to be O(Λ 4 ).We then extract a renormalon-free part from the perturbative expansion of E(t) by employing the technique of Refs.[12,13] (see also Ref. [14]).Subsequently, we argue that the perturbative expansion of E(t) is identified with the coefficient of the identity operator in the small flow time expansion [6] of E(t) (1.5), c ½ (t).This implies that other terms in the small flow time expansion of E(t) , such as the gluon condensate {F a µν F a µν } R (x) , are ambiguous when it is naively calculated as E(t) − c ½ (t).This suggests that one should separate the renormalon uncertainty from c ½ (t) in order to extract unambiguous condensates.On the basis of this reasoning, we give a crude estimate of the gluon condensate in the quenched QCD, using the above renormalon-free part of c ½ .
Note that the renormalon ambiguity is always problematic in determination of the gluon condensate, not only in our current case considering the small flow time expansion of E(t) but also in other cases, for instance, considering the operator product expansion (OPE) of vector currents.A significant effect of the renormalon ambiguity on the gluon condensate has already been emphasized in Ref. [15].See also Ref. [16].In view of this situation, determinations of the gluon condensate while avoiding the renormalon ambiguity have been performed in recent studies from the plaquette [15,[17][18][19].
Characteristics of the present work can be stated as follows.We first give an explicit definition of the gluon condensate where renormalon-free nature and observable independence are explicit.Namely, the gluon condensate is defined as a universal quantity, which is consistent with its original concept.Then, we will show that this renormalon-free condensate indeed appears in the small flow time expansion after absorbing the renormalon uncertainty of c ½ .(The renormalon uncertainty of c ½ is separated out in the above extraction of the renormalon-free part.)Therefore, a proper formalism to render the condensate free from renormalon is explicitly demonstrated, although our framework can be applied within the large-β 0 approximation at this stage.
Based on the above renormalon-free small flow time expansion, we give a crude estimate of the gluon condensate in the quenched QCD from the lattice data of E(t) .Unfortunately, it turns out that it is difficult to obtain a reliable estimate of the gluon condensate within the present work because it is not clear how to systematically improve the large-β 0 approximation.Nevertheless, we believe that it is worthwhile to report our analysis because we demonstrate how we can define an unambiguous condensate by properly subtracting renormalon ambiguity.

Renormalon ambiguity in the perturbative expansion of E(t) 2.1. Computation of E(t) in the large-β 0 approximation
We study the ambiguity induced by the IR renormalon in the perturbative expansion of Eq. (1.5).For this, we have to consider perturbation theory to infinite order in some way.A convenient framework is the so-called large-β 0 approximation [10,11] that proceeds in our present problem as follows.
First, to extract a gauge invariant subset of Feynman diagrams that gives renormalon, we consider the large-n f approximation (n f ≫ 1), where n f is the number of fermion flavors in a gauge theory, while g 2 0 n f is held fixed; g 0 is the bare gauge coupling.With our notational convention in Appendix A, the bare propagator of the gauge field is then O(1/n f ) and vertices among the gauge fields (including the ones induced by a fermion loop) are O(n f ).In this approximation, the leading O(1/n f ) gauge field propagator is then given by 3 where λ 0 is the bare gauge fixing parameter.In this expression, the factor [1 − ω(p)] −1 arises from the geometric sum of fermion loop chains in Fig. 1, where ω(p) is the vacuum polarization given by and γ E is the Euler constant.From Eqs. (2.2) and (2.3), we see that the renormalization in the MS scheme is accomplished by Then the renormalized O(1/n f ) gauge field propagator is given by (2.5) Note that there is no need of the wave function renormalization of the gauge field at the order we consider.A simple order counting also shows that an n-point correlation function of the gauge fields is O(1/n n−1 f ).
3 We adopt dimensional regularization in which the spacetime dimension is set to be D ≡ 4 − 2ǫ.We also use the abbreviation, The large-n f approximation can also be considered for correlation functions of the flowed gauge fields defined by Eq. (1.1).The formal solution of Eq. (1.1) is given by [5] where is the heat kernel and ).Thus, the leading O(1/n f ) flowed gauge field propagator B a µ (t, x)B b ν (s, y) is given by, after the substitutions z) and A b σ (w) by Eq. (2.2); the contribution of the non-linear term R µ (2.8) always increases the number of A and thus lowers the power of n f .In this way, the flowed gauge field propagator in the large-n f approximation is given by (2.9) The parameter α 0 does not receive the renormalization [6,7].The large-n f expression of E(t) (1.5) is then simply given by contracting two gauge fields in E(t, x) by the propagator (2.9); it is easy to see that the other Feynman diagrams that potentially contribute to E(t) always lower the powers of n f .The contraction yields This is the expression in the leading order of the large-n f approximation.Now, the large-β 0 approximation [10,11] is simply defined by replacing the factor − 43 T (R)n f in the above expression by the one-loop coefficient of the beta function, That is, in this large-β 0 approximation, E(t) is given by where in the last expression we have set τ ≡ p 2 and Equation (2.12) is our basic formula.It has been known that this large-β 0 approximation reproduces the leading logarithmic terms in the actual perturbation series; see Appendix B for this point.

Renormalon singularity in the large-β 0 approximation
By expanding Eq. (2.12) with respect to α, we have the perturbative series in the large-β 0 approximation, where k0 = 1 and k1 and k2 are explicitly given in Eqs.(B6) and (B7) in Appendix B, respectively.In Appendix B, these perturbative coefficients in the large-β 0 approximation are compared with the exact perturbative coefficients obtained in Refs.[5,20].We may consider the Borel transform corresponding to the perturbative series (2.14), The singularities of the Borel transform are located at u ≡ β 0 b/(4π) = 2, 3, 4, . . ., while the so-called ultraviolet (UV) renormalons (singularities at negative b) do not exist because of the UV behavior improved by the gradient flow.Due to the pole singularities at b > 0 in Eq. (2.15), the Borel sum for is ill-defined.This quantity is often regularized by deforming the integration contour in the complex-b plane as ∞ 0 db → ∞±iε 0±iε db so that it does not hit the singularities.Considering the singularity closest to the origin at u = 2, i.e., the leading IR renormalon, the difference between the two regularized integrals is (2.17)This is the u = 2 renormalon uncertainty.Here Λ is the one-loop Λ parameter in the MS scheme, Λ 2 ≡ µ 2 e −4π/(β0α) . (2.18)

Extraction of a renormalon-free part in the perturbative series of t 2 E(t)
Going back to Eq. (2.12), we write it as by introducing the function w(x), and the effective gauge coupling in the large-β 0 approximation, (2.21) The integral Eq. (2.19) as it stands is ill-defined because the integrand possesses a pole at the IR scale τ = e5/3 Λ 2 , at which the effective coupling α β0 (τ ) (2.21) diverges.We see that the difference between the integrals regularized in a parallel way ( ∞ 0 dτ → ∞±iε 0±iε dτ ) gives Eq. (2.17) as the leading contribution.Therefore, this ambiguity associated with the pole can be regarded as the effect of IR renormalon.
To render the integral well-defined, we introduce a "factorization scale" µ f such that 5 Λ 2 ≪ µ 2 f ≪ 1/t and define This integral over UV modes can be reliably evaluated by perturbation theory.Its dependence on µ f , however, indicates the limitation of perturbation theory because the original physical quantity E(t) is independent of the artificial cutoff µ f .In this respect, the µ f dependent part corresponds to a renormalon uncertainty.
To extract a µ f -independent part from Eq. (2.22), which is regarded as a renormalonfree part, we follow the procedure in Refs.[12,13].We introduce the function defined in a complex plane, Through an auxiliary function, where Γ (a, z) ≡ ∞ z dt t a−1 e −t is the incomplete Gamma function, we have Then using Eq.(2.24) in Eq. (2.22), we have where the integration contours C a and C b are displayed in Fig. 2. In the last line of Eq. (2.27), the first integral along C a is finite and independent of the cutoff µ f (see Fig. 2).We can further deform the contour where the function W + (x) is given in Eq. (2.25).
In the integral along C b in the last line of Eq. (2.27), we may use the expansion in Eq. (2.26) because tµ 2 f ≪ 1.The contribution of the first few terms in the expansion (2.26) turns out to be independent of µ f because their coefficients are real and the contribution can be evaluated by the integral along the contour surrounding the pole τ = e 5/3 Λ 2 [12,13], i.e.,

− C Im
4(e 5/3 tΛ 2 ) 2 ln(e 5/3 tΛ 2 ) + γ E + ln 2 . (2.29) The subsequent terms in the expansion (2.26) possess imaginary coefficients and their contribution cannot be independent of µ f ; the leading µ f -dependence for small t is estimated as Note that this is O(t 2 µ 4 f ).In this way, we separate the µ f dependence in Eq. (2.22) as where t 2 E(t) cut (µ f ) is given by Eq. (2.30) and we have obtained the renormalon-free part, that is independent of the cutoff µ f , by where (2.33) In Fig. 3, we plot the renormalon-free part t 2 E(t) RF (2.32) with Eq. (2.33) for G = SU (3) and n f = 0.

Small flow time expansion and the renormalon ambiguity
Equation (2.31) with Eq. (2.30) shows that the ambiguity in the perturbative expansion of t 2 E(t) is O(t 2 ) for small t.This is consistent with Eq. (2.17) that shows that the leading renormalon ambiguity in t 2 E(t) is O(t 2 ).The above analysis thus indicates that the O(t 2 ) term of t 2 E(t) suffers from the renormalon ambiguity and it cannot be determined by perturbation theory alone.
At this stage, we recall the small flow time expansion [6], an asymptotic expansion for small t of a composite operator of flowed fields in terms of un-flowed fields.For the energy density operator (1.4), it reads where ½ stands for the identity operator and {O} R (x) is a renormalized composite operator corresponding to a bare operator O(x).Adopting the MS scheme, the expansion coefficients Fig. 3: The renormalon-free part t 2 E(t) RF (2.32) with Eq. (2.33) for G = SU (3) and n f = 0 (green solid line).We also plot each contribution in Eq. (2.32): t 2 E(t) 0 (black dotted line), tΛ 2 -term (black dashed line), (tΛ 2 ) 2 -term (black dot-dashed line).
to the one-loop order are given by [21,22], ( The vacuum expectation value of Eq. (3.1) thus yields Now, we argue that the coefficient of the identity operator c ½ (t) is given by Eq. (2.31) within the large-β 0 approximation and thus suffers from the renormalon ambiguity.That is, we infer that First, this picture emerges from the similarity between the small flow time expansion and the OPE: the flow time t is analogous to (the square of) a separation r between two operators for which the OPE is considered.The expansion coefficients, c ½ (t), c F F (t), etc. in Eq. (3.1), are analogous to the Wilson coefficients in OPE.Here, let us recall that the integrationby-regions method [23,24], which is a technique to evaluate loop integrals by dividing loop momenta into hard and soft regions, gives the following correspondence in the context of OPE: The hard part carries the information of the external momentum (i.e., the separation r) and corresponds to the Wilson coefficient which can be evaluated by perturbation theory.The soft part, on the other hand, is independent of r and can be identified with the vacuum expectation value of an operator that is generally non-perturbative.Thus, the perturbative expression (2.22), which consists only of the "hard part," is naturally identified with the expansion coefficient c ½ (t). 6e can somewhat strengthen the above argument as follows: Going back to Eq. (2. 19), we write it as In the second integral over the soft momentum, if we include non-perturbative effects, the integrand V(tτ, tΛ 2 ) will be modified from that of the large-β 0 approximation (or some perturbative expression) so that the integral becomes well-defined.Although the explicit form of V(tτ, tΛ 2 ) is unknown, we may expand V(tτ, tΛ 2 ) in tτ (as well as tΛ 2 ), because τ is limited by µ f in the second integral.This expansion clearly produces the small flow time expansion of the form (3.5) in which each term of the expansion contains the integration over the soft region, which may be identified with the vacuum expectation value of an operator at t = 0 such as {F a µν F a µν } R (x).Although it is difficult to make this picture more quantitative, if we use the lowest order of perturbation theory, then V(tτ, tΛ 2 ) = 8πt 2 τ 2 e −2tτ α = w(tτ )α from Eqs. (2.12) and (2.20) and the second integral in Eq. (3.7) becomes This is precisely in the lowest perturbation theory, 7 provided that µ f is identified with the UV cutoff. 8This argument again supports the identification (3.6), that is, the first term of the right-hand side of Eq. (3.7) is identified with c ½ (t) in Eq. (3.5).
Another somewhat different support for Eq.(3.6) is given by applying the large-β 0 approximation to the both sides of Eq. (3.1).We denote the evaluation of the vacuum expectation value of an operator O in the large-β 0 approximation with the IR cutoff µ f as O (µ f ) like in Eq. (2.22).We assume that we set the IR cutoff for all integrals uniformly as Then, noting ½ (µ f ) = 1,9 we have Subtracting this from Eq. (3.5), ). (3.11)This shows that the soft part in t 2 E(t) , the second term of Eq. (3.7), does not include c ½ and is identified with O(t 2 ) terms (and higher order terms) in the small flow time expansion while µ f is regarded as the UV cutoff for vacuum expectation values.This is consistent with the statement (3.6), if we assume that the vacuum expectation values possess the UV cutoff µ f .From the arguments above, the identification (3.6) is quite plausible.We thus assume this in the following argument.
Equation (3.6) shows that the coefficient c ½ (t) suffers from the ambiguity in perturbation theory, i.e., the renormalon ambiguity, represented by the last term t 2 E(t) cut (µ f ).The lesson from this observation is that if one tries to extract the sum of the condensates, (3.12) from Eq. (3.5) by simply subtracting c ½ (t) from t 2 E(t) , which is non-perturbatively computed by the lattice simulation for example, the sum of the condensates also suffers from the renormalon ambiguity.We thus have to remove the renormalon ambiguities from both of c ½ (t) and the sum of the condensates.
We now explicitly examine a renormalon uncertainty of the condensates [the right-hand side of Eq. (3.11)].In our method using the cutoff, the renormalon uncertainty is exhibited as the UV cutoff dependence.Note that the UV cutoff should be introduced to the gluon loop momentum, corresponding to the definition of t 2 E(t) (µ f ).Then, the second and third condensates of Eq. (3.11), where corrections due to gluon loops are negligible from the viewpoint of the large-n f and thus large-β 0 approximation, do not have UV cutoff dependence at the order we consider.Hence, µ f dependence arises only from This cutoff dependence is regarded as the renormalon uncertainty of the gluon condensate.
Then, we propose a definition of the cutoff independent (renormalon-free) gluon condensate as Cb dτ iτ α β0 (τ ), (3.14)where the µ f dependence in Eq. (3.13) is exactly cancelled.Note that in this definition we do not use any information on observable.
Adopting the above definition, the renormalon-free condensates can be extracted from Here, we use the separation Eq. (3.6) in the small flow time expansion Eq. (3.11) and then we include t 2 E(t) cut (µ f ) of c ½ in the condensates, which reduces the condensate to the renormalon-free one.Note that c F F (t) ≃ 1.Note also that we can neglect the renormalons of the second and third condensates.This describes the procedure that one can find a renormalon-free condensate by appropriately absorbing the renormalon uncertainty of c ½ into the condensate.In the above equation, t 2 E(t) RF is given by Eqs.(2.32) and (2.33).We assume the limit t → 0 to get rid of the O(t 4 ) term in Eq. (3.5).In Eq. (3.15), the full t 2 E(t) should be computed non-perturbatively, for instance, by employing the lattice simulation.This formula tells how to extract the 4-dimensional condensates from the 5dimensional gluon condensate while taking into account the renormalon issue (under the assumption of the large-β 0 approximation).One might think that our definitions of the renormalon-free part given in Eqs.(2.31) and (3.14) are ad hoc: one may add an arbitrary function f (t) = O(t 2 ) that is independent of µ f to t 2 E(t) RF while simultaneously subtracting it from t 2 E(t) cut (µ f ) and from the definition of [ 1 4 {F a µν F a µν } R (x) ] RF .This redefinition changes the value of the condensate.This is certainly true.Namely, there are potentially many schemes to define renormalonfree parts.This originates from the fact that one can freely choose a way to specify what the renormalon uncertainty is.This is inevitable in defining a quantity which suffers from a renormalon uncertainty. 10In this sense, our definition (3.15) adopts a particular scheme.In our scheme, the gluon condensate is defined while the cutoff dependence is minimally subtracted based on the contour deformation.We note that the renormalon-free quantities defined in different schemes should not be compared directly.
We emphasize that our definition of the renormalon-free condensate has a universal meaning, once the scheme is fixed, in the following sense (within the large-β 0 approximation).Suppose that one considers the OPE for a dimensionless observable X with a typical scale Q, Then, within the large-β 0 approximation, the perturbative part c X ½ (Q) defined with the IR cutoff µ f can generally be separated into cutoff dependent part given by and the rest cutoff-independent part.Note that in Eq. (3.17), X-dependence appears in the overall factor of the integral whereas the integral itself is independent of X.This is because the µ f -dependence precisely cancels the UV cutoff µ f dependence of the gluon condensate irrespective of the chosen observable X. Equation (3.17) is certainly true for X = t 2 E(t) as we have seen above (c X ½ (Q) cut (µ f ) = t 2 E(t) cut (µ f ) in this case).For the static QCD potential X = rV QCD (r), Eq. (3.17) was shown in Ref. [25].Our renormalon-free definition of the condensate [Eq.(3.14)] is thus expected to give an identical condensate for any observable X within the large-β 0 approximation.

Gluon condensate in the quenched QCD
We try to estimate the gluon condensate by applying the formula (3.15) to the quenched QCD (i.e., n f = 0).In the quenched case, the equation becomes simple: where t 2 E(t) RF is given by Eqs.(2.32) and (2.33).This relation should be regarded as an asymptotic one holding in the limit t → 0. As we will see below, however, it turns out to be difficult to extract a reliable value of the gluon condensate directly employing the above formula.This is because our theoretical calculation heavily relies on the large-β 0 approximation, which includes only a partial set of Feynman diagrams.Nevertheless, we believe that it is useful to report our examination on the formula (4.1), as an attempt towards a possible systematic separation of the renormalon ambiguity (e.g., on the basis of the explicit perturbative expansion) in the gluon condensate.
We use lattice data obtained by the FlowQCD collaboration [26,27]. 11In Fig. 4, we show the lattice results for t 2 E(t) for the bare gauge couplings, β = 6.4,6.6, 6.8, 7.0, and 7.2.To show the lattice data in Λ 3-loop MS units, we used the relation between β and the lattice spacing a obtained in Ref. [27]. 12We see that the lattice data among different β's overlap each other in the region t(Λ 3-loop MS ) 2 0.01.Therefore, we use the lattice data with the finest lattice spacing (β = 7.2 and a 2 (Λ 3-loop MS ) 2 = 5.3 × 10 −4 ) shown by the black line in Fig. 4 for t(Λ 3-loop MS ) 2 ≥ 0.01 regarding it as the continuum limit. 13he lattice results in Fig. 4 and the renormalon-free part in Fig. 3 qualitatively look similar.To compare them quantitatively, we need the ratio r ≡ Λ 1-loop MS /Λ 3-loop MS , because our theoretical calculation on the basis of the large-β 0 approximation is given in Λ 1-loop MS units whereas the lattice results are shown in Λ 3-loop MS units.We determine this ratio by requiring the running couplings at 1-loop and 3-loop to have the same value at µ = a −1 , i.e., we impose α s,1-loop (a −1 /Λ 1-loop MS ) = α s,3-loop (a −1 /Λ 3-loop MS ) = 0.1214.(Note that the running coupling at kloop is a function of µ/Λ k-loop MS .)This condition assures that the calculation at leading-log (LL) matches well with the one at next-to-next-to-LL (NNLL) around the region t(Λ 3-loop MS ) 2 ∼ a −2 (Λ 3-loop MS ) 2 = 5.3 × 10 −4 .14This is legitimate because both predictions should be accurate in such a short distance region.The above condition yields r = 0.395.
Using the above r, in Fig. 5, we show the comparison between the lattice result and the renormalon-free part obtained in the large-β 0 approximation, t 2 E(t) RF in Eq. (2.32).The difference between them, shown in the upper right panel, is expected to have a linear behavior in t 2 (Λ 3-loop MS ) 4 according to Eq. (4.1).To investigate quantitatively if this is the case or not, in the lower panel, we plot an effective power of the difference defined by d ln f (x)/d(ln x), where f (x) is the difference and x ≡ t(Λ 3-loop MS ) 2 .From the lower panel, it seems that a component with the power smaller than 2 remains in the difference, that is, the right-hand side of Eq. (4.1) does not behave as O(t 2 (Λ MS ) 4 ).This failure is attributed to the large-β 0 approximation.
As a possible next best, we try to utilize the exact perturbative series for t 2 E(t) up to O(α 3 s ) [5,20] to improve the matching of t 2 E(t) RF with the actual lattice result.That is, we modify t 2 E(t) RF such that its asymptotic expansion coincides up to O(α 3 s ) as with δk i ≡ (k i − ki ) µ=1/ √ t for i = 1 and 2; k i are exact perturbative coefficients and ki are those in the large-β 0 approximation [Eq.(2.14)].See Appendix B for k i .In Fig. 6, we compare Fig. 4: Lattice results for t 2 E(t) .Different colored lines correspond to different β.The line width represents the statistical error.Fig. 5: Comparison of t 2 E(t) RF (2.32) with the lattice result.In the upper left panel, t 2 E(t) RF (green line) and the lattice data (black line) are shown together as functions of t(Λ 3-loop MS ) 2 .The difference between them is shown in the upper right panel, where the horizontal axis is taken as t 2 (Λ 3-loop MS ) 4 in order to examine if the linear behavior indicated by Eq. (4.1) is observed.In the lower panel, we show an effective power of the difference in t(Λ 3-loop MS ) 2 .Statistical error is not estimated in this panel.In the upper right and lower panels, we show only the data points in the region t(Λ 3-loop MS ) 2 0.01.
Eq. (4.2) with the lattice result.In the right panel, the difference between them looks linear in t 2 (Λ 3-loop MS ) 4 , which is now consistent with the expectation from the small-t expansion [Eq.(4.1)].We fit the data at 0.01 < t(Λ 3-loop MS ) 2 < 0.03 by the function At 2 (Λ 3-loop MS ) 4 and obtain the fitting parameter as A = 36.This yields, setting c F F (t) = 1 for simplicity, or, following the normalization in Ref. [15] (here −β(α)/(πβ These numbers should not be taken so seriously because these are obtained by an ad hoc mixture of the large-β approximation and the exact perturbation series.Moreover, if we adopt µ = 1/ √ 8t (instead of µ = 1/ √ t) as the renormalization scale in δk i as often done in the literature, the difference corresponding to the right panel of Fig. 6 takes a very different form and the coefficient of the quadratic curve becomes quite unstable; A 2 takes values A 2 = 10-39 if we change the order of polynomial by adding a t-independent constant or t 3 (Λ 3-loop MS ) 6 term as a test.
Fig. 6: Comparison of Eq. (4.2) with the lattice result.In the left panel, t 2 E(t) corr RF (blue line) and the lattice data (black line) are shown together as functions of t(Λ 3-loop MS ) 2 .For comparison, t 2 E(t) RF in the large-β 0 approximation is also shown by the green line.The difference between the lattice result and Eq.(4.2) is shown in the right panel, where the horizontal axis is taken as t 2 (Λ 3-loop MS ) 4 .The linear function determined by a fit is shown by the red line.The upper limit of the range used in this fit is shown by the dotted line.

Conclusion
In this paper, we studied the renormalon ambiguity in the perturbative expansion of t 2 E(t) , the vacuum expectation value of the "energy density operator" E(t, x) defined by the gradient flow.On the basis of the large-β 0 approximation, we obtained a renormalon-free part in the perturbative expansion, t 2 E(t) RF , by following the procedure in Refs.[12,13].We then identified the perturbative expansion of t 2 E(t) with c ½ (t), the coefficient of the identity operator in the small flow time expansion of E(t, x).Since c ½ (t) is the difference between the exact t 2 E(t) and the sum of vacuum condensates of dimension 4 operators (for small t), t 2 E(t) − t 2 E(t) RF provides a possible method to compute the condensate without suffering from the renormalon ambiguity.We emphasize that this formulation is not just a prescription for removing renormalon ambiguity, but based on a definition of the renormalon-free gluon condensate that is independent of physical processes (within the large-β 0 approximation).We applied this method for estimating the gluon condensate in the quenched QCD by using lattice data of E(t) .Unfortunately, we could not obtain a reliable number for the gluon condensate in this approach; this failure would stem from the fact that the large-β 0 approximation cannot sufficiently well-simulate the real behavior of the perturbative expansion.
In spite of this failure, we believe that it is worthwhile to stress the importance of the separation of the renormalon ambiguity from the perturbative expansion.The same problem exists in other theoretical approaches to the gluon condensate which are based on a combination of perturbation theory and non-perturbative calculation such as the lattice.To pursue our approach further, we need some systematic method to extract the renormalon ambiguity order by order by using the actual perturbative series for t 2 E(t) .Such a method is available for the static QCD potential [28] for which a representation similar to Eq. (2.19) holds exactly.We hope to come back this issue in future.and [20]  Here β 0 is given by Eq. (2.11) and β 1 is the two-loop coefficient of the beta function, The perturbative coefficients in the large-β 0 approximation (2.Comparing k 1 (B3) with k1 (B6) and k 2 (B4) with k2 (B7), we see that the leading logarithmic terms, i.e., the O(L) term in k 1 and the O(L 2 ) term in k 2 , are correctly reproduced in the large-β 0 approximation.Also, we see that the leading large-n f terms, the O(n f ) term in k 1 and the O(n 2 f ) term in k 2 , are correctly reproduced; this is expected because the large-β 0 approximation becomes exact in the large-n f limit.

Fig. 2 :
Fig. 2: The integration contours C a and C b used in Eq. (2.27).