The five-loop beta function of Yang-Mills theory with fermions

We have computed the five-loop corrections to the scale dependence of the renormalized coupling constant for Quantum Chromodynamics (QCD), its generalization to non-Abelian gauge theories with a simple compact Lie group, and for Quantum Electrodynamics (QED). Our analytical result, obtained using the background field method, infrared rearrangement via a new diagram-by-diagram implementation of the R* operation and the Forcer program for massless four-loop propagators, confirms the QCD and QED results obtained by only one group before. The numerical size of the five-loop corrections is briefly discussed in the standard MSbar scheme for QCD with n_f flavours and for pure SU(N) Yang-Mills theory. Their effect in QCD is much smaller than the four-loop contributions, even at rather low scales.


Introduction
The scale dependence ('running') of the renormalized coupling constant α i is a fundamental property of an interacting quantum field theory. In renormalization-group improved perturbation theory, the beta function governing this dependence can be written as da d ln µ 2 = β(a) = − ∞ n=0 β n a n+2 , a = α i (µ) 4π (1.1) where µ is the renormalization scale. The determination of the (sign of the) leading one-loop coefficient β 0 [1][2][3][4][5], soon followed by the calculation of the two-loop correction β 1 [6,7], led to the discovery of the asymptotic freedom of non-Abelian gauge theories and thus paved the way for establishing QCD as the theory of the strong interaction. The renormalizationscheme dependent three-loop (next-to-next-to-leading order, N 2 LO) and four-loop (nextto-next-to-next-to-leading order, N 3 LO) coefficients β 2 and β 3 were computed in refs. [9,10] and [11,12] in minimal subtraction schemes [13,14] of dimensional regularization [15,16].
In the past years, the N 2 LO accuracy has been reached for many processes at highenergy colliders. N 3 LO corrections have been determined for structure functions in inclusive deep-inelastic scattering (DIS) [17,18] and for the total cross section for Higgs-boson production at hadron colliders [19,20]. Some moments of coefficient functions for DIS have recently been computed at N 4 LO [21]. Reaching this order would virtually remove the uncertainty due to the truncation of the series of massless perturbative QCD in determinations of the strong coupling constant α s from the scaling violations of structure functions in DIS.
The corresponding five-loop contributions to the beta functions of QCD, with all colour factors 'hard-wired', and QED have already been computed in refs. [22,23]. Their leading large-n f contributions have long been known [24], and the sub-leading large-n f terms have been checked and generalized to a general simple gauge group in ref. [25]. The real tour de force of ref. [22] though, are the parts proportional to n 0 f , n 1 f and n 2 f which together JHEP02(2017)090 required more than a year of computations on a decent number of multi-core workstations in a highly non-trivial theoretical framework. These critical parts have neither been extended to a general gauge group nor validated by a second independent calculation so far.
In the present article we address this issue and present the five-loop beta function for a general simple gauge group. Unlike the calculations in refs. [4][5][6][7][8][9][10][11][12], we have employed the background field method [26,27], which we found to be more efficient -in validation calculations of the Forcer program [28][29][30] of the four-loop renormalization of Yang-Mills theories to all powers of the gauge parameter -than the computation of two propagators and a corresponding vertex. This method and other theoretical and calculational issues, in particular a new implementation [31] of the R * operation [32][33][34][35] for massless propagatorlike diagrams, are addressed in section 2; the details of the required tensor reduction can be found in the appendix. We present and discuss our result in section 3, and briefly summarize our findings in section 4.

Theoretical framework and calculations
In this section we briefly review the background-field formalism and the R * operation. We further define our notations for group invariants, and we give an overview of our calculation.

Background field method
A convenient and efficient method to extract the Yang-Mills beta function is to make use of the background field. We will briefly review this formalism. A convenient starting point is the Lagrangian of Yang-Mills theory coupled to fermions in a non-trivial (often the fundamental) representation of the gauge group, the theory for which we will present the 5-loop beta-function in the next section.
The Lagrangian of this theory can be decomposed as Here the classical Yang-Mills Lagrangian (CYM), a gauge-fixing term (GF), the Faddeev-Popov ghost term (FPG) and the fermion term (FER) are given by In the fermion term the sum goes over colours i, j, and n f flavours f , and we use the standard Feynman-slash notation. The field strength is given by and the covariant derivatives are defined as

