$H\to \gamma\gamma$ to all orders in $\alpha_s$ in the large-$\beta_0$ limit of QCD

In this note, we present the result for the QCD corrections to the decay of the Higgs boson into two photons in the large-$\beta_0$ limit of QCD, providing the first two terms in the heavy-top expansion. From our results, one can easily read off the exact leading-$n_f$ QCD contributions in analytic form to all orders in the strong coupling, $\alpha_s$, where $n_f$ is the number of massless quarks, and identify the leading renormalon singularities. We give explicit results for the leading-$n_f$ coefficients at 6 and 7 loops and use the large-$\beta_0$ result to speculate about the size of yet unknown (but small) higher-order contributions to the QCD series.


Introduction
The Higgs decay into two photons played a central role, ten years ago, in the discovery of the Higgs boson at the Large Hadron Collider (LHC).This channel remains crucial for precision studies of Higgs properties at the LHC as well as at future colliders.Given the expected increase in precision on the experimental side, the theoretical calculation of Higgs decays within the Standard Model (SM) must be performed at higher orders in perturbation theory.
Because the photon is massless, this decay is a loop-induced process at leading order (LO) within the SM.Its amplitude can be decomposed into a bosonic contribution, stemming from the W boson, and fermionic contributions.Here we are concerned with the fermionic contribution arising from the top quark, which dominates over the subleading b-quark and τ -lepton fermion loops (contributions from even lighter fermions are irrelevant).
The decay amplitude has been known at LO for many years [1][2][3][4].The next-toleading order (NLO) result was first obtained as an expansion in τ t = m 2 H 4m 2 t ≈ 0.13, with m H and m t being the Higgs and the top-quark masses, respectively, but later the full m t dependence of the NLO result was obtained [5][6][7][8][9][10][11][12].The NNLO QCD contribution, which is a three-loop calculation, was obtained about ten years ago, again as an asymptotic expansion in the top mass; numerical results with full top-quark mass dependence were obtained only much more recently [13][14][15].The NNNLO (four-loop), O(α 3 s ), corrections were first partially computed in Ref. [16] and later completed in Ref. [17].Partial results for the N 4 LO contributions where both photons couple to a massive top-quark loop can also be found in Ref. [16].
In this work we present, for the first time, the result for the QCD corrections in the large-β 0 limit of QCD within the heavy-top limit.From our results, the exact leading-n f QCD contributions, i.e. coefficients of terms proportional to α n s n n−1 f , with n f being the number of active (massless) quark flavors, can be directly read off in analytic form to all orders in the strong coupling.In our calculation, we are concerned with non-singlet diagrams only, i.e. diagrams in which the Higgs boson and the two photons in the final state couple to the same top-quark loop.Contributions from the singlet diagrams are subleading in the large-β 0 power counting, and do not contribute to the terms O(α n s n n−1 f ). 1  In the large-β 0 limit of QCD [18], one considers first the limit of large n f , keeping α s n f ∼ O (1).The fermionic corrections to the gluon propagator acquire a special character since they count as O(1), and a dressed gluon propagator, known to all orders in α s , can be obtained.When no external gluons are present, such as in the case of this work, the leading-n f contributions to a given process can be obtained, to all orders, by replacing the gluon propagator in the calculation of the O(α s ) QCD correction by the Borel transform of the dressed gluon propagator.The large-β 0 limit consists then in the replacement of the fermionic contribution to the QCD β function, β 0f , proportional to n f , by the complete one-loop QCD beta-function coefficient, β 0 .This procedure, known as naive non-abelianization [19,20], effectively resums a set of non-Abelian graphs related to the running of the coupling, thereby restoring the non-Abelian character of the results.
Our calculation is performed with a modified version of the publicly available package MATAD [21], which is written in FORM [22].Working with the Borel-transformed dressed gluon propagator in the large-β 0 limit amounts, essentially, to a perturbative calculation with an analytically regularized gluon propagator.This requires that the exponent of the squared momentum in the gluon propagator be generalized to a real number [18].This was implemented through modifications of the original source code of MATAD, since this generalization is not supported by the standard version.
The importance of calculations in the large-β 0 limit should not be overstated.Albeit realistic, it is a toy model for higher orders, but one that can lead to insights about the physics of renormalons, non-perturbative corrections, and even, in some cases, allow for estimates of unknown higher-order contributions.The fact that the leading-n f terms are exactly obtained can also be useful as an independent cross-check of future calculations.
We present, in analytical form, the Borel transform of the decay amplitude in the large-β 0 limit.The amplitude is computed as an expansion in the heavy-top limit and we provide the first two terms of the corresponding series in τ t .From this result, it is simple to extract in analytical form the exact leading-n f contributions to the QCD series to all orders in α s , with n f being the number of massless quarks.We reproduce the known leading-n f results up to 5 loops, and give for the first time the explicit expression at 6 and 7 loops.The large-β 0 result allows for a discussion of the different renormalon singularities, and we show that the leading ultraviolet (UV) renormalon is very prominent, leading to a sign alternating perturbative series (for not-so-large renormalization scale choices).The large-β 0 series reproduces the approximate magnitude of the known QCD perturbative coefficients, although with incorrect signs.This allows for a speculation about the size of yet unknown contributions to the QCD series.
This work is organized as follows.In Sec. 2 we set up the theoretical framework and the notation.Our main results are given in Sec. 3 together with the discussion of some of their potential implications.In Sec. 4 we present our conclusions.We relegate to App.A a brief description of the modifications to MATAD that are required to work in the large-β 0 limit, while in App.B we present the full results for the leading-n f coefficients, up to O(τ t ), at 6 and 7 loops.

