The perturbative QCD gradient flow to three loops

The gradient flow in QCD is treated perturbatively through next-to-next-to-leading order in the strong coupling constant. The evaluation of the relevant momentum and flow-time integrals is described, including various means of validation. For the vacuum expectation value of the action density, which turns out to be a useful quantity in lattice calculations, we find a very well-behaved perturbative series through NNLO. Quark mass effects are taken into account through NLO. The theoretical uncertainty due to renormalization scale variation is significantly reduced with respect to LO and NLO, as long as the flow time is smaller than about 0.1 fm.


Introduction
QCD is a remarkable theory in many respects. It has been unchallenged in the description of strong interactions since its original formulation of more than 40 years ago [1]. It crucially impacts theoretical descriptions of a vast range of observations, reaching from the hadron mass spectrum to cross sections at particle colliders. When quarks are neglected (quenched QCD), the only fundamental parameter of QCD is the strong coupling constant α s , which simultaneously defines a mass scale Λ QCD ∼ O(100 MeV) due to its momentum dependence implied by quantum field theory (dimensional transmutation). Also, due to asymptotic freedom [2,3], QCD does not have an ultra-violet cut-off; it remains consistent up to arbitrarily high energies.
A very successful approach for calculations based on the QCD Lagrangian is perturbation theory. It corresponds to an expansion in the strong coupling constant α s and is applicable for processes where the typical mass scale Q is much larger than Λ QCD . The natural input quantity is therefore α s , whose numerical value at a reference scale Q 0 Λ QCD is determined by comparing theoretical predictions with measurements for known processes. Any dependence of physical observables on Λ QCD is only implicit through α s : where β 0 is the leading order coefficient of the QCD β function and will be defined below.
Perturbation theory is obviously inadequate for calculating observables like the hadron mass spectrum or the pion decay constant which are strongly dependent on Λ QCD . Such problems are accessible in lattice gauge theory, however [4]. In this approach, space-time is discretized by a characteristic lattice spacing a which serves as a UV regulator.
Unfortunately, lattice gauge theory and perturbation theory are not only complementary approaches to QCD, but there is practically no overlap region where both would yield competitive results. There is a certain amount of cross-fertilization though, in particular in the context of renormalization [5][6][7][8] or flavor physics (for a review, see Ref. [9], for example).
A particularly promising theoretical quantity which is accessible both perturbatively and on the lattice is the so-called Yang-Mills gradient flow [10][11][12][13]. For a particular gauge invariant quantity (the so-called QCD action density, to be defined in more detail below), it was shown by Lüscher [12] that it exhibits some welcome features on the lattice which allows its efficient evaluation with rather high precision. He also explicitly calculated this quantity perturbatively through next-to-leading order (NLO) and showed that the standard QCD renormalization of the gauge coupling constant is sufficient in order to obtain a finite result [12]. This property was later proven to all orders in perturbation theory [13]. The perturbative and the lattice result were found to be compatible over a significant interval of the so-called flow-time parameter t, thus opening up a wide range of possibilities for cross-fertilization in both fields.
In lattice QCD, the benefits and the appeal of the gradient flow have already been established, for example in its use for determining the absolute mass scale of a lattice calculation [12,14] ("scale setting", see Ref. [15], for example). From the perturbative point of view, on the other hand, the concept has received almost no attention since the original works of Refs. [12,13]. However, since one may expect rather precise results on the lattice in this framework, this may pose a challenge for perturbative calculations as well, possibly leading to interesting first-principle results for QCD.
With this motivation in mind, we are going to study the QCD action density in the framework of the gradient flow up to next-to-NLO (NNLO) in a perturbative approach. The perturbative expansion will be obtained via Wick contractions of the original field operators. This results in three D-dimensional momentum and up to four flow-time integrations over products of massless Feynman propagators times exponential factors involving loop momenta and flow-time integration variables. They are solved by sector decomposition and suitable numerical integration routines. While quark-mass effects will be neglected in the NNLO calculation, we show that they can be included through a simple one-dimensional integral at NLO.
The remainder of this paper is organized as follows. In the next section, after briefly recalling the flow-field formalism, the generation of the perturbative series and the evaluation of the resulting integrals is described. Additionally, we give a list of checks performed on our calculation, validating the numerical as well as the conceptional steps. In Section 3, we present the numerical value for the NNLO coefficient of the action density, which is the main result of this paper, and provide a brief analysis of the numerical effects. We conclude and give a short outlook on possible extensions and applications of this work in Section 4.

