Higgs boson decay into photons at four loops

Future precision measurements of Higgs boson decays will determine the branching fraction for the decay into two photons with a precision at the one percent level. To fully exploit such measurements, equally precise theoretical predictions need to be available. To this end we compute four-loop QCD corrections in the large top quark mass expansion to the Higgs boson--photon form factor, which enter the two-photon decay width at next-to-next-to-next-to-leading order. Furthermore we obtain corrections to the two-photon decay width stemming from the emission of additional gluons, which contribute for the first time at next-to-next-to-leading order. Finally, we combine our results with other available perturbative corrections and estimate the residual uncertainty due to missing higher-order contributions.


Introduction
The decay of the Higgs boson into two photons was among the main discovery channels at the Large Hadron Collider (LHC) [1,2] and plays a crucial role in precision studies of its properties. Among these are measurements of its mass, see e.g. [3], as well as the study of the interference between the di-photon signal with the continuum background, which leads to a shift in the di-photon mass spectrum and in production rates, thus allowing for studies of the total Higgs boson width [4][5][6][7][8][9]. In addition, future colliders will allow the measurement of ratios of branching fractions at the sub-percent level (see e.g. [10]), thus demanding theoretical predictions of partial decay widths at the same level of precision.
While the next-to-next-to-leading order (NNLO) top quark-induced corrections amount to less than one percent of the full decay width, it is possible that their smallness is accidental and higher order corrections are of the same size. Furthermore, the scale and scheme dependence need to be quantified for a serious estimate of the theoretical uncertainty. This work tries to cover both of these points by performing a four-loop computation of the top quark-induced QCD corrections, combining them with all available higher-order corrections and finally quantifying the theoretical uncertainty on the partial decay width into photons.
The amplitude for this decay can be written as where p i and i are the momenta and polarization vectors of the photons, with p 2 i = 0 and i · p i = 0. The Lorentz-scalar function A(s) only depends on the centre-of-mass energy s = 2p 1 · p 2 , as well as the masses of internal particles. It enters the h → γγ decay width as where M h is the mass of the Higgs boson. The leading-order (LO) (one-loop) contributions to A(s) have been known for a long time [11,12]. They can be decomposed into W-Boson contributions, A W , as well as fermionic contributions, A t,b,τ , due to the top quark, bottom quark and tau lepton.
Next-to-leading-order (NLO) QCD corrections to A t have been computed in the limit of a large top quark mass [13][14][15][16][17] and later with the full mass dependence [18][19][20][21], thus also providing NLO QCD corrections to A b . NNLO QCD corrections in an expansion for a large top quark mass have been obtained in [22,23] and recently, exact numerical results became available [24]. Partial results in the limit of an infinitely heavy top quark are known at next-to-next-to-next-to-leading order (N 3 LO) [25]. NLO electroweak corrections, for which the above decomposition is not possible, have been obtained in [26][27][28][29].
Virtual QCD corrections to A t,b do not lead to infrared singularities, so one does not need to include real-radiation contributions to obtain a finite result. Nonetheless, these realradiation contributions should be considered since in principle the emission of additional, possibly soft, gluons leads to a shift in the di-photon mass spectrum. At NLO, such contributions vanish due to colour conservation; the first non-vanishing contribution arises from the process h → γγgg at NNLO. To our knowledge, the size of these contributions has not been investigated in the literature.
In this work, we concentrate on A t and compute the N 3 LO virtual contributions in an asymptotic expansion for a large top quark mass. We additionally consider the NNLO realradiation contributions in the same approximation. The rest of this article is structured as follows: in Section 2 we present the method by which we compute the N 3 LO virtual corrections, followed by a discussion of the NNLO real-radiation contributions in Section 3. Finally we discuss the numerical impact of both of these contributions in Section 4.

Calculation of N 3 LO virtual corrections
Our computation of the virtual corrections follows [30], in which we computed the top quark contributions to the gluon-gluon-Higgs form factor. In the following, we discuss the computational setup and present analytical results for the h → γγ form factor. These corrections can not be obtained from the results of [30], since starting from three loops it is not possible to simply substitute colour factors.

