Four-loop QCD MOM beta functions from the three-loop vertices at the symmetric point

For the first time, we compute three-loop contributions to all triple vertices in QCD at the symmetric point. The analytic results are obtained in massless QCD with an arbitrary color group in the Landau gauge. All new loop integrals are expressed in terms of harmonic polylogarithms at the sixth root of unity. These corrections allow us to derive expressions for the four-loop QCD beta function in a set of momentum-subtraction schemes.


INTRODUCTION
The strong coupling constant α s ≡ g 2 s /(4π) is a fundamental parameter of QCD. It enters the predictions of many Standard Model (SM) observables, e.g., the Higgs production cross-section [1] at the LHC, and the uncertainty of α s significantly influences the comparison between theory and experiment. In itself, α s is not a physical observable, its value depends on the scale Q 2 and the utilized renormalization prescription (or scheme). The freedom in choosing a normalization point manifests itself in "running" governed by the renormalization group (RG) equations dα s (Q 2 ) d ln Q 2 = β(α s (Q 2 )). (1) Given the renormalization scheme (RS), one can calculate the β-function order by order in perturbation theory (PT). To parametrize the strength of strong interactions, one usually uses the MSscheme and quotes the value of α (5) s (M Z ) = 0.1179 (10) [2] defined in effective QCD with n f = 5 massless flavors and evaluated at the Z-boson mass. This definition is of great convenience when studying inclusive observables dominated by short-distance effects. For almost 20 years the four-loop β MS [3,4] in the MS scheme has been utilized in state-of-the-art QCD calculations, and only recently the five-loop result [5][6][7][8] became available in the literature. The latter allows one to significantly improve the agreement between α s determinations (see, e.g., Ref. [2] for details) in precision measurements carried out at different energy scales.
At lower energies other α s definitions can be more convenient. For example, in lattice QCD one can introduce the strong coupling as a particular RG-invariant combination (so-called invariant charge) of the vertex and two-point Green functions evaluated at fixed external momenta in a fixed (usually Landau) gauge (see, e.g., Refs. [9,10]). Since lattice results [11] give rise to one of the most precise determinations of α (5) s (M Z ), it is important to study the strong coupling and its running in the momentum-subtraction (MOM) schemes. In these RSs the vertex functions and propagators are normalized in such a way that at a certain kinematic point there are no corrections beyond the tree level. As a consequence, the renormalized couplings coincide with the corresponding invariant charges and can be directly compared to the nonperturbative lattice results. FIG. 1. For symmetric point integrals p 2 1 = p 2 2 = q 2 and for auxiliary integrals we use p 2 1 = p 2 2 , q 2 = x p 2 1 Contrary to the MS scheme, in the MOM scheme one needs to know the Green functions beyond the divergent terms. Some choices of kinematics can make the calculation quite challenging , especially of vertices. In this paper, we consider three-point functions that depend on momenta p 1 , p 2 , and q = p 1 + p 2 (see Fig. 1) and utilize the symmetric normalization point (SMOM) in the Euclidean region p 2 1 = p 2 2 = q 2 = −Q 2 . One-loop QCD renormalization in the SMOM scheme in a general linear gauge has been known for quite a long time, following the pioneering work arXiv:2002.02875v2 [hep-ph] 10 Apr 2020 in Ref. [12]. At the moment, only two-loop expressions for the QCD vertices in this kinematics are available in the literature in numerical [13] and analytical [14] forms, allowing one to find three-loop [13,15] strong-coupling beta functions in SMOM RS. In our work, we exploit modern Feynman-integral evaluation methods to improve these results by one more order of PT. To simplify our calculation, we routinely use the Landau gauge, in which the gluon propagator is transversal and there is no need to consider the gauge-parameter renormalization.
In the context of perturbative QCD (pQCD) one can relate renormalized couplings a ren ≡ α ren s /(4π) defined in different RSs. In what follows, we consider a ren = {a ggg , a ccg , a qqg }, which satisfy with a bare being the bare coupling defined in the dimensionally regularized theory with d = 4 − 2ε. The bare coupling is divergent and is related to the renormalized ones a ren via divergent factors Z aren . The latter are combinations of the renormalization constants of the three-gluon (Z ggg ), ghost-gluon (Z ccg ), and quark-gluon (Z qqg ) vertices together with those of the gluon (Z gg ), ghost (Z cc ) and quark (Z qq ) fields. The vertex renormalization constants denoted collectively by Z V are obtained order by order in perturbation theory from the bare Γ V bare vertices by making sure that the renormalized counterpart Γ V is finite for all external momenta q i and satisfies certain normalization conditions. In the MSscheme only divergent terms are subtracted, and one can prove that a MS = a MS ggg = a MS ccg = a MS qqg . On the contrary, in the SMOM schemes one also subtracts finite terms and requires that there are no O(a ren ) corrections to a particular Γ V ren at the symmetric point, characterized by momentum Q 2 . To avoid the appearance of log(Q 2 /µ 2 ) in the renormalization constants, it is convenient to choose Q 2 = µ 2 .
The vertex functions of our interest have color and space-time indices and one can decompose them in terms of basis tensors with Lorentz-invariant coefficients (form-factors). The choice of the basis is not unique and we make use of a decomposition [14] taht is valid for the symmetric point. To save space, we present here only terms [13] that are relevant for the definition of the required vertex renormalization constants (3): Here all momenta are assumed to be outgoing, p 1 in Eq.(4) corresponds to the antighost, and T µνρ = g µν (p 1 −p 2 ) ρ +. . . represents the tensor that enters into the tree-level three-gluon vertex. The SU(N) generators T a ij in the quark-gluon vertex (6) ij with structure constants f abc . The expressions for the bare form factors Γ V bare with V = {ggg, ccg, qqg} are extracted from the tensor vertices by means of the projectors given in Ref. [14]. To define a ggg in the respective SMOM scheme (MOMggg) via Eq. (2), we require that at the symmetric point Γ ggg ren (−µ 2 ) = 1, i.e. Z −1 ggg = Γ ggg bare (−µ 2 ). In the same way one can relate the bare coupling a bare to a ccg or a qqg by requiring that either Γ ccg ren (−µ 2 ) = 1 (MOMh), or Γ qqg ren (−µ 2 ) = 1 (MOMq). Since field renormalization constants also enter into Eq.(2), we have to impose conditions on Z gg , Z cc , and Z qq in the SMOM scheme: there should be no corrections to the corresponding tree-level propagator for the external momentum q 2 = −µ 2 .
From Eq.
(2) one can deduce that two renormalization prescriptions for a ren , say a R and a MS , are related via finite correction factors X R : Given X R at L loops, one can determine the (L + 1) MOM-scheme beta functions from β MS via a relation that is valid in the Landau gauge, where in the final step we invert Eq. (7) to express β R in terms of a R . The main aim of this paper is to calculate the three-loop corrections X R to the relations (7) between a ggg , a ccg , a qqg and a MS . As one of the applications of our result, we use the four-loop beta function in the MSscheme [3, 4] to find the corresponding beta functions β ggg , β ccg , and β qqg in the considered SMOM schemes.

