On Higgs decays to hadrons and the R-ratio at N^4LO

We present the first determination of Higgs-boson decay to hadrons at the next-to-next-to-next-to-next-to-leading order of perturbative QCD in the limit of a heavy top quark and massless light flavours. This result has been obtained by computing the absorptive parts of the relevant five-loop self-energy for a general gauge group and combining the outcome with the corresponding coefficient function already known to this order in QCD. Our new result reduces the uncertainty due to the truncation of the perturbation series to a fraction of the uncertainty due to the present error of the strong coupling constant. We have also performed the corresponding but technically simpler computations for direct Higgs decay to bottom quarks and for the electromagnetic R-ratio in e^+ e^- ->hadrons, thus verifying important fifth-order results obtained so far only by one group.


Introduction
The production and decay processes of the Higgs boson, discovered five years ago at CERN [1,2] with a mass M H of 125 GeV, are among the most important research topics in collider physics. The dominant standard-model decay is that to bottom quarks, H →bb (+ hadrons). The QCD calculations of this decay mode have been completed up to the fourth order in the strong coupling α s , see ref. [3] and references therein. A crucial component of this high accuracy is the next-to-next-to-next-to-next-to-leading order (N 4 LO) computation [4] of the decay to quarks via their direct coupling to the Higgs. This calculation, in which the quark mass can be neglected (except in the Yukawa coupling), has not been repeated so far.
The second important hadronic decay channel arises via H → gg, where the coupling of the Higgs to gluons is predominantly mediated, in the standard model, by the top quark. Due to M H ≪ 2m t , high-order QCD corrections to this process can be evaluated in an effective theory in which the top quark has been integrated out [5]; for 1/m t corrections up to NNLO see refs. [6,7]. The resulting coefficient function for the effective Higgs coupling to gluons is known to N 4 LO [8][9][10][11][12][13]. The absorptive part of the corresponding vacuum polarization is not yet known at this order. The N 3 LO corrections have been computed in ref. [14] and checked in ref. [15] (except for their kinematic π 2 terms) and very recently in ref. [16]; see refs. [5,17,18] for the previous orders.
In this article we present this hitherto missing fourth-order correction, thus completing the N 4 LO corrections to Higgs decay into hadrons in the limit of a heavy top quark and any number of massless flavours. We have also performed the computationally far simpler N 4 LO determination of the H →bb decay rate and verified the result of ref. [4]. A somewhat more demanding but closely related computation is that of the N 4 LO corrections to the electromagnetic R-ratio for the process e + e − → hadrons. So far these corrections were determined only by one group [19][20][21][22]. We have also re-calculated this quantity to the fourth order, and find complete agreement for both the non-singlet and singlet contributions.
Our N 4 LO computations employ the same overall strategy as those by the Karlsruhe-Moscow group mentioned above (see also ref. [23]): the pole terms of the relevant correlation function are calculated in dimensional regularization [24,25] at five loops, and subsequently the absorptive part is extracted. Our use of this approach has been made possible by the development of (a) Forcer [26][27][28], a Form [29][30][31] program for the parametric reduction of four-loop self-energy integrals, and (b) a program [32] efficiently implementing the R * operation, see refs. [33][34][35][36][37], locally for the evaluation of L -loop pole terms in terms of (L −1)-loop integrals. In order to cope with the computations for H → gg, which are far more demanding than those required to determine the five-loop beta function [38,39], the latter program has undergone substantial modifications and extensions.
The remainder of this article is organized as follows: in section 2 we define our notations and briefly address some computational details. Our new N 4 LO result for the H → gg decay width is presented and discussed in section 3. Due to the rather large higher-order coefficients in the expansion in α s , it is interesting to compare the results in the standard MS scheme [40,41] to those in a fairly common (and, in certain contexts, more physical) alternative, the miniMOM scheme [42,43]. The transformation to this scheme, in contrast to other MOM schemes, is known to N 4 LO [44]; it has argued to be preferable to MS for H → gg in a recent N 3 LO study [45].
In section 4 we briefly address the decay H →bb. We present the N 4 LO correction for a general gauge group and a general renormalization scale which has not been written down in the literature so far. The N 4 LO results in QCD, which show a far less problematic behaviour at the only relevant case of α s ≈ 0.1 than their H → gg counterparts, have been known and discussed for more than ten years. Hence there is no need to go into more detail in this case. This is somewhat different for the R-ratio addressed in section 5, despite its even smaller coefficients in the expansion in α s , since this quantity is of physical relevance down to rather low scales and correspondingly high valus of α s . Hence the size and scale (in-) stabilily of this quantity is illustrated in both the MS and the miniMOM scheme at two low-scale reference points. We briefly summarize our results in section 6.