Computational setup
We generate 2 one-loop, 12 two-loop, 206 three-loop and 5062 four-loop diagrams contributing to A t using QGRAF [31]. We then generate FORM [32] code for each diagram using q2e [33,34], compute the colour factors using COLOR [35] and perform an expansion-bysubgraph [36] in the limit M 2 t s = M 2 h using exp [33,34]. As a result, all diagrams are mapped onto one-to four-loop massive tadpole integrals and one-to three-loop massless form factor integrals, separating the two scales M t and s.
Both types of integrals have received a lot of attention in the literature; all required master integrals are known analytically [37][38][39][40][41][42][43][44][45]. Due to the asymptotic expansion, we need to deal with tensor integrals, i.e. integrals for which we can not re-write all numerator structures in terms of inverse propagators of the integral family. As a consequence, we have to perform a tensor reduction of the numerator structures. Since in the fully hard region of the asymptotic expansion only tadpole integrals can appear, and in all other regions both types of integrals appear, it is advantageous to perform the tensor reduction for the tadpole integrals. For one-and two-loop integral families, there are general algorithms for treating tensor tadpole integrals of an arbitrary rank [46]. These are implemented in MATAD [47]. At three and four loops, we have implemented routines to reduce tensors up to rank 10. This is a sufficient rank to expand the form factor up to M −6 t . After tensor reduction we perform an integration-by-parts reduction of the three and four loop tadpole integrals, as well as the massless form factor integrals, in order to obtain the h → γγ form factor in terms of master integrals. For this we use LiteRed [48,49] and FIRE6 [50].
We then proceed to renormalize the bare strong coupling constant α 0 s and the top quark mass m 0 t in the MS scheme. For this we require the strong coupling renormalization constant at two loops, (since the leading-order contribution to the form factor does not depend on α s ) and the quark mass renormalization constant at three loops 1 . This renders the form factor finite, as there are no infrared divergences present. This is in contrast to the gluon-gluon-Higgs form factor discussed in Ref. [30]. Since the top quark does not appear as a dynamical degree of freedom at energy scales of the order of the Higgs boson mass, we transform α s from the six-to the five-flavour scheme, using the two-loop decoupling constant 2 . Furthermore, we also produce the form factor with the top quark mass renormalized in the on-shell scheme, which we obtain from the MS result by applying the MS-OS relation at three loops [53][54][55][56]. from four loops, there are diagrams such as the last of the second line, in which the photons couple to different light-quark loops; these also sum to zero for the same reason. For the three contributions of Eq. (3) we obtain, for A t,f = α A Here, ζ n is the Riemann ζ-function evaluated at n, a n = Li n (1/2), To reduce the size of the expressions the colour factors have been set to the following values: however the full expressions, including the colour factors and for an arbitrary scale µ 2 , are provided in the ancillary files of this paper.
This result passes a number of cross-checks. Firstly our result is finite after renormalization and, in the abelian limit, agrees with [30]. Secondly, the part stemming from the fully hard region in the asymptotic expansion agrees with [25].

Calculation of NNLO real corrections
Starting from NNLO in QCD, corrections with two gluons in the final state can also be taken into account 3 . While these corrections do not lead to infrared divergences and do not lead to large logarithms, they can potentially impact the di-photon invariant mass spectrum, thus possibly influencing extractions of the total width and Higgs boson mass.
In the following we compute the contributions to the inclusive Higgs boson decay width.
We employ the method of reverse unitarity [58] to express the phase-space integral over the squared decay amplitude for h → γγgg as cut two-point forward-scattering diagrams.
In each of these diagrams the Higgs bosons couple to different top-quark loops, which are connected to each other by the cut gluons and photons. A sample diagram is shown in Fig. 2. A direct large-mass expansion of such five-loop diagrams is computationally expensive, thus we adopt the building-block approach described in [59,60]; we pre-compute the 24 one-loop diagrams contributing to the hγγgg vertex in the large-mass expansion.
Since there are only hard subgraphs, the result is polynomial in all external momenta and can be used as an effective vertex. As a consequence, instead of 120 five-loop diagrams, we only expand 24 one-loop diagrams which then are inserted into the cut three-loop twopoint diagram shown in Fig. 2. The resulting diagram can be computed using MINCER [61]. While each of the one-loop five-point diagrams starts at M 0 t in the large-mass expansion, their sum starts at M −4 t . Consequently, the expansion of each five-loop diagram also starts at M 0 t but their sum starts at M −8 t . This can be understood by the fact that gauge-invariant operators leading to a γγgg interaction need to contain two gluonic and two photonic field-strength tensors and thus are dimension-8 operators.
Computing the individual cut five-loop diagrams would require the computation of five terms in the large-mass expansion, the first four of which ultimately cancel in the sum. Using instead the building-block approach allows us to compute contributions to Γ h→γγgg to order M −12 We obtain We note that by using FORCER [62] to compute four-loop two-point phase-space integrals, it would be a relatively straightforward task to compute also the N 3 LO real and realvirtual contributions. As demonstrated later in Eq. (17) the numerical impact of the NNLO contributions is very small, so we do not perform this computation here.