The flow field
The theoretical framework of our calculation is defined by the equations [12,13] where B a µ (t, x) is the flow field with space-time index µ and color index a, g 0 is the bare QCD coupling constant, f abc are the SU(3) structure constants, and A a µ (x) is the fundamental gauge field of QCD. The derivative ∂ µ is understood w.r.t. the D-dimensional Euclidean 1 space-time variable x, while t denotes the so-called flow-time. It is easily seen [12] that solutions of Eq. (2) for different values of λ are related by a t-dependent gauge transformation of the flow field B µ .
It was shown in Ref. [13] that the flow equation (2) can be written as the following integral equation in momentum space: whereK and × X (n,0) (q 1 , . . . , q n ) ab 1 ···bn The vertices X (n,0) read The fact that the first term in Eq. (3) is proportional to g 0 allows an iterative solution of that equation, which leads to an asymptotic series for B a µ , B a µ = n≥1 g n 0 B a n,µ .
With each power of g 0 , the number of fundamental gauge fields A a µ increases by one. Furthermore, B a n,µ involves terms with n/2 , n/2 + 1, . . . , n − 1 flow-time integrations, where n/2 denotes the greatest integer less than or equal to n/2.
Note that Eqs. (4) and (6) become particularly simple for λ = 0; for example, the lowestorder solution of the flow-field equation is simplyB a µ (t, p) = e −tp 2Ã a µ (p) in this case.

Calculation of the action density
The quantity to be computed in this paper is the vacuum expectation value of the action density, Since E(t, x) is gauge invariant, we are allowed to set λ = 0 in our calculation, which minimizes the number of integrals to be evaluated. The case λ = 0 will be considered as an important check of our calculation in Sect. 2.4.
The perturbative expansion of the vacuum expectation value is obtained by inserting the asymptotic expansion of the flow field B a µ as obtained in the previous section, and including higher orders of the fundamental perturbative vacuum, i.e.
where S QCD is the interaction part of the fundamental QCD action which depends on the fundamental gauge fields A a µ . From Eqs. (9) and (3) it follows that E = O(g 2 0 ); only the first term on the r.h.s. of Eq. (9) contributes at lowest order, while the second and the third term are of order g 4 0 (odd powers in g 0 vanish due to an odd number of fields in the matrix elements). However, all three terms on the r.h.s. of Eq. (9) contribute to higher orders as well, either through higher orders in the expansion of the B-fields, see Eq. (3), or through the perturbative expansion of the exponential in Eq. (10). The former case generally leads to an increase in the number of flow-time integrations, while the latter corresponds to corrections due to fundamental QCD. The general form of a matrix element to be evaluated at order g n 0 can therefore be symbolized by where B m i is the m th i coefficient of the asymptotic series in Eq. (7). This classification turns out useful with respect to the way we subsequently simplify the matrix elements in the sense that individual terms cannot be combined among different classes. Note that 0 ≤ m − k is the maximum number of flow-time integrations in M n (k, m), and since m ≤ n, the maximum number of flow-time integrations at order g n 0 is 2 n − 2. One particularly simple class when calculating E is M n (2, 2), which is fully determined by the (n − 2)-loop self-energy of the fundamental gluon field; this will serve as a welcome check of our calculation, see Sect. 2.4. At LO, M 2 (2, 2) is in fact the only class that contributes. At order g n 0 for n ≥ 4, one needs to evaluate 3(n − 2) classes, namely M n (k, m) with k ∈ {2, 3, 4} and k ≤ m ≤ n. Thus, at NNLO, there are twelve classes that contribute to E , and the maximum number of flow-time integrations is four. For comparison, at NLO there are six classes and at most two flow-time integrations.