DETAILS OF THE CALCULATION
We generate Feynman diagrams with DIANA [16], and obtain 8, 106, and 2382 graphs for the Γ ggg ren three-gluon vertex at one, two, and three loops, respectively. Both ghost-gluon and quark-gluon vertices give rise to 2, 33, and 688 diagrams at the same loop levels. After the application of projectors [14] and taking fermion and color traces [17], we are left with scalar Feynman integrals, which we reduce to a set of two one-loop, eight two-loop, and 51 three-loop master integrals by means of Reduze 2 [18] and FIRE6 [19].
The main challenge is calculating of the full set of three-loop three-point integrals in SMOM kinematics, which was not available in the literature. To compute master integrals, we rely on the linear reducibility of massless vertex-type integrals with arbitrary off-shell momenta, which is proven to take place up to three-loop order [20,21]. Due to the latter property, for our integrals in more restricted kinematics we choose two strategies of evaluation. First of all, we derive a new basis of finite master integrals [22] and try to compute them directly by means of the HyperInt package [23] in terms of generalized polylogarithms (GPLs). Unfortunately, we are unable to calculate the most complicated integrals and instead use another strategy based on the solution of the system of differential equations (DEs). The method of DEs cannot be applied directly to the calculation of single scale integrals we are interested in, and we construct a set of auxiliary integrals with arbitrary external q 2 = x p 2 1 (Fig.1). After switching to a new variable z, x = 2 − z − 1/z, we reduce the original DE system to the so called ε-form [24]. Due to the presence of singularities at complex points, we make use of the epsilon package [25], which is capable of dealing with the latter. The obtained system in the ε-form is easily solved order by order in the ε expansion. The solution is given by linear combinations of GPLs G(a 1 , . . . , a n ; z), where a i correspond to different sixth roots of unity, with a number of unknown constants to be fixed from boundary conditions. The latter can be obtained by matching the expansion in the limit q 2 → 0, corresponding to z → 1, with the explicit result of large-momentum asymptotic expansion. Since z = 1 is a singular point of the DEs, naive Taylor expansion is not sufficient to fix all of the constants and we exploit the EXP package [26,27] to generate the series in terms of massless propagators. We compute the latter by the MINCER package [28,29] keeping exact dependence 1 on the space-time dimension variable.
In this way, we obtain the analytic results for auxiliary integrals depending on z. Taking the limit z → e iπ/3 , which is regular and corresponds to SMOM kinematics, we compute the required single-scale master integrals. The correctness of the analytic calculation is verified numerically using pySecDec [30]. Substituting the integrals in the expressions for the three-loop Γ V bare and expanding in ε up to the necessary order, we see that the maximal transcendental weight in the final result is, as expected, 2L = 6 with L being the loop order. Using the basis and the reduction rules from Ref. [31], we are able to simplify the expressions for Γ V bare significantly. As a consequence, we obtain the SMOM renormalization constants, conversion factors (7), and the beta functions (8) in a rather compact form.