Theoretical framework and calculations
Inclusive Higgs-boson decay to gluons In the limit of a large top-quark mass and n f effectively massless flavours, the decay of the Higgs boson to hadrons can be calculated using the effective Lagrangian [5,18] Here H is the Higgs field, and G µν a the renormalized gluon field-strength tensor for QCD with n f flavours and the Lagrangian L QCD(n f ) . The renormalized coefficient function C 1 includes the top-mass dependence. G F ≃ 1.1664 · 10 −5 GeV −2 denotes the Fermi constant.
At the leading order (LO) of perturbative QCD, eq. (2.1) implies that the Higgs decays to hadrons only via H → gg. At the (next-to-) n -leading order, N n LO, up to n additional partons occur in the final state. As usual, we will refer to the inclusive decay induced by eq. (2.1) as H → gg also beyond LO. The corresponding partial decay width Γ H→ gg can be related, via the optical theorem, to the imaginary part of the Higgs-boson self energy: where M H is the Higgs boson mass, δ is an infinitesimally small positive real parameter and Π GG denotes the contribution to the self energy of the Higgs boson which is induced by its effective gluonic couplings as produced by eq. (2.1).
The Wilson coefficient C 1 can be extracted via a low-energy theorem from a decoupling relation [9] which relates the value of the strong coupling in a theory with n f light flavours, to its value α (n f +1) s in the corresponding theory with n f light flavours and one heavy flavour. The analytic QCD expression for C 1 up to N 4 LO has been provided in ref. [13] as a function of α (n f +1) s at the renormalization scale µ = µ t , where µ t = m t (µ t ) is the scale invariant (SI) top quark mass, i.e., the MS mass evaluated at scale µ t . Using the decoupling relation [10][11][12][13], the renormalization group and the three-loop relation between the MS mass and the onshell (OS) mass [46,47], we have rewritten the four-loop Wilson coefficient as a function of α (n f ) s (µ 2 ) at an arbitrary renormalization scale µ for the SI, MS and OS top quark masses. The same has recently been done to three loops for the OS scheme in ref. [3].
For the convenience of the reader we include the resulting analytic expressions for C 1 . These are presented in the form Here X labels the mass scheme employed, and we have indicated the reduced coupling a s that we employ for all analytic expressions. The first two coefficients are the same in the above top-mass schemes up to the different definitions of masses entering L t = ln(µ 2 /m 2 t ), where ζ n denotes the values of the Riemann ζ-function and a n = Li n ( 1 2 ) = ∞ k=1 (2 k k n ) −1 . The corresponding expressions for the MS and OS masses are given by Our first calculation of the second component of eq. (2.2), Im Π GG , to N 4 LO is addressed below; for a typical Feynman diagram see the left part of Figure 1. The results are presented and combined with C 1 to N 4 LO results for Γ H→ gg in section 3.