Evaluation of the perturbative series
Except for the final numerical integration, all stages of the calculation were performed with the help of Mathematica [16]. Rather than following the diagrammatic method developed in Ref. [13], we directly implemented the Wick contractions of the gauge and quark fields after the iterative expansion of the flow fields according to Eqs. (3) and (5), the perturbative expansion of exp(−S QCD ) in Eq. (10), and the insertion of QCD-Feynman rules. The Dirac algebra is performed with the functionalities of FeynCalc [17] and color factors are calculated using ColorMath [18]; vacuum diagrams are discarded as required by the normalization factor in Eq. (10).
After these algebraic and symbolic manipulations, one ends up with integrals of the general form (at O(g 6 0 )) where are sets of integers, N ≤ 4, t 0 ≡ t, and the upper limits for the flow-time integrations are linear combinations of the other flow-time variables, t up . The momenta p 4 , p 5 , p 6 are linear combinations of the integration momenta p 1 , p 2 , p 3 . Quark-mass effects have been neglected in Eq. (12); at NLO, we will take them into account in Sect. 2.5.
Needless to say that in order to minimize computer time, it is important to identify integrals which differ by linear transformations of the loop momenta and flow-time integration variables at this stage, to cancel numerators with denominators in the integrals as far as possible, and to discard scale-less integrals which vanish in dimensional regularization. After these simplifications, the number of integrals of the form given in Eq. (12) is listed in Table 1, both split according to the classification defined in Eq. (11), and according to the number of flow-time integrations. For comparison, we also give the corresponding numbers for the NLO case in Table 2.
When quark masses are neglected, the only mass scale in the problem is the flow time t, and therefore The numbers may not strictly be minimal; they are to be understood as a reference, in particular in comparison to the NLO numbers given in Table 2.
The numbers may not strictly be minimal; they are to be understood as a reference, in particular in comparison to the NNLO numbers given in Table 1.
Introducing Schwinger parameters as where n ∈ N, the momentum integration reduces to a Gaussian integral: where p = (p 1 , p 2 , p 3 ), and A(s, t) is a coefficient matrix which is linear in the Schwinger parameters s = {s 1 , . . . , s 6 } and the flow-time variables t = {t 0 , . . . , t N }.
Through simple rescaling of the flow-time variables and the Schwinger parameters, one ends up with integrals of the form where M > 0, the P i are polynomials in x 1 , . . . , x M , and the exponents a i can be Ddependent. In the limit 4 − D = 2 → 0, the integrals develop divergences. The integration over the x n can be carried out analytically only for a few simple cases, which is why one needs to resort to numerical integration. 3 This requires the isolation of the terms that become singular as → 0, which can be achieved algorithmically through sector decomposition [20]. In our calculation, we apply this method through the Mathematica package FIESTA [21], which provides us with the result in the form where the ellipsis denotes higher order terms in = (4 − D)/2, and the J n are convergent integrals over rational functions times logarithms of the parameters x 1 , . . . , x M . They can thus be evaluated numerically. We prevented FIESTA from performing this integration, and rather used a fully symmetric integration rule of order 13 [22]. All parts of the integration are performed with high precision arithmetics using the MPFR library. 4 We checked that this algorithm provides us with a reliable estimate of the numerical accuracy.
Let us give the explicit result for one particular non-trivial integral of the type in Eq. (12) which occurs in the calculation of t 2 E(t) . It has four flow-time integrations and thus belongs to the class M 6 (2, 6). Furthermore, from the flow-time integration limits, we see that it originates from the iterated insertion of four 3-point flow-time vertices X (2,0) : The numerical result in the last line is obtained by following the evaluation procedure described above. The numbers in brackets indicate the integration error; for the 1/ 2 -terms we were able to derive an analytical result, for which we simply quote the first few digits of its numerical value. The precision of order 10 −9 as quoted in Eq. (20) for the 1/ 0 -term corresponds to about 250 CPU minutes on an 3 GHz AMD A8 processor; a precision of 10 −6 (10 −4 ) could be achieved within about ten (two) minutes. The CPU time for the 1/ -term is typically several orders of magnitude smaller.