RESULTS AND CONCLUSION
All of the necessary renormalization constants up to three-loop order were calculated iteratively via Eq. (3). We also reproduced the well-known three-loop MS expressions, required to derive the relations (7) between the SMOM couplings and a MS evaluated at the same scale Q 2 = µ 2 . To save space, we present the results in numerical form with all QCD color factors substituted explicitly: The analytic results for the general gauge group are available as a Supplementary Material and have several remarkable properties. We managed to simplify the two-loop part [15] and expressed it solely in terms of powers of π, odd ζ values, and two polygamma functions ψ (1) (1/3) and ψ (3) (1/3). It turns out that to write down the three-loop contribution we only need to introduce three additional constants: ψ (5) (1/3), and two combinations H 5 and H 6 of GPLs with uniform transcendental weights five and six, respectively. Evaluating H 5,6 numerically [32] with high precision and using the PSLQ algorithm [33], we reconstructed the new constants 2 through a more restricted basis [35] of real parts of harmonic polylogarithms of the argument e iπ/3 .
As a first application of Eqs.(9),(10), and (11) we obtain a set of SMOM beta functions via Eq. (8), thus extending the results of Ref. [15] to four-loop order: with β uni (a) = −a 2 (11 − 2/3 n f ) − a 3 (102 − 38/3 n f ) being the universal two-loop scheme-independent contribution. The obtained formulas can be used in a number of ways. For example, one can improve the precision of matching between the lattice and pQCD results. In addition, the possibility to switch from one RS to another in truncated PT series for observables, or more generally RG-invariant quantities (see, e.g. Refs. [36,37]), provides us with an additional handle on theoretical uncertainties beyond simple scale variation. Moreover, the computed integrals can also be used in SMOM calculations of more complicated three-point Green functions with operator insertions, e.g., for studies of the light-quark masses as in Refs. [38,39]. We