Higgs decay to bottom quarks and the R-ratio
As ref. [4], we compute the inclusive Higgs decay to bottom quarks at N 4 LO in the limit of a small bottom mass, keeping only the leading term proportional to the Yukawa coupling. The corresponding partial decay width can be extracted, again via the optical theorem, from the imaginary part of the bottom-Yukawa induced Higgs-boson self energy Π BB , A diagram contributing to this process is shown in the right part of figure 1.  The third observable we consider is the hadronic R-ratio, see refs. [19][20][21][22] and references therein, defined as Away from the Z-pole, the most important contribution to R(s) is given by the partial decay width of an off-shell photon into massless quarks. Here we re-compute the N 4 LO QCD corrections to this electromagnetic contribution. Analogous to the Higgs decay, this quantity can be extracted from the imaginary part of the photon self energy Π µν (q 2 ) = (−g µν q 2 + q µ q ν ) Π(q 2 ) (2.14) via (2.15) with N R = 3 in QCD. The sum runs over n f quark flavours f with electromagnetic charges e f . The functions r(s) and r S (s) represent the respective non-singlet and singlet contributions to the R-ratio. Example diagrams for these two contributions are shown in Figure 2.

Calculations
For all three observables under consideration, we are interested in the imaginary parts of self energies. These can be readily obtained by analytic continuation, is the dimensional regulator and L the number of loops. The crucial point is now that the imaginary part of the self energy is suppressed by a factor of ε: Consequently the finite part of Im Π(−q 2 ) can be obtained from the 1/ε term of Π(q 2 ).
To compute the single poles we employ the R * -operation for Feynman diagrams with arbitrary numerators [32] to express the poles of five-loop diagrams in terms of four-loop diagrams. The R * -operation thus allows us to compute all ingredients required here using the Forcer program [26][27][28], which automates the reduction and calculation of massless four-loop self energy diagrams. The same approach was used in ref. [39] to compute the five-loop beta function for an arbitrary simple compact gauge group.
However, the Higgs decay to gluons poses a much greater computational challenge: the diagrams are all quartically divergent. In order to infrared rearrange the diagrams, the superficial degree of divergence of the diagrams must be logarithmic. We achieve this by computing the fourth order coefficient of the Taylor expansion in the external momentum q about the point q = 0, i.e. we apply the differential operator to all Feynman diagrams. As a result, an 'explosion' of terms, with complicated numerator structures, is created. To deal with this complexity, we have significantly improved our algorithms, in particular for the reduction of high rank tensor vacuum graphs.
The Feynman diagrams for all three cases have been generated using QGRAF [48] and were then processed by a Form [29][30][31] program that assigns the topology and determines the colour factor using the program of ref. [49]. Diagrams of the same topology, colour factor, and maximal power of n ℓ have been combined to meta diagrams for computational efficiency. Lower-order self-energy insertions have been treated as described in ref. [50].
In the case of Π GG , this procedure leads to 1 one-loop, 5 two-loop, 38 three-loop, 394 four-loop and 6405 five-loop meta diagrams. These are fewer meta diagrams than for our calculation of the 5-loop beta function using the background-field method (by a factor 0.68 to 0.69 beyond two loops), but the present diagrams are much harder, as discussed above. The computations were performed on the same set of a modern and somewhat dated machines as used for our five-loop beta function [39], and required an order of magnitude more time.
In the much more modest cases of Π BB and Π in eqs. (2.12) and (2.15), for which we can use the same diagram set which different external vertices and projections, we computed 1 one-loop, 2 two-loop, 9 three-loop, 64 four-loop and 804 five-loop meta diagrams.
We have checked our results by computing all diagrams by at least two different infrared rearrangements. A different rearrangement results in the computation of a different set of counterterms, but should give the same result in the end. This therefore constitutes a highly non-trivial consistency check of our setup.
The first strategy of IR rearrangement consists of attaching external momenta around the line with the worst IR divergence, for example The resulting integral is an L-loop 'carpet' integral, which can be reduced to a L −1 loop propagator integral [28]. By attaching the external momenta around the worst IR divergent line, the number of counterterms that include this line is limited.
The second IR rearrangement consists of inserting a mass into the worst IR divergent propagator: (2.20) The resulting counterterm diagrams can always be split up into a massive one-loop vacuum bubble and an L−1 loop massless propagator integral which can be computed using Forcer.
The advantage of this method is that the massive line cannot be part of any IR counterterm and that the 'carpet' rule reduction is avoided. Overall, this rearrangement is about 20% to 50% faster than attaching external momenta.