Numerical results
In the following we study the numerical impact of the N 3 LO top quark contributions, including renormalization scale and scheme dependence, before combining them with the NNLO QCD corrections for bottom and charm quark loops and the NLO electroweak corrections. Finally, we discuss sources of uncertainty on the partial decay width.

N 3 LO top quark contributions
In the following we discuss the numerical impact of the corrections obtained in the previous sections on the Higgs boson decay width into photons. To this end, we employ the exact LO result for A(s), neglecting light fermion loops, and combine it with the large-mass expansion results for A t starting from NLO.
The parameters entering the evaluation are given by [63]: We use the program RunDec [64,65] to evolve the strong coupling constant to different renormalization scales and to obtain the top quark mass in the MS scheme.
For the renormalization scale µ = M h we obtain the decay width with the top quark mass renormalized in both, the on-shell and the MS scheme Γ OS h→γγ × 10 6 GeV −1 = 9.1322 + 0.1558 + 0.0029 − 0.0031 = 9.2878 , Γ MS h→γγ × 10 6 GeV −1 = 9.1188 + 0.1639 + 0.0070 − 0.0023 = 9.2874 , where we have separated the LO, NLO, NNLO and N 3 LO contributions. While the NNLO corrections are two orders of magnitude smaller than the NLO corrections, the N 3 LO corrections are larger than the NNLO corrections in the on-shell scheme, but negative.  In the MS scheme they are also negative and the same order of magnitude, partially cancelling the effect of the NNLO corrections. In both cases, the effect on the decay width is small, resulting in a 0.034% and a 0.024% decrease, respectively. This apparently bad perturbative convergence raises the question of whether this behaviour is an artifact of the renormalization scale choice. In Fig. 3 we show the on-shell results and in Fig. 4 the MS results, normalized to the LO at µ = M h . In both cases, the renormalization scale dependence of the N 3 LO correction is slightly more stable than at NNLO. In the on-shell scheme the renormalization scale dependence decreases from 0.013% to 0.009%, whereas in the MS scheme the scale dependence decreases from 0.013% to 0.006%. The near-exact cancellation between the NNLO and N 3 LO corrections in the on-shell scheme only occurs for µ ≈ M h . However, in the whole range of renormalization scales considered, the N 3 LO correction is of the same order of magnitude as the NNLO correction, and negative.
It is instructive to study the subsequent terms of the large-mass expansion entering the N 3 LO correction to the decay width. In the case of the logarithm of the gluon-gluon-Higgs form factor it was found at N 3 LO, in the on-shell scheme, that the first mass-suppressed term amounts to roughly 10% of the leading term, and the second mass-suppressed term amounts to 1% [30]. This is in contrast to lower orders, where the convergence is better. However, the gluon-gluon-Higgs form factor is infrared divergent and thus not a physical observable, unlike the decay width which we consider here. Normalizing the NNLO and N 3 LO contributions to Γ OS h→γγ to their leading term in the large-mass expansion, we obtain  for µ = M h : where we have separated the ρ 0 , ρ 1 , ρ 2 and ρ 3 contributions. Thus, both for NNLO and N 3 LO the mass-suppressed terms are larger than for the gluon-gluon-Higgs form factor. However, the mass corrections at N 3 LO are only slightly larger than their NNLO counterparts. Taking µ = √ 2M h as a renormalization scale, we obtain In all of the above cases, the mass-suppressed terms amount to at least 25% of the leading term, and in some cases as much as 400%. However, in all cases, the corrections at order ρ 3 are sufficiently small. Also in all of the above cases, the MS scheme displays better convergence. For brevity we do not show them here.
In the case of the gluon-gluon-Higgs form factor, it was observed that the mass-suppressed terms become increasingly important at higher perturbative orders. However in this case, no such clear statement can be made; the convergence of the large-mass expansion at different perturbative orders depends strongly on the chosen value of the renormalization scale. The scheme difference decreases from 0.013% to 0.004% for µ = M h and the decrease is approximately constant across renormalization scales in the range A final remark concerning the NNLO real-radiation corrections, computed in Section 3, is in order. They amount to Γ OS h→γγgg × 10 16 GeV −1 = 1.68 + 0.48 + 0.09 = 2.25 , where we show the individual terms in the large-mass expansion. They converge well, but are clearly negligible 4 .

