Wilson coefficients for Higgs boson production and decoupling relations to Oαs4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^4\right) $$\end{document}

An important ingredient for the calculation of Higgs boson properties in the infinite top quark mass limit is the knowledge of the effective coupling between the Higgs bosons and gluons, i.e. the Wilson coefficients CH and CHH for one and two Higgs bosons, respectively. In this work we calculate for the first time CHH to four loops in a direct, diagrammatic way, discussing in detail all issues arising due to the renormalization of operator products. Furthermore, we also calculate the Wilson coefficient CH for the coupling of a single Higgs boson to gluons as well as all four loop decoupling relations in QCD with general SU(Nc) colour factors. The latter are related to CH and CHH via low-energy theorems, which are used to obtain five-loop results for the Wilson coefficients.


Introduction
In the coming years the production and decay of Higgs bosons will play a central role in many analyses performed at the Large Hadron Collider (LHC). A crucial ingredient is often provided by the matching coefficients which govern the coupling of Higgs bosons to gluons. The corresponding effective Lagrange density is valid in the heavy top quark limit which provides a good approximation for Higgs boson decays to gluons and the total production cross section of a single Higgs boson. For less inclusive processes the applicability of the effective theory approach is limited to parts of the phase space. This is also true for Higgs boson pair production.
The straightforward way to compute the matching coefficients for the coupling of gluons to one and two Higgs boson (denoted by C H and C HH , respectively) considers the corresponding Green's functions with two gluons and one or two Higgs bosons as external particles. However, C H and C HH can also be computed using low-energy theorems (LETs). They relate the matching coefficients to the decoupling constant relating the strong coupling constant defined in five-and six-flavour QCD, ζ αs . For the case of C H the LET has been derived in ref. [1], where an essential ingredient is the effective Lagrangian containing a complete basis of all scalar dimension-four operators. On the other hand, for C HH a LET has only been proposed very recently in ref. [2] without providing a strict derivation. Note that both for C H and C HH the corresponding LET can be used to obtain (n + 1)-loop results from the n-loop expression for ζ αs .