Validation of the calculation
Since this is the first three-loop calculation in the gradient-flow formalism, we considered it of utmost importance to validate our setup. We successfully completed the following checks.
Lower order results. It is important to note that our calculation does not rely on any of the results of Refs. [12,13]. The fact that we reproduced the NLO results evaluated in these papers is therefore an important check of the setup in general. Since the NLO result is known analytically, we can use it also to cross check the numerical accuracy claimed by our integration routine, and we find rather conservative estimates. Specifically, our numerical result agrees with the analytical expression through 10 −15 .
UV-poles at NNLO. The terms of order 1/ 2 and 1/ obtained in our three-loop calculation need to be cancelled by the corresponding terms due to the renormalization of the strong coupling constant at lower orders. We verify this cancellation by analytical integration for the 1/ 2 terms, and numerically through one part in 10 10 for the 1/ terms. Note that the number and complexity of the integrals is typically smaller for higher order poles. However, even though this means that we cannot expect the same numerical accuracy for the finite terms, it should still be sufficient for any foreseeable practical application.
We note in passing that, in the case of the quantity under consideration, the cancellation of the poles is equivalent to the renormalization group (RG) invariance of the final result: where µ is the renormalization scale. The quantity E depends on µ implicitly through α s (µ), and explicitly through terms of the form ln µ 2 t. Knowing the logarithmic dependence in t is thus equivalent to knowing the one in µ. The former is directly obtained from expanding Eq. (14) for → 0, while the latter follows from RG-invariance and can be derived from lower order terms through the perturbative solution of the QCD renormalization group equation: with the first two coefficients of the β function given by 5 where n f is the number of active quark flavors.
Two-loop gluon propagator. As already pointed out above (see the discussion after Eq. (11)), the class M n (2, 2), where in the first term on the r.h.s. of Eq. (9) the flow fields B a µ in Eq. (9) are replaced by their lowest-order terms B a 1,µ , is fully determined by the fundamental gluon self-energy. In fact, using Feynman gauge and adopting the notation of Ref. [12], we may write with the gluon self-energy Using the perturbative expansion of Eq. (25) can be calculated analytically. The coefficientsω i can be taken from the literature. 6 In Feynman gauge, they read where with d given in Eq. (14). We have confirmed the equivalence of both approaches in our setup for some of the most complicated integrals at the level of one part in 10 10 .
Gauge parameter independence. Our setup allows us in principle to perform the calculation for arbitrary gauge parameter λ = 0, see Eq. (2). We have confirmed general λ-independence at NLO, where the number of terms to be evaluated increases by about a factor of ten compared to the case λ = 0. At NNLO, however, the sheer volume of integrals when allowing for general λ makes it impossible to evaluate all of them with meaningful precision in reasonable time. A much more practical though still powerful way is to perform an expansion around λ = 0 and consider only the terms linear in λ. The most significant simplification following from this is that instead of Eq. (4), we obtaiñ In this way, the number of integrals increases again only by a factor of O(10) relative to the case λ = 0. We find gauge parameter independence of the NNLO result for E at O(λ) through 10 −3 for the finite term, and 10 −10 for the 1/ pole terms.

NLO quark-mass effects
Quark loops occur first through the one-loop gluon self-energy, Eq. (26). Quark-mass effects can therefore be taken into account along the lines of Eqs. (25)-(28) by replacing where the sum runs over all active quark flavors q, and the quark-mass (m q ) terms are given by [30] We thus find where The function Ω 1q depends only on 8m 2 q t ≡ m 2 q /q 2 8 ; its numerical size is displayed in Fig. 1. In the limits of small and large quark mass, one finds 3 Results