Theoretical framework
The decay width of the process H → γγ starts at one-loop order in the SM and can be written in terms of bosonic and fermionic contributions as where the first amplitude, A W , is due to W boson diagrams, while A f stems from the decay mediated by charged fermions with mass m f [23].In the above expression, , where m W is the W -boson mass.The fermionic contributions are strongly dominated by the top-quark loop, with small subleading components due to the bottom quark and the tau lepton.Here, we focus exclusively on the top-quark contribution, A t .We work in the heavy-top limit and compute A t as an expansion around τ t → 0, keeping the first two terms of the corresponding series.We expect the heavy-top limit (i.e. the leading term) to be sufficient to estimate effects due to higher orders in perturbation theory, since the first τ t correction is of about 10%.We compute the first subleading term to crosscheck this assumption.We write the amplitude A t as an expansion in powers of the strong coupling, α s , as where with α being the fine-structure constant, Q t is the electromagnetic top charge, N c = 3 is the number of colors, and v = 1/ √ 2G F is the vacuum expectation value of the Higgs field with the Fermi constant as input.With this convention (which is the same as Ref. [16], for example), the leading order contribution to the amplitude A t , Eq. ( 2 Using expansion by regions (see for example Ref. [24]), this series can be obtained from the computation of the loop integrals in the hard region, which leads to so-called massive tadpole integrals.This remains true at higher orders for non-singlet diagrams.
The NLO term due to the top loop is obtained through the exchange of a virtual gluon between the internal top propagators [23].The topology of the main diagrams for the NLO process is shown in Fig. 1.There are a total of 12 diagrams at the two-loop level, and the final result for the amplitude up to corrections of O(τ t ) is and m t ≡ m t (µ) is the MS top-quark mass.Here, we have reproduced the NLO calculation using qgraf [25] and MATAD [21], since this forms the basis for the computation in the large-β 0 limit, as we discuss further below.Analytical expressions for higher order corrections, A (2) t and A (3) t , up to four loops, can be found in Refs.[16,17].
We turn now to the setup of the large-β 0 limit calculation.In the large-β 0 limit we deal with the perturbative series to all orders in α s .As is well known, in QCD, series of this type have coefficients that diverge factorially and the perturbative series is, at best, asymptotic.In this context it is convenient to work with the Borel transform of the series, which suppresses the factorial divergence of the coefficients and can have a finite radius of convergence.For the perturbative expansion of a quantity R starting at O(α s ) (without any loss of generality), we define its Borel transform in the following way [18] B The Borel transform is the inverse Laplace transform of the original series.The procedure can be inverted and a 'true value' of the asymptotic series can be assigned (in the Borel sense) through the Laplace transform as provided the integral exists.In the Borel t-plane, singularities known as renormalons appear.They arise from infra-red (IR) and ultra-violet (UV) regions in the loop subgraphs.The IR renormalons are particularly important here, since they appear on the right-hand side and obstruct the integration in Eq. ( 8).Circumventing these singularities generates an imaginary ambiguity in the Borel integral.This ambiguity is related to nonperturbative terms from condensate matrix elements in the operator product expansion (OPE).Here, since the typical scale of the problem is quite high compared with Λ QCD , these OPE corrections -and accordingly the imaginary ambiguities of IR renormalons -are tiny and can, for all practical purposes, be neglected.Using the definition in Eq. ( 7), the dressed gluon propagator with a four-momentum k obtained after resumming the massless-quark bubble-loop corrections in Borel space reads [18] where we introduced the variable u, defined as (We recall that β 0f are the fermionic contributions to β 0 .)In Eq. ( 9), µ2 is the renormalization scale and ξ is the gauge parameter.The constant C sets the renormalization scheme: in the MS scheme, C = −5/3; whereas in the MS scheme, for instance, C = −5/3 + γ E − ln 4π.Notice that, in the definition of the transformation, Eq. ( 9), we multiplied the gluon propagator by α s ; thus, the lowest order term in the u expansion corresponds already to the NLO QCD correction.In Eq. ( 9), the factor (−µ 2 e −C ) u ensures scheme and scale invariance of the Borel integral result, Eq. ( 8).We use the following definition for the QCD β function and its coefficients: with where , and n f is the number of massless quark flavors.The second contribution on the right-hand side of Eq. ( 12) is the fermionic contribution, Working in Landau gauge, with ξ = 0, the Borel transform of the dressed propagator is essentially the original propagator with a modification in the exponent of the denominator momentum, k 2 → (k 2 ) 1+u -which amounts to working with an analytically regularized gluon propagator.In our case, using Eq. ( 9) as the gluon propagator in the NLO QCD calculation is sufficient to generate the exact leading-n f terms to all orders in α s .This is depicted in Fig. 1, where the main topologies involved in the diagrammatic calculation are shown and the dashed line represents the resummed, dressed, gluon propagator.The calculation of the diagrams of Fig. 1 yields the main result of this paper.We implemented this calculation in MATAD with suitable modifications to the source code to account, essentially, for the new gluon propagator.Some technical details of this implementation can be found in Appendix A. The source is also publicly available on GitHub. 2 In the next section we discuss the final results.

Results
The Borel transform of the amplitude A t in the large- where Ât is defined in Eq. ( 3) and P (u) = 126 + 155u + 180u 2 + 100u 3 + 9u 4 .The last line in Eq. ( 13) stems from the relation between the MS and pole mass in the large-β 0 limit [26,27].Removing this line effectively converts the mass scheme from the former to the latter.The function G0 is defined as a power series in u with coefficients g n /n!, where the g n are obtained from the generating function The expression in Eq. ( 13) exhibits the renormalon singularities arising from UV and IR regions of loop sub-graphs, which are encoded in the Γ functions.In the leading term, corresponding to an infinitely heavy top quark, the UV renormalons, which are located at negative integer values of the variable u, are all double poles, with the sole exception of the leading UV pole at u = −1, which is simple.This happens due to a partial cancellation with the (u 2 − 1)/Γ(1 + 2u) term.Note also that there is no singularity at u = −1/2 and the function is regular at this point, as expected.The IR renormalons, on the other hand, are all simple poles with the leading singularity at u = 2 being a quartic IR sensitivity, connected with dimension four corrections in the Operator Product Expansion.Regarding the dominant renormalon singularities another observation is that the residue of the leading UV pole and that of the leading IR pole are of the same order.Their ratio being of order 1 implies that the UV renormalon, being the closest to the origin, should take over rather earlier and lead to a sign alternating series.This, however, depends on the renormalization scale that is chosen: larger values of µ enhance the IR poles and suppress the UV contributions, which postpones the dominance of the UV renormalon, reducing its residue, and consequently delays the sign alternation.The structure of the term proportional to τ t is similar.The main difference here is the appearence of a term without the prefactor µ 2u e 5u/3 , which ensures scheme and scale invariance of the Borel integral.This happens because the last term in Eq. ( 13) is a consequence of the renormalization scheme (and scale) dependence of the quark mass.Finally, we note that the expression is regular at u = 0.
From the Taylor expansion of Eq. ( 13) one recovers the perturbative expansion in α s in the large-β 0 limit.Of course, the leading term must correspond to the exact NLO QCD result, and it does lead to A (1) t of Eq. ( 5), as expected.Then, further expanding in u and performing the substitution u = −β 0,f t, one can re-obtain the leading-n f contributions, i.e. the coefficients of terms proportional to α n s n n−1 f , that we denote of the exactly known QCD corrections in perturbation theory up to O(τ t ).At 3 and 4 loops we get where ζ i is the Riemann Zeta Function evaluated at i and we recall that µ ≡ ln µ 2 m 2 t .These results are in agreement with Refs.[14,16].At 5 loops, we reproduce the result for the leading-n f term in the heavy-top limit [16] and obtain for the first time the leading-n f terms of the O(τ t ) correction as A new result of this work is the leading n f -terms to all orders in perturbation theory.
As an example, in the heavy-top limit, the N 5 LO leading-n f coefficient, A The corresponding O(τ t ) corrections are lengthy and are given in App.B. The higherorder coefficients can be easily generated from Eq. ( 13) following the procedure outlined above.
Let us now discuss the series in the large-β 0 limit, after naive non-abelianization, with the replacement β 0f → β 0 .This transformation preserves the leading n f coefficients A (n) t,n n−1 f and generates approximate results for all sub-leading powers of n f .The perturbative series thus obtained is, at best, a good qualitative approximation to the full QCD results.With this caveat in mind, we show in Fig. 2a the perturbative series for four different choices of the renormalization scale µ, order by order in perturbation theory.The horizontal line is the Borel sum of the series, given by the integral of Eq. ( 8). 3 In Fig. 2b we display the scale dependence of the result, order by order, varying the renormalization scale µ in the logarithms µ (see Eqs. ( 14)-( 18)) and in ) is performed at one loop, for consistency with the large-β 0 limit, using as input α (n f =5) s (m Z ) = 0.1179 while for the top mass the value of m t (m t ) = 163.643GeV [28] is kept fixed, since the MS quark-mass running is of order 1/β 0 itself, which would then generate 1/β 2 0 terms not included in the perturbative coefficients [29,30].The stabilization with respect to scale variations is clearly observed when the order is increased, starting from N 4 LO, or O(α 4 s ), with the N 5 LO result being extremely stable and in excellent agreement with the Borel sum.
Before comparing with the QCD calculation, an observation regarding the singlet diagrams is in order.In Ref. [16], it is shown that with the scale choice µ = m t (m t ) (which resums the logarithms µ ) at the 3-loop level the singlet diagrams are approximately of the same magnitude as the non-singlet diagrams, while at a lower scale, µ = m H , the singlet contribution is significantly enhanced and is approximately three times larger than the non-singlet contribution.It is argued that one might expect the same to happen at 4-and 5-loop orders.We have checked that this is generally the case at the 4-loop level for the singlet diagrams where both photons couple to massless quarks, using the result from Ref. [17].However, here the singlet contribution is suppressed at µ = m t (m t ), before dominating around µ ≈ 270 GeV (due to the non-singlet contribution becoming zero) and becoming relatively smaller again at even higher scales.Thus, if the large-β 0 limit represents well quantitatively the non-singlet diagrams in full QCD, at higher values of µ, the difference between the result in QCD and in large-β 0 should decrease.Below, the results we discuss are for µ = m t (m t ), which reduces the singlet contributions.
After naive nonabelianization, our result for A t,large-β 0 at the scale µ = m t (m t ) reads (a s ≡ α s (m t (m t ))/π) This result shows that in the large-β 0 limit the pattern of sign alternation is already present, systematically, starting at O(α 2 s ).This is a consequence of the dominant behavior of the leading UV singularity.(This pattern disappears if even larger values of µ are used.For example, for µ = 280 GeV the coefficients have fixed signs, as can be seen in Fig. 2a for the leading term in τ t .)As expected, we find that the subleading terms in τ t are generally about an order of magnitude smaller than the leading ones (the O(α 2 s ) coefficient is an exception).We can now compare our result for the leading term in τ t with results from Refs.[16,17].In the same conventions, and again for µ = m t (m t ) and a s = α 25 GeV, we find for the QCD result At O(a 4 s ) only the contributions where the photons couple to a massive top-quark loop are known from Ref. [16] (we have re-expressed the result of Ref. [16] in terms of a s = α (n f =5) s /π) and the constant c 4 has yet to be calculated.But these contributions do dominate the real part of the O(α 3 s ) coefficient, as can be seen comparing the results of Refs.[16,17].The singlet contributions where the photons couple to massless quarks at O(a 2 s ) and O(a 3 s ) are indicated by 'si'.The imaginary part in the singlet contributions is required by unitarity.In the decay width, however, the first contribution arising from an imaginary part, namely from Im(A s ) and can be expected to be suppressed, since it competes with the real parts from higher-order coefficients, which tend to be significantly larger, as seen in the QCD result.We therefore compare the real parts of the coefficients of Eqs.(19) and (20).From these results we see that the magnitude of the coefficients obtained in the large-β 0 limit are not too far from the magnitude of the known QCD coefficients, but the sign is wrong in all cases.In QCD, the systematic sign alternation is not yet established at order O(α 2 s ) although the signs do alternate starting from O(α 3 s ) for the known terms.This is an indication that the UV renormalon in QCD should be less prominent than in the large-β 0 limit, in agreement with other processes where the large-β 0 result is known [30][31][32].
Based on the observation that the large-β 0 result roughly reproduces the size (but not the sign) of the QCD contributions, we can speculate that the α 5 s non-singlet coefficient would have the same magnitude as the one shown in Eq. (19).Given the important contribution of the singlet diagrams, it is safer to attach a 100% uncertainty to this result, which leads to the following estimate for the O(a 5 s ) contribution to A t for µ = m t (m t ): A t ≈ 300 ± 300.
Using this value and the QCD series of Eq. ( 20), in the decay width, this implies a correction of, Γ(H → γγ) which is about as large as the contribution from the partial α 4 s correction calculated in Ref. [16], which amounted to mere 0.02‰ of the total decay width.(For this estimate we have added the LO W contribution and top-mass corrections in the 1-and 2-loop QCD results -see, e.g., the expressions given in Ref. [14].)If significantly lower renormalization scales are used one finds, for µ = m H /2 for example, Γ(H → γγ) which is, at this scale, 2.5 times smaller than the (partial) α 4 s correction.In all cases, we find that the O(α 5 s ) correction is at most as large as the O(α 4 s ) contribution, which means that estimating the truncation error from the last included term is safe.This clearly indicates that the QCD contributions are under very good control here, with an uncertainty from the truncation which is an order of magnitude smaller than the uncertainty arising parametrically from the value of α s itself, for example, and much smaller than the uncertainty from m H which amounts to about 0.04 keV.