Higgs decay to gluons
After the calculation of the Feynman diagrams, the extraction of the absorptive part and its renormalization, the coefficients g n up to N 4 LO in with N A = 8 in QCD, are found to be and Eqs. (3.2) -(3.4) agree with the previous results in refs. [5,[14][15][16][17][18] (n ℓ instead of n f is often used for the number of light flavours) eq. (3.5) represents the main new result of the present article. In all these equations T F = 1/2 has been inserted; this factor can be re-instated by substituting n f → 2 T F n f in all terms that do not involve quartic group invariants.
The coefficients (3.2) -(3.5) are valid for the standard choice µ 2 = q 2 of the renormalization scale. The additional terms for µ 2 = q 2 can be obtained from the scale invariance of (β(a s )/a s ) 2 Im Π GG (q 2 ) [5,14]. This can be done, e.g., by inserting the expansion of a s (q 2 ) in terms of a s (µ 2 ) which can be read off to the order required here, for example, from eq. (2.9) and footnote 2 of ref. [51]. The resulting generalizations of eqs. (3.2) -(3.5) read in terms of the above coefficients g n , the coefficients β n of the MS beta function up to N 3 LO [52,53] and L q ≡ ln(q 2 /µ 2 ). The resulting explicit coefficients up to g 3 agree with eq. (26) of ref. [3], where the definitions of L and a s are slightly different.
At the scale µ 2 = q 2 the numerical expansion of the function G(q 2 ) in eq. (3.1) is given by  for QCD with up to 5 quark families, i.e., n f = 1, . . . , 9 light flavours. In the only physically relevant case of n f = 5 the effect of the fourth-order correction is larger than that of the previous order for α s > ∼ 0.1. It is clear from eqs. (3.7), though, that this is not a generic feature of the QCD perturbation series, but a consequence of the 'accidentally' small size of the third-order term for this number of flavours. A similar situation has been observed for Higgs decay to bottom quarks, see eq. (8) of ref. [4] and eq. (4.6) below.
The fourth-order coefficient g 4 in eq. (3.5) is the first to receive contributions from quartic group invariants. The overall effect of these terms is small in the range of n f considered above; nullifying all these terms changes the coefficients of α 4 s in eqs. (3.7) by about 5% or less for n f = 3. For n f = 3, the relative effect is larger since the coefficient is atypically small.
As discussed in ref. [14], the ζ 2 = π 2 /6 contributions in eqs. (3.3) -(3.4) only arise from the analytic continuation (2.16) and are predictable from lower-order results. The same holds for the terms linear in ζ 2 in eq. (3.5). However, the 'genuine' five-loop contributions from the functions Π GG (q 2 ) include terms with ζ 2 2 , so not all powers of π 2 are 'kinematical' from this order onwards. The numerical decomposition of the expansion in eq. The cancellations between the genuine and kinematic contributions are somewhat less striking than for the corresponding contribution to H →bb, see eq. (7) of ref. [4], yet the conclusion remains the same: it is not possible to obtain reliable results without computing the genuine contributions.
The decay rate Γ H→ gg in the limit of a heavy top quark and n f effectively massless flavours is obtained by combining eqs. The mass-and scale-dependent expansion coefficients γ n for the physical case n f = 5 in where the top-mass logarithms are now given by L tH = ln(M 2 t /M 2 H ). The size of the higher order corrections and the improvement of the renormalization scale dependence from NLO to N 4 LO is illustrated in Fig. 3 for (β(a s )/a s ) 2 Im Π GG (M 2 H ), recall the discussion above eq. (3.6), and for Γ H→ gg in eq. (3.12). The first term in the expansion has been normalized for both quantities, i.e., Γ 0 in the figure is given by  The uncertainty due to the truncation of the perturbation series at N 4 LO is definitely much smaller than the uncertainty due to that of α s (M Z ) which may exceed the value of 1% quoted by the Particle Data Group [55]; see ref. [54] for a recent deviating analysis.
We conclude our discussion of Γ H→ gg by re-expressing its perturbative expansion in another renormalization scheme, the miniMOM scheme [42,43]. The transformation from MS to miniMOM and the beta function in this scheme have been derived at N 4 LO in ref. [44].
The decay width (3.12)  The resulting perturbative expansion of Γ H→ gg in the Landau-gauge miniMOM scheme is shown in fig. 4. The general pattern in miniMOM is somewhat different from that in MS -qualitatively the curves appear shifted to the right. Yet the overall scale range for the interval in µ displayed in the figure is very similar to (if slightly wider than) that in the MS scheme and covered by eq. (3.16). Given this small uncertainty, further investigations of 'optimized scale settings', as performed at N 3 LO in ref. [45] are not warranted.

Higgs decay to bottom quarks
We denote the perturbative expansion of the function R(q 2 ) in eq. (2.12) by 1 N R R(q 2 ) = 1 + n=1r n a n s (q 2 ) (4.1) in terms of the reduced coupling defined in eq. (2.4). The coefficients up to order a 4 s read, for QCD and its generalization to any simple compact gauge group,  Only the case of n f = 5 is phenomenologically relevant, yet it is instructive to briefly consider the numerical dependence of R on the number of light flavours n f over a wide range, The main trend is similar to that of the larger coefficients for Γ H→ gg in eqs. (3.10) and (3.11): the n f -dependent coefficients decrease with increasing n f . The main difference is that the fourth-order term is already negative at n f = 1. The α 3 s -term changes sign close to n f = 7, leading to the large fourth-order / third-order ratio at n f = 5 already observed in ref. [4].
The break-up of the coefficientsr for QCD into 'genuine' and 'kinematic' contributions can be found in eq. (7) of ref. [4]. The numerical scale dependence of R has been included in the comprehensive study of Higgs decays to hadrons to order α 4 s in ref. [3]. However, the scale dependence of the coefficientsr n is available in the literature only to order α 3 s [57]. For the convenience of the reader, we therefore conclude our brief account of Γ H→bb by writing down the generalization of the coefficients (4.2) -(4.5) to a general scale µ 2 : in terms of L q = ln(q 2 /µ 2 ),r n in eqs. (4.2) -(4.3), the coefficients β n of the beta function, and the coefficients γ n of the mass anomalous dimension in the MS scheme up to N 3 LO, see refs. [58,59] and references therein. The coefficients tor 3 (L q ) agree with eq. (17) of ref. [57].

The electromagnetic R-ratio
The non-singlet and singlet contributions to the electromagnetic R-ratio, r(q 2 ) and r S (q 2 ) in eq. (2.15), can be expanded in the same manner as N −1 R R(q 2 ) in eq. (4.1). At µ 2 = q 2 the coefficients for the dominant non-singlet part read Additional singlet contributions enter from the third order in α s , viz with d abc F d abc F /N R = 5/18 in QCD; for the 'time-dependent' normalization of this colour factor see the discussion below eq. (30) of ref. [62]. The above results are in complete agreement with previous calculations, see refs. [19-22, 60, 61] and references therein. The fourth-order coefficients (5.4) and (5.6) had not been verified by a second calculation before.
The numerical expansion of the non-singlet contribution r(q 2 ) in QCD is given by The physically relevant numbers of effectively massless flavours are n f = 3, . . . , 6. The overall effect of the quartic group invariants on the fourth-order coefficient is between 3% and 5% for these values of n f .
The α s -corrections in eq. (5.7) are much smaller than their counterparts for H → gg in eq. (3.7) and H →bb in eq. (4.6); as in the latter case the n f -dependent coefficients decrease with increasing n f . The fourth-order correction is largest for low values of n f : r 4 amounts to 5.6 times r 3 at n f = 1. For three flavours the α 4 s correction contributes as much as the α 3 s terms at α s ≃ 0.3. This situation is at least exacerbated by the kinematic π 2 contributions, as can be seen from the example decompositions where, as above, those contributions have been underlined. For the full n f -dependence of this decomposition see eq. (7) of ref. [19].
The generalization of the coefficients in eqs. (5.1) -(5.4) to µ 2 = q 2 can be obtained from eqs. (4.7) by dropping the terms with γ n which arise from the Yukawa-coupling prefactor in agreement with the corresponding parts of eq. (3.5) -(3.10) in ref. [63] (see also Ref. [64]), where the expansion has been written down in terms of a s = α s /(4π). The qualitative pattern in eq. (5.9) is rather different from that in eq. (5.7): here the ratios of the third-order and second-order coefficients are large. If the fourth-order results were not known, one might by tempted to expect a further rapid growth of the coefficients at this order. Yet, the actual numbers are much smaller than their third-order counterparts for the physical values of n f .
The generalization of eq. (5.9) to µ 2 = q 2 is again given by eqs. (4.7) with γ n = 0, but with the MS beta function replaced by its miniMOM counterpart [42,43]. The µ-dependence of the R-ratio up to order α 4 s is shown for both schemes in figs. 5 and 6 at two low-q 2 reference   points. The first is above thecc resonances but below the Υ threshold, where an analysis with n f = 4 is appropriate [21]. The second is below the J/ψ resonance with n f = 3.
The respective scales are specified via (order-independent) MS values of α s (q 2 ), the corresponding miniMOM values of α s are rounded results of the N 4 LO conversion of ref. [44]. The MS scale variation in fig. 6 amounts to R ns −1 = (6.51±0.11)·10 −2 at N 4 LO. The corresponding miniMOM band is consistent with this result, and only slightly wider if the small-µ spike is not taken into account. For a very low q 2 with α s (q 2 ) = 0.3, the results become unstable below about µ = q in MS and µ = 2q in miniMOM. Disregarding these regions, the N 4 LO results are fairly stable with a 3% uncertainty and R ns −1 = (9.5 ± 0.3) · 10 −2 in the MS scheme.

Summary
We have completed the N 4 LO corrections, i.e., the contributions of order α 6 s , for the decay of the Higgs boson to hadrons via H → gg at the leading order in the heavy-top limit. This correction is slightly smaller than the 1/m top effects at NNLO [7] but, in all likelihood, larger than the presently unknown 1/m top correction at N 3 LO. Hence our result provides an improvement of the overall accuracy of Γ H→ gg . The remaining uncertainty due to the truncation of the perturbation series can be estimated, rather conservatively, as ±0.6%. This is much smaller than the uncertainty of 2.5% induced by a 1% uncertainty of α s (M 2 Z ). An experimental uncertainty of 1% is, of course, not feasible at the LHC. However, a future linear e + e − collider may be able to address the coupling of the Higgs boson to gluons at this level, see section 2.3 of ref. [65].
Furthermore we have calculated, also for a general gauge group, the fourth-order corrections to H →bb in the massless limit and to the electromagnetic R-ratio for e + e − → hadrons. These corrections have been computed by one group before in refs. [4] (where the result is presented only for the gauge group SU(3)) and [19][20][21][22], respectively; we find complete agreement with those results.
Our calculations have been performed using the Forcer program [28] and, at five loops, a Form implementation of the R * -methods introduced very recently in ref. [32]. These methods differ substantially from the global R * -method used in refs. [4,[19][20][21][22]38]. Up to four loops we were easily able to keep all powers of the gauge parameter. This was prohibitively expensive at the five-loop level where we worked in the Feynman gauge. These results have been checked in two ways. The first is verifying the correctness of the higher poles in ε, which have to cancel against the effect of lower-order diagrams in the renormalization procedure. The second check is the computation of all five-loop diagrams in at least two different ways.
We have illustrated the size and renormalization-scale dependence of Γ H→ gg and the R-ratio in both the standard MS scheme and the miniMOM scheme. The N 4 LO results and their stability in these two schemes are comparable for Γ H→ gg and R-ratio at high scales q 2 ; the miniMOM results for the R-ratio at renormalization scale µ 2 ≃ q 2 become unstable at higher q 2 than their MS counterparts. Overall the miniMOM scheme does not appear to be preferable over MS in cases, such as the ones considered here, where the perturbation series is known to a high order.
Form files with our results can be obtained from the preprint server http://arXiv.org by downloading the source of this article. They are also available from the authors upon request.