Action density at three-loop level
We write the result for the vacuum expectation value of the action density as with the NNLO correction factor where α s ≡ α (n f ) s (µ) is the strong coupling renormalized at the scale µ with n f active quark flavors (assumed massless), and N A is the dimension of the adjoint representation of the underlying gauge group (N A = 8 in QCD). Setting µ = 1/ √ 8t, the perturbative coefficients read The NLO coefficient k 1 has been obtained analytically for m q = 0 in Ref. [12]; we add mass effects Ω 1q obtained in Eq. (34). However, for most of our analysis, we find that these terms are numerically irrelevant, and we will neglect them unless stated otherwise. The NNLO coefficient k 2 is the main result of our paper. Similar to Eq. (20), the numbers in brackets denote the numerical uncertainty. The n 2 f -term in k 2 is completely determined by the two-loop gluon propagator, given analytically in Eq. (26). Similar to k 1 , we simply quote the first few digits of its numerical value. Although our main focus is on QCD, we expressed the result of Eqs. (37,38) in terms of "color" factors of a general simple Lie group (see above). For illustration, we also inserted their QCD values and find very well-behaved perturbative coefficients for any realistic value of n f . The expression of t 2 E(t) for general values of the renormalization scale µ is easily reconstructed using Eq. (23). Fig. 4 shows the variation with this unphysical scale for various values of From the input value α  s (see Refs. [31,32] for more details, for example). These start values are then further evolved for fixed n f at the corresponding loop order in order to produce the plots: for the LO/NLO/NNLO result, we apply one/two/three-loop running of α s . Fig. 3 shows the dependence of t 2 E(t) as a function of µ/q 8 for q 8 = 100 GeV and q 8 = 2 GeV. In both cases, one observes a sound perturbative behavior in the interval µ ∈ [q 8 , 3q 8 ]. In addition, the µ-dependence decreases significantly with increasing loop order. These features quickly fade away when going to lower values of µ. Our conclusion is that the best prediction for t 2 E(t) is obtained within the µ-interval [q 8 , 3q 8 ]; its variation within this interval will be used as an estimate of the theoretical uncertainty. Values of µ outside this interval will be disregarded in what follows. Fig. 4 shows t 2 E(t) within this interval for a few values of q 8 ≤ 1 GeV. It is interesting to note that for q 8 = 1/1.5 GeV, corresponding to √ t ≈ 0.1 fm, we may still make quantitative predictions when focussing on the µ-interval identified above. For lower energies, the uncertainty at NNLO becomes of the order of 100%, and the NLO and NNLO correction are of the same order of magnitude.
A common feature of all the plots in Figs. 3 and 4 (except the one at q 8 = 1/1.6 GeV, a value which we will not consider any further in this paper) is that, within µ ∈ [q 8 , 3q 8 ], the maximum is quite precisely at µ = 1.15 q 8 , while the minimum is at µ = 3q 8 . Therefore, the error interval of t 2 E(t) as defined above is given to a very good approximation by its values at µ = µ − ≡ 3q 8 and µ = µ + ≡ 1.15q 8 .    three-loop matching at the quark thresholds), and subsequently at the pertinent order from α s (µ), with µ = 1.15q 8 and µ = 3q 8 for the upper and lower edge of the uncertainty band, respectively. One observes that the resulting NLO and the NNLO bands nicely overlap, which gives confidence in using these bands as measures of the theoretical uncertainty. There is hardly any overlap of these curves with the LO band though.