Conclusions
In this work we presented the result for the QCD corrections to H → γγ in the large-β 0 limit, providing the first two terms in the heavy-top expansion.The analytical result was obtained with a modified version of MATAD [21].From the Borel transformed amplitude, upon re-expanding in the Borel variable, we can reconstruct the perturbative series to all orders in α s in the large-β 0 limit.In particular, the exact analytical results for the leading-n f contribution at each order can be extracted and, apart from reproducing the known results, we have given explicit expressions for the previously unknown contributions at 6 and 7 loops, in Eqs.(17) and (18).The higher-order leading-n f coefficients can easily be extracted from the analytical result for the Borel transformed amplitude, Eq. ( 13).These results can serve as a partial cross-check for future calculations, should the full QCD corrections be independently calculated by other groups.
Furthermore, the exact knowledge of the Borel transform in the large-β 0 limit allows for a study of the renormalon singularities.We find the usual and expected towers of UV and IR renormalons: the leading UV renormalon, located at u = −1 in the Borel plane, being the closest to the origin and dominates the series at high orders, while the leading IR renormalon is related to a quartic IR sensitivity, and appears at u = 2.In the large-β 0 limit, the UV renormalon has a somewhat large residue and it dominates the series for values of the renormalization scale below 200 GeV or so.The perturbative series has a systematic sign alternation which is not observed in QCD.In QCD, the coefficients can have different signs but no systematic alternation is observed, which indicates a weaker UV renormalon, as observed in other processes [30][31][32], and a more complicated interplay between UV and IR contributions.
The perturbative series in the large-β 0 limit has coefficients of roughly the same order of magnitude as in full QCD up to O(α 3 s ), or even O(α 4 s ) if we compare with the partial results of Ref. [16], although the signs of the coefficients are opposite.Assuming this observation survives at even higher orders, we estimate the six-loop coefficient to be A (5) t ≈ 300 ± 300, for µ = m t (m t ).This, in turn, leads to a contribution to the decay width of at most 0.00015 keV, which is significantly below the parametric uncertainty induced by m H and α s and much smaller than the available estimates of the O(α 4 s ) contribution [16].