Including light quarks and electroweak corrections
We are now in the position to combine the N 3 LO corrections computed in the previous section with the known NLO electroweak corrections [29], as well as the NNLO QCD corrections with massive bottom and charm quark loops [24]. To this end, we convert the results of [24], which are given in terms of the on-shell quark mass, into the MS scheme and transform α s to the five-flavour scheme. A detailed discussion of the conversion can be found in Appendix A. Note that here the top quark mass is in the on-shell scheme, to be compatible with the NLO electroweak corrections of [29]. Furthermore, we include τ loops at LO.
For the light fermion masses we take [63] m We employ the program RunDec to obtain the two light quark masses in the five-flavour scheme at µ = M h .
As result, we obtain Γ h→γγ × 10 6 GeV −1 = 9.2581| LO − 0.1502| NLO,EW + 0.1569| NLO,t + 0.0157| NLO,bc + 0.0030| NNLO,t + 0.0037| NNLO,bc − 0.0031| N 3 LO,t = 9.2840 , (19) where we split the higher-order corrections into electroweak (EW), top-induced (t) and bottom/charm-induced (bc) contributions 5 . We note that the NNLO corrections coming from bottom and charm quark loops are of the same order as the NNLO top quark contributions, as are the N 3 LO top quark corrections. This seems to indicate that the NNLO top quark contributions are "accidentally small", since the bottom and charm contributions themselves seem to converge well.

Theoretical uncertainties
We can now comment on the sources of theoretical uncertainties on the partial decay width; missing higher-order QCD and electroweak corrections, as well as parametric uncertainties.
The N 3 LO QCD corrections involving top quark loops and the NNLO QCD corrections involving bottom and charm quark loops are both smaller than 0.05% while the scale and scheme dependence is even smaller. Furthermore, the mixed quark-flavour contributions should be well approximated by our treatment (see Appendix A). Thus, we assign an overall uncertainty of 0.1% for missing higher-order QCD corrections.
In contrast to the pure QCD corrections, the effects of three-loop electroweak and mixed QCD-electroweak corrections are harder to estimate. While naive arguments, such as the O (α s ) suppression w.r.t. the NLO corrections would suggest that they are smaller than 1%, counter-examples for similar quantities exist in the literature. For example in the case of the Higgs decay into gluons, the NNLO O (G F M 2 t α 3 s ) corrections are of the same order of magnitude as the NLO O (G F M 2 t α 2 s ) corrections [66]. Furthermore, there is a significant cancellation between fermionic contributions at NLO and purely bosonic contributions 6 , which are not necessarily the case at NNLO. As a consequence, we conservatively assign an uncertainty of 1.6% due to missing higher-order electroweak corrections, which is the size of the NLO electroweak corrections.
To obtain the parametric uncertainties, we vary the input parameters within the experimental uncertainties. The dominant parametric uncertainty stems from the Higgs boson mass and amounts to ±0.5%. All other parametric uncertainties are negligible.

Conclusion and outlook
In this paper we investigate the top quark contributions to the Higgs boson decay into photons at N 3 LO in QCD and compute real-radiation corrections at NNLO. While the corrections are small and will not play a role in the precision programme of the LHC, we find that the N 3 LO corrections are of the same order of magnitude as the NNLO corrections and that the mass-suppressed terms are sizeable relative to the leading term in the large-mass expansion; they must be included for a reasonable description of the top quark mass effects.
We combine our results with the recently-obtained NNLO QCD corrections due to bottom and charm loops, as well as NLO electroweak corrections and, for the first time, quantify the residual uncertainty due to higher-order QCD corrections. We obtain where the three different sources of uncertainties are discussed in Sec. 4.3. We find that the largest source of theoretical uncertainty on Γ h→γγ are the missing NNLO electroweak and mixed QCD-electroweak corrections, which would be required to reduce the total uncertainty to 1% and below.

A Light quark contributions
In this section we briefly describe the steps necessary to adapt the results of [24] to the conventions used in this paper. The results of [24] are given for µ = M h , the quark mass renormalized in the on-shell scheme and α s in the four-flavour scheme for bottom quark and three-flavour scheme for charm quark. Thus several steps have to be taken to obtain results for arbitrary µ and the quark masses renormalized in the MS scheme.
We write the bottom quark contribution to the amplitude in terms of the on-shell quark mass M b as We first express α (4) s (M h ) in the five-flavour scheme and for an arbitrary renormalization scale µ, α (4) Here β 0 = 11 3 C A − 4 3 T f n l is the one-loop QCD β-function and ζ (1) is the one-loop contribution to the decoupling constant of α s .
Next, we convert the on-shell mass M b to the MS mass m (5) b (µ). We proceed by replacing the on-shell mass in the exact LO and NLO expressions in terms of the MS mass, according to and then expand in α s . Here c (l) are the l-loop coefficients of the OS-MS relation and depend on µ and m For the NLO contribution only the first line is necessary. Note that here the charm quark is treated as a massless quark.
For the charm quark we start out in the three-flavour scheme, where in addition to the above steps a decoupling of the bottom quark in both α s and m c has to be performed.