JHEP02(2017)090
The conventions associated to the generators T a and structure constants f abc of the gauge group will be explained in section 2.2. The gauge-fixing term depends on making a suitable choice for G a , which is usually taken as G a = ∂ µ A a µ . The background-field Lagrangian is derived by decomposing the gauge field as where B a µ (x) is the classical background field whileÂ a µ (x) contains the quantum degrees of freedom of the gauge field A a µ (x). The background-field Lagrangian is then written as L BYM+FER = L BCYM + L BGF + L BFPG + L BFER . (2.6) L BCYM and L BFER are derived simply by substituting eq. (2.5) into the corresponding terms in the Yang-Mills Lagrangian. However a clever choice exists [26,27] for the ghost and gauge fixing terms, which allows this Lagrangian to maintain explicit gauge invariance for the background field B a µ (x), while fixing only the gauge freedom of the quantum field A a µ (x). The gauge fixing then uses instead while the ghost term is given by The Lagrangian L BYM+FER then gives rise to additional interactions which are different from the normal QCD interactions of the quantum fieldÂ a µ (x) also contain interactions of B a µ (x) with all other fields. A remarkable fact is found when considering the renormalization of this Lagrangian. Indeed it turns out, see e.g., [26,27], that the coupling renormalization, g → Z g g, which determines the beta function, is directly related to the renormalization of the background field, B → BZ B , via the identity: When working in the Landau gauge, the only anomalous dimension needed in the background field gauge formalism is then the beta function. However in the Feynman gauge the gauge parameter ξ requires the renormalization constant Z ξ -which equals the gluon field renormalization constant -but only to one loop lower. In turn this allows one to extract the beta function from the single equation where Π µν B (Q 2 ; ξ, g) is the bare self energy of the background field. This self-energy is computed by keeping the fields B external while the only propagating fields areÂ, η and ψ. A typical diagram which contributes to Π B (Q 2 ; ξ, g) is given in figure 1.
Obtaining the beta function through the background field gauge is faster and simpler than the traditional method of computing the gluon propagator, ghost propagator and ghost-ghost-gluon vertex due to a lower total number of diagrams and the above reduction to a scalar renormalization.