JHEP11(2018)141
Two-loop results for C H have been known since the beginning of the eighties [3,4] and at three-loop order C H was obtained for the first time from a direct calculation of the Higgs-gluon vertex in the large-m t limit in ref. [5] (see also ref. [6]). Later the result was confirmed with the help of the LET [5] (see also ref. [7]). Using the three-loop decoupling constant for α s , the LET in combination with the four-loop beta function [8,9] even leads to the four-loop result for C H . The same reasoning has been applied in refs. [10][11][12] to obtain the five-loop prediction for C H , where an important input is provided by (the fermionic part of) the five-loop beta function which was computed in refs. [13][14][15]. To date there is no direct calculation of C H at four loops. Note that the existing four-and five-loop results are only available for SU (3).
For C HH the situation is as follows: at one-and two-loop order C HH and C H are equal. At three-loop order a direct calculation has been performed in ref. [16] by matching the full theory to the effective theory (see eq. (2.5) below). The result has been confirmed via the LET from ref. [2], which can be used to derive the four-loop result for C HH . 1 The main purpose of this paper is the direct calculation of C HH to four-loop order, thus providing an independent check for the LET from ref. [2]. In order to gain confidence in our computer programs we first compute the matching coefficient C H at four-loop order by considering the Higgs-gluon-gluon vertex. We compare our result to the LET expressions in the literature [5,[10][11][12]. All of our results are expressed in general SU(N c ) colour factors.
The essential ingredient for the LETs is the QCD decoupling constant for the strong coupling. For this reason we re-visit the calculation of all four-loop decoupling constants and provide results for a generic SU(N c ) gauge group. Currently the four-loop expression for ζ αs is only available for N c = 3 from refs. [10,11]. The result for ζ αs can be used to obtain C H and C HH to five-loop order, again for general SU(N c ) colour factors.
The remainder of this paper is organized as follows: in section 2 we fix our notation and introduce the decoupling constants and the effective Lagrange density for the Higgs-gluon coupling. In sections 3 and 4 we present our results for the decoupling relations and Wilson coefficients, respectively. We discuss in detail the extraction of the coupling of two Higgs bosons to gluons (C HH ) and in particular the subtleties in the matching procedure due to the renormalization of products of operators. Our findings are summarized in section 5. In the appendix we collect analytic results for the decoupling constants.

Technicalities
For convenience of the reader and to fix our notation we repeat in this section the definition of the decoupling constants in QCD and the Wilson coefficients in the effective Lagrange density describing Higgs-gluon couplings. For a detailed discussion we refer to ref. [1]. We will work in the MS scheme throughout this paper, except for the heavy quark mass which we renormalize both in the MS and on-shell scheme. The MS counterterms are needed up to four-loop order (see, e.g., ref. [17]) and the renormalization constant for the MS to on-shell conversion for the heavy quark mass to three loops [18][19][20][21].

Decoupling constants
The bare and renormalized parameters and fields of the QCD Lagrangian are connected by renormalization constants defined through Here g s is the QCD gauge coupling with α s = g 2 s /(4π) being the strong coupling constant, µ is the renormalization scale, D = 4 − 2 the space-time dimension and ξ the gauge parameter with ξ = 0 corresponding to the Feynman and ξ = 1 to the Landau gauge. The gluon field is given by A a µ , ψ q is the quark field of flavour q with mass m q and c a is the ghost field. Bare quantities are denoted by the superscript "0". The renormalization constants Z X are needed up to O(α 4 s ) [8,9,17,22,23] for our purposes. In the following we assume a strong hierarchy in the quark masses and integrate out a heavy quark with mass m h from QCD with n f active quark flavours. 2 The resulting effective Lagrangian has the same form as the original QCD Lagrangian. However, it only has n l = n f −1 active quark flavours and thus only depends on the light degrees of freedom. The parameters and fields in the effective n l -flavour and full n f -flavour theory are related via the so-called (bare) decoupling constants where the bracketed superscripts denote the number of active quark flavours. For simplicity we refrain from showing the colour indices for the fields. The different decoupling constants ζ X contain the radiative effects of the heavy quark and can be computed in a perturbative series in α s . One obtains the renormalized decoupling constants by replacing the bare parameters and fields in eq. (2.2) by renormalized counterparts using eq. (2.1). As an example, consider the gauge coupling where the renormalized decoupling constant is given by Note that Z (n l ) g depends on g (n l ) s which has to be transformed to g (n f ) s using eq. (2.3). Thus, it is natural to apply eq. (2.3) iteratively, order by order in g s , to arrive at four loops. Note that the loop corrections to the renormalization constants only contain poles whereas the decoupling constants also contain positive powers in . Thus, Z

Wilson coefficients for Higgs boson production and decay
In the Standard Model, the coupling of a Higgs boson to gluons is mainly mediated by top quark loops and thus in the following we have n f = 6 and n l = 5 for the full and effective theories, respectively. The effective Lagrange density which describes the coupling of one or two Higgs boson to gluons is obtained after integrating out the top quark and is given by where O 1 = G a µν G µν,a /4, with G a µν being the gluon field strength tensor, is the only physical operator one has to consider. It is defined in the (n l = 5)-flavour effective theory. The Wilson coefficients C 0 H and C 0 HH comprise of the radiative effects of the top quark, which is in analogy to the decoupling constants introduced in eq. (2.2).
The renormalization of O 1 has been discussed in detail in ref. [25] (see also ref. [26]). In fact, the renormalization constant Z O 1 can be expressed through the QCD beta function through all orders in perturbation theory [25] Z O 1 can be used to obtain the renormalized Wilson coefficients via the relation with X ∈ {H, HH}.

Low energy theorems
There is a close connection between the decoupling constants from subsection 2.1 and the Wilson coefficients from 2.2 which is established by the so-called LETs. In [1] a LET relating ζ αs and C H has been derived

JHEP11(2018)141
where we have adapted the prefactors to match our conventions. Beyond three loops C H has only been obtained by using eq. (2.9). In this work we perform an explicit calculation of C H for general SU(N c ) colour factors, to four-loop order.
Recently, in ref. [2] a LET has been proposed for C HH which reads (2.10) It provides the correct result for C HH at three-loop order [16]; in section 4 we perform an explicit calculation of C HH and show that eq. (2.10) also works at four loops. Note that in QCD, ζ αs depends on m t only via logarithms of the form log(µ 2 /m 2 t ). Thus, it is possible to reconstruct the m t dependence of C HH at (n + 1)-loop order from the n-loop result, by using the renormalization group equations. Using the LETs this immediately leads to the (n + 1)-loop results for C H and C HH .

Computational setup
For our calculation we use a well tested, automated setup, starting with the generation of Feynman diagrams using qgraf [27]. The output is processed by q2e and exp [28][29][30], which generate FORM [31] code for the amplitudes and map the diagrams onto individual integral families. We then compute the colour factors of the diagrams using color [32] and combine amplitudes with the same colour factor and integral family from so-called superdiagrams, so that we can process them together.
After processing Lorentz structures and expanding in the external momenta, we are left with single-scale tensor tadpole integrals. We perform a tensor decomposition and reduce the remaining, scalar integrals to master integrals, using LiteRed [33,34] and FIRE5 [35]. With the help of the FindRules command of FIRE5 we identify equivalent master integrals from different integral families.
The master integrals are all known to sufficiently high order in [36] (see also [37,38]). The missing 3 term of the integral J 6,2 (in the notation of [36]) was provided in [39].
As a cross-check, we also computed the ggH amplitude and the decoupling constants to three loops using MATAD [40].

Calculation of decoupling constants
We aim to calculate all QCD decoupling constants up to four-loop order with general SU(N c ) colour factors. They are obtained from ζ 0 3 ,ζ 0 3 , ζ 0 2 and ζ 0 m as introduced in eq. (2.2) and the decoupling constant of the ghost-gluon vertex,ζ 0 1 . The decoupling constant for the gauge coupling is then given by Γ Gcc -5 228 10 118 Table 1. Number of diagrams contributing to the decoupling constants up to four loops.
The remaining decoupling constants for the gluon-quark vertex (ζ 0 1 ), the three-gluon vertex (ζ 0 3g ) and four-gluon-vertex (ζ 0 4g ) are obtained with the help of the Ward identities The bare decoupling constants ζ 0 3 ,ζ 0 3 , ζ 0 2 and ζ 0 m are obtained from the hard part of the gluon and ghost vacuum polarizations Π G (p 2 ) and Π c (p 2 ), as well as the vector and scalar parts of the light quark self-energy Σ V (p 2 ) and Σ S (p 2 ), as [1] To obtain the decoupling constants we only need the leading term in the limit m h → ∞. Thus, we can Taylor expand in the external momenta and set them to zero after factoring out the tree-level tensor structure. This reduces the integrals to single-scale tadpole integrals. In analogy to eq. (3.3),ζ 0 1 is obtained from the ghost-gluon vertex where p and q are the four-momenta of the ghost and gluon, respectively. After projecting out the tree-level contribution both p and q are set to zero. In table 1 we present the number of diagrams generated by qgraf for the individual Green's functions. Sample four-loop Feynman diagrams are shown in figure 1. We perform the calculation keeping the full dependence on the gauge parameter ξ which drops out for ζ 0 αs and ζ 0 m , as expected on general grounds. All other decoupling constants have an explicit ξ dependence. At three-loop order our results agree with those of ref. [1] and at four loops we reproduce the results for ζ αs from refs. [10,11] and ζ m from [41] after specifying N c = 3. The results for ζ αs and ζ m as well as the leading terms of the other decoupling constants can be found in appendix A. We provide the results for all renormalized decoupling constants in computer readable form in the supplementary material attached to this paper. For convenience we offer several options concerning the renormalization scheme of the heavy quark (MS vs. on-shell) and α s (n f vs. n l active flavours).

Direct calculation of the matching coefficients
This section is devoted to the direct calculation of C H and C HH defined in the effective Lagrange density in eq. (2.5). We use the notation for the matching equations introduced in ref. [16] and compute the ggH and ggHH amplitudes in the limit where both the effective and the full theory are valid, i.e. for small external momenta as compared to the top quark mass. This leads again to single-scale vacuum integrals up to four-loop order. In the following we use α s ≡ α (5) s (µ) if not otherwise indicated.

C H
The Wilson coefficient is obtained by comparing the ggH amplitude in the effective and full theory which leads to the following matching formula On the full-theory side A h denotes the hard part of the amplitude, which is obtained from a Taylor expansion in the two external momenta. It is assumed that the top quark mass and α s are renormalized using standard counterterms up to three loops and the factor 1/ζ 0 3 takes care of the non-vanishing part of the gluon wave function renormalization. Due to our choice of the kinematic variables there are only tree-level contributions on the effectivetheory side. Additionally, we have the renormalization constant of the effective operator, Z O 1 , and the sought-after (renormalized) matching coefficient C H , which is obtained by dividing eq. (4.1) by s whereas the quantities on the r.h.s. depend on α (6) s . Before combining the various parts we use the decoupling constant to transform the strong coupling constant on the r.h.s. to α (5) s . We first renormalize the top quark mass in the MS scheme and transform to the on-shell scheme afterwards.  The number of diagrams generated by qgraf for A h is shown in table 2 and sample Feynman diagrams are shown in figure 2. We begin by applying the projector where q µ 1 and q ν 2 are the incoming four-momenta of the external gluons with polarization vectors ε µ (q 1 ) and ε ν (q 2 ). After tensor reduction we obtain the same kind of integral families as for the decoupling constants of the previous section.
As before, we perform the calculation for generic SU(N c ) colour factors and full dependence on the gauge parameter ξ, which drops out after summing all contributions to A h . We cast the final result for the Wilson coefficient C H in the form where the C (i) are given by (4.7) ζ(n) is the Riemann ζ-function, evaluated at n, M t is the on-shell top quark mass and the SU(N c ) colour factors are given by

C HH
The matching procedure to obtain C HH is more involved than that of C H . First of all there are three contributions on the effective-theory side which are shown in figure 3: a one-particle irreducible (1PI) term proportional to C HH , a one-particle reducible (1PR) term, which involves C 2 H , and a term mediated by a virtual Higgs boson which splits into a Higgs boson pair via the Higgs boson self-coupling λ. The latter is similar in nature to the effective amplitude in the matching formula for C H . In fact, also on the full-theory side this contribution involves diagrams which we already encountered in the computation of C H . As mentioned in ref. [16] it is easy to see, that these diagrams exactly cancel between the full and effective theory. Thus, the contributions relevant to extract C HH are the 1PI and 1PR contribution with λ = 0.
The effective-theory side of the matching formula is obtained after renormalizing the operators in the various contributions of figure 3. Whereas the left and right contributions of figure 3 are both renormalized with a factor Z O 1 , the term in the middle needs special care. In fact, a naive renormalization with (Z O 1 ) 2 leads to uncanceled poles as has already been observed in refs. [42,43]. A careful analysis of the renormalization of the product of two operators O 1 has been performed in ref. [26] along the lines of [25]. It has been observed that apart from the naive multiplicative renormalization a further term is needed which is proportional to a single O 1 . Adapting the findings of ref.
[26] to our notation one has where A eff It has its first non-vanishing contribution at order α 2 s . As we will see below, in our calculation we need the combination Z L 11 /Z O 1 up to order α 2 s which is given by   We are now in the position to write down the matching formula for C HH . Equating the effective-theory side, which is basically given by figure 3, with the corresponding full-theory amplitudes and taking into account eq. (4.9) leads to 3 where sample Feynman diagrams contributing to A h 1PI and A h 1PR,λ=0 can be found in figures 4 and 5, respectively. As already mentioned above, the contributions with λ = 0 cancel in eq. (4.12). Note that our matching formula differs from the one of ref. [16] by the term proportional to Z L 11 which contributes for the first time at four-loop order, since both C 2 H and Z L 11 are of order α 2 s . Let us in the following discuss some features of the matching procedure. At one-loop order the only non-zero contribution on the r.h.s. of eq. (4.12) is A h 1PI and one obtains C (1) H . This also holds at two-loops where the 1PR contributions on effectiveand full-theory side match exactly. A non-trivial interplay between A h 1PI and A h 1PR,λ=0 is observed for the first time at three-loop order [16]. In fact the 1PI and 1PR contributions are not separately finite any more and the poles only cancel in the sum. Starting from this order C HH is different from C H . While the 1PI and 1PR contributions are separately ξ-independent at three loops, for the four-loop colour structure C 2 A T F it only drops out in the proper combination.
We computed the 1PI and 1PR amplitudes in the full theory separately and keep in both cases terms linear in the gauge parameter ξ. For both contributions it is important 3 When applying eq. (4.9) to Higgs boson pair production we have A eff

JHEP11(2018)141
to keep the three external momenta different from zero and different from each other in order to avoid the mixing with unphysical operators [25]. The external momenta can be set to zero after projection to the matching coefficient which is done with the help of where q ij = q i · q j and q 2 T = 2q 13 q 23 /q 12 − q 33 . q µ 1 and q ν 2 are the incoming four-momenta of the external gluons with polarization vectors ε µ (q 1 ) and ε ν (q 2 ) and q 3 is the incoming four-momentum of one of the Higgs bosons.
The number of diagrams for the 1PI amplitude can be found in table 2 and sample diagrams are shown in figure 4. Once the projector of eq. (4.13) is applied one obtains scalar expressions which still contain scalar products of q 1 , q 2 and q 3 and loop momenta in the numerator. After solving the corresponding tensor vacuum integrals the resulting scalar products q ij cancel against the corresponding contributions with negative powers from the projector and all external momenta can be set to zero.
The 1PR amplitude has been obtained in two different ways. First, we computed the 1PR diagrams up to four-loop order (see table 2 for the number of diagrams and figure 5 for typical Feynman diagrams) in analogy to the 1PI contribution. As a cross-check we computed the 1PI parts of the 1PR contributions separately and constructed the n-loop 1PR ggHH amplitude from ggH amplitudes computed up to (n−1) loops. In this approach one of the gluons in the ggH amplitude has to be off-shell, which leads to more nonvanishing Lorentz structures. In practice, we computed the 1PI ggH amplitudes with open Lorentz indices up to three loops. Full agreement has been found between the two methods.
We cast the final result for the Wilson coefficient C HH in the form , (4.14) where the C (i) H are given in eq. (4.7) and the differences are given by ∆ (1) (4.15)

JHEP11(2018)141
The three-loop result can be found in ref. [16]. Our four-loop result ∆ (4) HH agrees with the expression from eq. (2.10) [2]. We can thus confirm the validity of the LET for C HH [2] through four loops. In analogy to C H also for C HH it is possible to construct the fiveloop approximation for general colour structure. The corresponding results can be found in computer readable form in the supplementary material attached to this paper. After setting N c = 3 we agree with the numerical results given in ref. [2], both for MS and on-shell top quark mass.

Conclusions
We perform for the first time a direct four-loop computation of the Wilson coefficients C H and C HH of the effective operators which couple gluons to one and two Higgs bosons, respectively. C H and C HH enter, as building blocks, various physical quantities, e.g., the next-to-next-to-next-to-leading order predictions for single [44,45] and double Higgs boson production. 4 Our results for C H and C HH agree with the expression obtained by means of LETs. Furthermore, we compute all QCD decoupling constants up to four-loop order. Where possible we compare with the literature and find agreement after specifying the colour factors. All our results are expressed for general SU(N c ) colour factors whereas the four-loop expressions in the literature are only available for N c = 3.
A major result of this paper is the derivation of the matching equation (4.12) which receives a non-trivial renormalization contribution from the effective-theory amplitude with two insertions of the operator O 1 . The new term contributes for the first time at four-loop order and is essential to obtain a finite result.
For the convenience of the reader we collect all analytic results obtained in this paper in the supplementary material attached to this article.

A Decoupling constants
In this appendix we collect the results for the decoupling constants for general SU(N c ) colour factors. We provide results for ζ αs and ζ m up to four loops and show for all other ζ constants the expressions for the first non-vanishing loop-order. Computer readable expressions up to four loops can be found in the supplementary material attached to this paper. Our results read

JHEP11(2018)141
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.