A Modifications to MATAD's source code
In this appendix we briefly describe the technical implementation in MATAD of the calculation of the diagrams of Fig. 1.The implementation of the large-β 0 limit calculation in MATAD requires some modifications to the source code.The first modification deals with the Borel transformed dressed gluon propagator, where the variable u of Eq. ( 9), which is not supported in the original package, must be included.This basically generalizes the exponent of the k 2 term in the denominator of the gluon propagator to the real domain.The second modification takes into account the inclusion of the variable u into the 2-loop tadpole integrals, where the modified gluon propagator enters the calculation.
The first modification is to the gluon propagator to account for Eq. ( 9).We arranged the diagrams such that the gluon propagators always carried momentum p 3 , without any external momenta.With this, the modification is rather simple, and we only need to multiply the original propagator by a new function that we call Denu, which accounts for the factors of u in the modified propagator.In standard FORM notation, this amounts to id Dg(?x) = Dg(?x)*Denu(u); In the next step, one identifies the powers in the quark and gluon propagators and unifies the scalar integrals into a single function with the relevant coefficients in order to make contact with the analytic calculation of the relevant two-loop integrals.As discussed previously, p 3 was reserved for the gluon propagator; the momenta p 1 and p 2 were assigned to the top propagators.In the following expression, s1m and s2m represent the quark propagators with momentum p 1 and p 2 , respectively.In FORM code, the unification of the scalar integrals into a function with the relevant coefficients reads id s1m^a1?*s2m^a2?/p3.p3^a3?= f(a1,a2,a3); After this, we simply exclude the scaleless integrals which integrate to zero in dimensional regularization.
Since we are interested in the expansion in τ t for non-singlet diagramsx, we need to consider only the hard region of the loop integrals.In this region the integrals become massive tadpoles.We can then perform the corresponding scalar loop integrals using The u variable is attached to the massless propagator, i.e., in our case it appears together with the a 3 variable.We implement this with the following id statement: and analogously for the inverse of the Γ function.
B Leading top-mass corrections in the large-β 0 limit Here we provide the explicit results for the leading-n f terms up to O(τ t ) at 6 and 7 loops, which were omitted in Eq. ( 17) and (18) for the sake of brevity.The complete expressions are given below:

Figure 1 :
Figure 1: Sample of diagrams for the calculation of the leading-n f terms for the decay H → γγ.The internal dashed lines represent the resummed gluon propagator, Eq (9).

Figure 2 :
Figure 2: (a) Perturbative series, Eq. (2), in the large-β 0 limit for four different values of the renormalization scale µ, order by order.The Borel sum of the series (see text) is shown as the horizontal dashed line.(b) Perturbative series as a function of the renormalization scale up to N 5 LO (or O(α 5 s )).The scale variation range 50 GeV ≤ µ ≤ 350 GeV contains the interval m H /2 ≤ µ ≤ 2m t .