Group notations
In this section we introduce our notations for the group invariants appearing in the results of the next section. T a are the generators of the representation of the fermions, and f abc are the structure constants of the Lie algebra of a compact simple Lie group, (2.12) The quadratic Casimir operators C F and C A of the N -dimensional fermion and the N Adimensional adjoint representation are given by [T a T a ] ik = C F δ ik and f acd f bcd = C A δ ab , respectively. The trace normalization of the fermion representation is Tr( At L ≥ 4 loops also quartic group invariants enter the beta function. These can be expressed in terms of contractions of the totally symmetric tensors Here the matrices [ C a ] bc = −if abc are the generators of the adjoint representation. It should be noted that in QCD-like theories without particles that are colour neutral, Furry's theorem [36] prevents the occurrence of symmetric tensors with an odd number of indices. For the fermions transforming according to the fundamental representation and the standard normalization of the SU(N ) generators, these 'colour factors' have the values (2.14) The results for QED (i.e., the group U(1)) are obtained for C A = 0, d abcd For a discussion of other gauge groups the reader is referred to ref. [11].

The R * -operation
As outlined above, it is possible to extract the five-loop beta function from the poles (in the dimensional regulator ǫ) of the bare background field self-energy Π B (Q). At present JHEP02(2017)090 p Figure 2. One external line is moved to create a topology that can be integrated. Here we do this for the diagram of figure 1. One should take into account that there can be up to 5 powers of dot products in the numerator, causing many subdivergences. Furthermore, the double propagator that remains on the right can introduce infrared divergences. After the subdivergences have been subtracted, the integral over p can be performed and the remaining four-loop topology can be handled by the Forcer program.
it is beyond current computational capabilities to calculate the required five-loop propagator integrals directly. The main obstacle preventing such an attempt is the difficulty of performing the required integration-by-parts (IBP) reductions.
Fortunately the problem can be simplified via the use of the R * -operation. The R *operation [32][33][34][35] is a subtraction operation capable of rendering any propagator integral finite by adding to it a number of suitable subtraction terms. The subtraction terms are built from potentially high rank tensor subgraphs of the complete graph, whose tensor reduction requires involved methods which we present in appendix A. Via the procedure of IR-rearrangement, these subtraction terms can subsequently be related to simpler propagator integrals. The IR-rearranged integral is, in general, any other propagator integral obtained from the original one by rerouting the external momentum in the diagram. This is illustrated in figure 2.
For integrals whose superficial degree of divergence (SDD) is higher than logarithmic, the SDD is reduced by differentiating it sufficiently many times with respect to its external momenta, before IR-rearranging it.
The upshot of this procedure is that the IR-rearranged propagator integrals can be chosen to be carpet integrals, which correspond to graphs where the external lines are connected only by a single propagator. A carpet integral of L loops can be evaluated as a product of an (L − 1) loop tensor propagator integral times a known one-loop tensor integral. In the case of the five-loop beta function this means that we can effectively evaluate the poles of all five-loop propagator integrals from the knowledge of propagator integrals with no more than four loops. A sketch of the R * -operation to compute the superficial divergence of a 3-loop diagram is shown below:

JHEP02(2017)090
where sup denotes the superficial divergence, and K isolates the pole of a Laurent series in ǫ. As can be seen, the R * -operation is recursive, since the same procedure needs to be applied to compute the superficial divergence of each counterterm. The Forcer program [29,30], written in the Form language, is capable to efficiently compute the subtraction terms. It reduces four-loop propagator integrals to simpler known ones by integrating two-point functions, and by applying parametrically solved IBP reduction rules to eliminate propagators. We have automated the R * -operation in a fast Form program, capable of performing the subtraction of propagator integrals with arbitrary tensorial rank. Having interfaced the Forcer program with the R * program we were able to compute the poles of all integrals entering the five-loop background field self-energy. The algorithms and details of our implementation of the R * -operation follow to some degree the ideas which were presented in the literature (see e.g., [32-35, 37, 38]), however we have generalized certain notions in order to deal with arbitrary tensor integrals and their associated ultraviolet and infrared divergences. These generalizations are subtle and will be presented elsewhere [31].

Diagram computations and analysis
The Feynman diagrams for the background propagator up to five loops have been generated using QGRAF [39]. They have then been heavily manipulated by a Form [40][41][42] program that determines the topology and calculates the colour factor using the program of ref. [43]. Additionally, it merges diagrams of the same topology, colour factor, and maximal power of n f into meta diagrams for computational efficiency. Integrals containing massless tadpoles or symmetric colour tensors with an odd number of indices have been filtered out from the beginning. Lower-order self-energy insertions have been treated as described in ref. [44]. In this manner we arrive at 2 one-loop, 9 two-loop, 55 three-loop, 572 four-loop and 9414 five-loop meta diagrams.
The diagrams up to four loops have been computed earlier to all powers of the gauge parameter using the Forcer program [28][29][30]. For the time being, our five-loop computation has been restricted to the Feynman gauge, ξ F = 1 − ξ = 0. An extension to the first power in ξ F would be considerably slower; the five-loop computation for a general ξ would be impossible without substantial further optimizations of our code. Instead of by varying ξ, we have checked our computations by verifying the relation Q µ Q ν Π µν B = 0 required by eq. (2.11). This check took considerably more time than the actual determination of β 4 .
The five-loop diagrams have been calculated on computers with a combined total of more than 500 cores, 80% of which are older and slower by a factor of almost three than the latest workstations. One core of the latter performs a 'raw-speed' Form benchmark, a fourdimensional trace of 14 Dirac matrices, in about 0.02 seconds which corresponds to 50 'form units' (fu) per hour. The total CPU time for the five-loop diagrams was 3.8 · 10 7 seconds which corresponds to about 2.6 · 10 5 fu on the computers used. The TForm parallelization efficiency for single meta diagrams run with 8 or 16 cores was roughly 0.5; the whole calculation of β 4 , distributed 'by hand' over the available machines, finished in three days.
For comparison, the corresponding R * computation for ξ F = 0 at four loops required about 10 3 fu, which is roughly the same as for the first computation of the four-loop beta JHEP02(2017)090 function to order ξ 1 F by a totally different method in ref. [11]. The computation with the Forcer program at four and fewer loops is much faster, in fact fast enough to comfortably demonstrate the full three-loop renormalization of QCD in 10 minutes on a laptop during a seminar talk [45].
The determination of Z B from the unrenormalized background propagator is performed by imposing, order by order, the finiteness of its renormalized counterpart. The beta function can simply be read off from the 1/ε coefficients of Z B . If the calculation is performed in the Landau gauge, the gauge parameter does not have to be renormalized. In a k-th order expansion about the Feynman gauge at five loops, the L < 5 loop contributions are needed up to ξ 5−L F . The four-loop renormalization constant for the gauge parameter is not determined in the background field and has to be 'imported'. In the present k = 0 case, the terms already specified in ref. [12] would have been sufficient had we not performed the four-loop calculation to all powers of ξ F anyway.

Results and discussion
Before we present our new result, it may be convenient to recall the beta function (1.1) up to four loops [4][5][6][7][8][9][10][11][12] in terms of the colour factors defined in section 2, where n f is the number of fermion (in QCD, quark) flavours. β n are the same in all MS-like schemes [13,14], i.e. within the class of renormalization schemes which differ only by a shift

JHEP02(2017)090
of the scale µ. In the same notation and scheme, the five-loop contribution reads

JHEP02(2017)090
ζ denotes the Riemann zeta function with ζ 3 ∼ = 1.202056903, ζ 4 = π 4 /90 ∼ = 1.08232323 and ζ 5 ∼ = 1.036927755. As expected from the lower-order and QED results, higher values of the zeta function do not occur despite their occurrence in the results for individual diagrams; for further discussions see refs. [23,46]. Inserting the group factors of SU(3) as given in eq. (2.14) leads to the QCD results In truncated numerical form β 3 and β 4 are given by

JHEP02(2017)090
The (corresponding parts of the) results (3.5), (3.7) and (3.11) are in complete agreement with the findings of refs. [22][23][24][25]. Consequently, eq. (3.11) also agrees with the result for QED at n f = 1, which was obtained in ref. [47] somewhat earlier than the general result [23]. As already noted in ref. [22], the five-loop QCD coefficient of the beta function is rather small [ recall that we use a convenient but very small expansion parameter in eq. (1. where β ≡ −β(a s )/(a 2 s β 0 ) has been re-expanded in powers of α s = 4π a s . Clearly there is no sign so far of a possible divergence of the perturbation series for this quantity.
In order to further illustrate the n f -dependent convergence (or the lack thereof) of the beta function of QCD, we introduce the quantity Recalling the normalization (1.1) of our expansion parameter, α s (n f ) represents the value of α s for which the n-th order correction is 1/4 of that of the previous order. Therefore, s (n f ) defines (somewhat arbitrarily due to the choice of a factor of 1/4) a region of fast convergence of β(α s , n f ). Obviously, the absolute size of the n-th and (n−1)-th order effects are equal for α s = 4 α (n) (n f ). Thus the quantity (3.13) also indicates where the expansion appears not to be reliable anymore, α s > ∼ 4 α (n) s (n f ), for a given value of n f that is not too close to zeros or minima of the coefficients β n−1 and β n .
It is interesting to briefly study the N -dependence of the convergence behaviour for the case of SU(N ) gauge theories. For our brief illustration we confine ourselves to pure Yang-Mills theory, n f = 0, and consider α (n) , (3.14) where the factor N compensates the leading large-N dependence N n+1 of β n , i.e., the parameter that needs to be small in SU(N ) Yang-Mills theory is not α YM but N α YM . The quantities (3.13) and (3.14) are displayed in the left and right panel of figure 3, respectively. The behaviour of α (n) s at the upper end of the n f range shown in the figure is affected by the zeros and minima of the coefficients β n > 0 mentioned below eq. (3.9). The N -dependence of α YM for pure Yang-Mills theory, where only terms with N n+1 and N n−1 enter β n (the latter only at n ≥ 4 via d abcd A d abcd A /N A , cf. eq. (2.14) above), is rather weak. With only the curves up to four loops, one might be tempted to draw conclusions from the shrinking of the 'stable' α s region from NLO to N 2 LO and from N 2 LO to N 3 LO that are not supported by the N 4 LO (five-loop) results of ref. [22] and the present article.   figure 4 are as small as 0.08% (0.4%) at µ 2 = 3 GeV 2 (1 GeV 2 ); the corresponding N 3 LO corrections are 0.5% (2%). Of course these results do not preclude sizeable purely non-perturbative corrections, but it appears that the perturbative running of α s is now fully under control for all practical purposes.

Summary and outlook
The five-loop (next-to-next-to-next-to-next-to-leading order, N 4 LO) coefficient β 4 of the renormalization-group beta function has been computed in MS-like schemes for Yang-Mills theories with a simple compact Lie group and one set of n f spin-1/2 fermions. This computation confirms and extends the QCD and QED results first obtained, respectively, JHEP02(2017)090 in ref. [22] -where also some direct phenomenological applications to α s determinations from, e.g., τ -lepton decays and Higgs-boson decay have already been discussed -and ref. [23]. It also agrees with the high-n f partial results of refs. [24,25]. We have illustrated the size of the resulting N 4 LO corrections to the scale dependence of the coupling constant for α s -values relevant to MS, the default scheme for higher-order calculations and analyses in perturbative QCD. For physical values of n f , the N 4 LO corrections to the beta function are much smaller than the N 3 LO contributions and amount to 1% or less, even for α s -values as large as 0.4. More generally, there is no evidence of any increase of the coefficients indicative of a non-convergent perturbative expansion for the beta functions of QCD and SU(N ) gauge theories. Our computation has been made possible by the development of a refined algorithm [31], implemented in Form [40][41][42], for the determination of the ultraviolet and infrared divergences of arbitrary tensor self-energy integrals via the R * operation [32][33][34][35] -for another recent diagrammatic implementation of R * for scalar integrals and its application to ϕ 4 theory at six loops, see refs. [48,49] -and the Forcer program [28][29][30] for the parametric reduction of four-loop self-energy integrals. It should be noted that this approach is quite different from those taken in refs. [22] and [25]. In the former the R * operation has been carried out 'globally', the latter uses a five-loop extension of the method of fully massive vacuum diagrams as applied for the determination of the four-loop beta function in refs. [11,12]; see also ref. [50].

JHEP02(2017)090
One may expect that the present implementation of the R * operation will be useful for other multi-loop calculations, at least after further optimizations. An example is the computation of the fifth-order contributions to the anomalous dimensions of twist-2 spin-N operators in the light-cone operator product expansion, which now represent the only missing piece for full N 4 LO analyses of low-N moments of the structure functions F 2 and F 3 in inclusive deep-inelastic scattering.
A Form file with our result for the coefficient β 4 and its lower-order counterparts can be obtained from the preprint server http://arXiv.org by downloading the source of this article. It will also be available from the authors upon request.

JHEP02(2017)090
metric tensors can be utilized to reduce the size of the system. From eq. (A.3) it follows that the projector P σ is in the same symmetry group (the group of permutations which leave it invariant) as T σ . For example, given a permutation σ 1 = (123 . . . (r − 1)r), T µ 1 ... µr σ 1 = g µ 1 µ 2 g µ 3 µ 4 . . . g µ r−1 µr . (A.4) The corresponding projector P µ 1 ... µr σ 1 must be symmetric under interchanges of indices such as µ 1 ↔ µ 2 , (µ 1 , µ 2 ) ↔ (µ 3 , µ 4 ) and so on. Grouping the metric tensors by the symmetry leads to the fact that P σ is actually written in a linear combination of a small number of m tensors instead of n (m ≤ n), (A.5) The m sets of permutations A σ k=1...m must therefore each be closed under the permutations which leaves T σ invariant and at the same time their union must cover once the set 2 S n . Contracting P σ with T τ s where we choose a representative permutation τ from each A σ k , i.e one permutation from A σ 1 , one permutation from A σ 2 etc, gives an m × m matrix which can be inverted to yield the coefficients b k . The number of unknowns m is, for example m = 5 for r = 8 and m = 22 for r = 16, which are compared to n = 105 for r = 8 and n = 2027025 for r = 16. The comparison of these numbers illustrates that the exploitation of the symmetry of the projectors makes it possible to find the tensor reduction even for very large values of r, which could never have been obtained by solving the n × n matrix.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.