Extracting α s (m Z )
One of the most interesting applications of our results would be the derivation of a numerical value of α s (m Z ) ≡ α (5) s (m Z ) using lattice data as input. This will be most promising, of course, if the lattice calculation for t 2 E(t) could be extended to the perturbative regime, which seems to have become a realistic perspective [33].
Assume that a lattice value e(t) for t 2 E(t) is known, evaluated at t = 1/(8q 2 8 ) and for n f active quark flavors. Using the perturbative result of Eqs. (36) and (37) through order l (including its µ-dependence), one can derive an l-loop value for α (n f ) s (µ), which can then be converted into a value for α (5) s (m Z ) through four-loop RG evolution and three-loop matching to the n f = 5 theory. Table 3 Table 3: Numerical values for 10 4 · t 2 E(t) corresponding to various α s (m Z ) ≡ α (5) s (m Z ). Given a numerical result for t 2 E(t) (e.g., from a lattice calculation), this table lets one deduce the corresponding value of α s (m Z ). The associated perturbative uncertainty for n f = 5 and n f = 3 can be read off from Fig. 6. the center of the error band, i.e., they are the arithmetic means of t 2 E(t) evaluated at µ = 1.15 q 8 and µ = 3 q 8 . These numbers take into account the NLO quark effects given in Eq. (34), whereupon the lightest three quark flavors are taken massless, while m c = 1.67 GeV and m b = 4.78 GeV. The mass effects therefore only affect the columns with n f ≥ 4. At q 8 = 2 GeV, their effect on t 2 E(t) is about 0.8%, at q 8 = 10 GeV it is less than 0.3% both for n f = 4 and n f = 5, while at q 8 = m Z , they have no effect on the digits given in the table.
In accordance with our previous considerations, we estimate the theoretical accuracy of this extraction by considering t 2 E(t) at µ = 1.15 q 8 and 3 q 8 when deriving α The result for n f = 3 is shown in Fig. 6. In lack of a precise value of e(t) at sufficiently large q 8 , we substitute it by the perturbative NNLO expression for t 2 E(t) at µ = q 8 , where the numerical value for α The lower part of the figure shows the theoretical accuracy that could be achieved by such an analysis, derived by taking the relative width of the bands of the upper part of the plot, For example, if e(t) is given only at t = 1/(8GeV 2 ), the NNLO uncertainty on α s (m Z ) would be around 2.5%. On the other hand, knowning e(t) at t = 1/(8m 2 Z ) would allow one to derive α s (m Z ) to 0.5% accuracy which is at the same level as the current world average on this quantity [34]. Also shown in the lower plot is the uncertainty which results from knowing e(t) for n f = 5 active flavors (lower dotted red line). In this case, the numbers above decrease to ∼ 1.1% and ∼ 0.3%, respectively, because of the lower value of the QCD β function.
s (m Z ) derived at LO (gray), NLO (orange), and NNLO (red) from a hypothetical exact value of t 2 E(t) | n f =3 (see main text for details). Lower plot: corresponding theoretical uncertainty (see Eq. (40)). The red dotted line in the lower plot shows the uncertainty when the analysis is based on t 2 E(t) | n f =5 .
Again, the µ-dependent terms can be easily reconstructed using renormalization group invariance.
Performing a similar analysis for W (t) as done in the preceeding sections for t 2 E(t) , we see no improvement concerning the precision for the extraction of α s relative to the one based on t 2 E(t) .

Conclusions and Outlook
The action density for QCD gradient flow fields has been evaluated at three-loop level. The perturbative expansion has been derived by standard Wick contractions, and the resulting integrals have been solved by sector decomposition supplied by a suitable numerical integration algorithm. A number of strong checks on the result has been performed. In addition, quark-mass effect have been included at NLO.
Our NNLO coefficient indicates a very well-behaved perturbative series for the action density down to energy scales of about q 8 ∼ 0.65 GeV, corresponding to √ t ∼ 0.11 fm. This seems well within reach of a direct comparison to a lattice evaluation of t 2 E(t) . Given that t 2 E(t) can be evaluated independently (e.g. by a lattice calculation) at sufficiently large values of the flow time with high precision, one may derive a numerical value for α s (m Z ) by comparison to the perturbative result. We provide an estimate of the resulting uncertainty and find that it could be competitive with the current world average.
On the perturbative side, further steps could be the development of more efficient tools for the evaluation of the integrals, the consideration of other observables, or the application of the flow-field formalism to quark fields as introduced in Ref. [35].
Finally, it should be noted that there is no conceptual limitation of the calculational method described in this paper which would restrict it to the three-loop level. In the current implementation, however, an extension to four loops would require a significant increase in the computing resources.