Quark mass effects in two-loop Higgs amplitudes

We provide two two-loop amplitudes relevant for precision Higgs physics. The first is the two-loop amplitude for Higgs boson production through gluon fusion with exact dependence on the top quark mass up to squared order in the dimensional regulator $\epsilon$. The second result we provide is the two-loop amplitude for the decay of a Higgs boson into a pair of massive bottom quarks through the Higgs-to-gluon coupling in the infinite top mass limit. Both amplitudes are computed by finding canonical bases of master integrals, which we evaluate explicitly in terms of harmonic polylogarithms. We obtain the bare, renormalized and IR-subtracted amplitude and provide the results in terms of building blocks suitable to changing renormalization schemes.


Introduction
The Run II of the LHC has allowed experimental collaborations to probe the Higgs boson to unprecendented levels of precision. Recent combination results by CMS [1] and ATLAS [2] show a 60% reduction in the global signal strength compared to the historic Run I combinations. Run II has also seen impressive developments in differential observables [3,4] which can provide rich information on the dynamics of the Higgs boson. As a result of this steady experimental progress, making precise theoretical predictions of the relevant observables is highly important as their comparison to measurements will allow us to test the Standard Model and highlight possible new physics through any discrepancy.
Intense theoretical effort has been devoted to making such highly precise predictions of Higgs observables. The inclusive cross section was obtained at next-to-next-to-next-to-leading order (N 3 LO) in the QCD coupling using the Higgs Effective Theory (HEFT) where the top quark is infinitely massive [5]. Approximated Higgs-boson rapidity distributions were also obtained at N 3 LO [6,7].
The infinite top mass approximation has a ∼ 6% effect on the cross section, estimated from the NLO [8] prediction, which is applied to the state-of-the-art N 3 LO through a multiplicative correction factor. The finite-mass corrections mostly factorize from the perturbative corrections [5], so that rescaling by known exact results induces only an estimated ∼ 1% uncertainty on the prediction [9,10]. At the inclusive level, the uncertainty associated to this rescaling consitutes therefore a sizeable portion of the ∼ 5% theoretical error. At the differential level, mass effect are all the more important in the high energy region where the HEFT has been shown to fail by the first exact NLO prediction of the Higgs boson tranvsere momentum (p T ) [11]. While this work also highlighted that more refined approaches such as the so-called FT approx description can provide a reasonable description within 10% up to high energies, the projected ∼ 5% uncertainty of future HL-LHC transverse momentum spectrum measurements [12,13] warrants turning our sights toward a better control of mass effects in Higgs physics predictions. This situation could be improved by computing the NNLO hadronic Higgs boson cross section including exact top-mass effects.
This goal is becoming realistic thanks to the recent derivation of the three-loop doublevirtual contribution, which started with the light-fermion contributions [14] and was completed by the numerical derivation of the complete result [15]. Combined with the knowledge of all integrals of the two-loop real-virtual contributions [16], a full prediction is now within reach.
Motivated by this situation, the first part of this paper provides the analytic result for the two-loop amplitude for the process gg → H to order O ε 2 , which is required to build the infrared (IR) subtraction terms of the double-virtual contribution. These were used in [15] to obtain a finite remainder but are not publicly available.
Probing rare production channels beyond the dominant top-loop mediated process is in-strumental to a comprehensive study of the interactions of the Higgs boson with other Standard Model particles. Weak production modes such as Higgstrahlung (VH) provide key insights to test our understanding of electroweak symmetry breaking as well as a window into the Higgs to bottom quark decay channel, which is otherwise dominated by QCD backgrounds [17,18]. The current uncertainties do not qualify this process as a precision observable, but statistics will significantly improve the situation [19] to a point where theory uncertainties are expected to dominate. The current theoretical state-of-the-art predictions [20,21] combine NNLO predictions for the production [22,23] and the decay [24,25] fully differentially. This work has shown that even the pure NNLO correction to the decay can have large effects on differential observables, which motivates improving our description of the decay of the Higgs boson to a pair of bottom quarks (H → bb).
The current state-of-the-art prediction for the H → bb decay is N 4 LO at the fully inclusive level [26] and N 3 LO at the differential level [27]. These predictions are made in the limit where the bottom quark is massless and therefore neglect contributions from top-quark loops induced by the top-quark Yukawa coupling. These appear at NNLO and generate difficult to treat infrared divergences in the massless bottom-quark limit [21]. This difficulty means that the state-of-the-art predictions for V H, H → bb [20,21], which rely on massless bottom NNLO calculations [24,25] also miss these contributions. Top-induced effects are currently untractable in massless bottom calculations, so that they can only be obtained by performing the calculation of the Higgs decay into massive bottom quarks. This was completed at NNLO recently [28] and included a HEFT description 1 of the first non-zero top quark effects. This work highlighted that the top-induced contributions have an impact of around 2% on the Higgs width through their interference with the leading process and therefore contribute about 25% of the pure NNLO effects. In order to further improve our control of the H → bb decay, it is desirable compute the top-induced N 3 LO effects in the HEFT. This is the first order at which squared top-induced processes occur, which we can straightforwardly compute using automated tools such as MG5 aMC@NLO [30]. At the inclusive level, these squared contributions have an effect of about 1% on the decay width, making them dominant over the existing N 3 LO prediction, which are of around 0.2% [26], motivating the derivation complete the N 3 LO topinduced contributions. The missing piece of this calculation is the two-loop amplitude for the decay of a Higgs boson to a pair of massive bottom quarks mediated by the Higgs-to-gluon coupling in the HEFT, which we provide in the second part of this paper. This paper is organized in the following way. Section 2 presents the calculation of the one and two-loop amplitudes for gluon-fusion Higgs production to high order in the dimensional regulator and provides the results and their expansions in two kinematic limits. Section 3 presents the same results for the one and two-loop amplitudes that contribute to the top- 1 The HEFT description of top-induced H → bb was found to be extremely accurate by comparing to the exact calculation [29] quark-Yukawa-induced Higgs to bottom decay. We subsequently discuss the analytic continuation of the result in section 4 and finally discuss the details of the computation of the master integrals (MIs) in section 5.
2 Amplitudes and results for gg → H

Notation for bare amplitudes
The bare amplitude A 0 gg→H of the process g(p 1 )g(p 2 ) → H can be written as and v 0 denotes the bare vacuum expectation value. The gluons are in physical gauge with ε i (p i ) · p i = 0. In order to compute the form factors M 0 X , we first generate all contributing diagrams with QGraf [31] and perform the color-, Dirac-and Lorentz algebra in Mathematica. Traces of γ matrix chains are performed with FORM [32]. The integration-by-parts (IBP) reductions [33,34] to scalar master integrals (MIs) are done with the programs AIR [35] and Kira [36].
We separate M 0 N LO according to where infrared singularities are isolated in M 0 IR and ultraviolet poles are contained in M 0 UV and M 0 UV,m , respectively harboring terms renormalized by coupling and mass counterterms. Of the two remaining regular terms, M 0 fin,scale contains the complete dependence on the renormalization scale µ 2 while M 0 fin corresponds to the case of µ 2 = m 2 H . We have

5)
We checked the bare LO and NLO against [8] by inserting our expressions into eq (7.4) therein. We furthermore compared the ε 0 of M 0 NLO against [37] and the analytic expression implemented in the program iHixs 2 [38] 2 . We find full agreement in all cases 3 . The higher orders in ε of the amplitude (2.1) are very cumbersome. We therefore only report their expansion in kinematics limits where they simplify in this manuscript and provide the exact results as ancillary files.

Renormalization and IR-subtraction
The renormalized amplitude reads where Z g denotes the gluon wave-function renormalization function, µ the renormalization scale and the superscript 0 indicates bare quantities. The renormalized parameters are related to the bare ones by: and we define the bare Yukawa coupling by its relation to other parameters: y 0 t = m 0 t /v 0 . We renormalize the strong coupling and the gluon field in a mixed scheme with N f = 5 light flavours, whose contributions are subtracted in MS while contributions involving the top-quark are renormalized on-shell, at zero momentum [41]. This yields and (2.11) For the sake of more compact expressions we renormalize the top mass in MS: (2.12) 2 see: https://github.com/dulatf/ihixs/blob/master/src/higgs/exact_qcd_corrections/nlo_exact_ matrix_elements.cpp and function ggf exact virtual ep0 therein. 3 The comparison with iHixs 2 requires the subtraction of IR divergences with β 0 ε + 2Nc ε 2 (− s µ 2 ) −ε as defined in e.g. [39] or [40] removing the β0 dependence in M 0 fin .
Note that these choices are such that the counterterms generated by eq. (2.11) and eq. (2.12) cancel the UV terms in the bare amplitude of eq. (2.3) up to neglected orders in α s but to all orders in ε: where we exploit the fact that m 0 t − m t = O (α s ) can be neglected at this order. The renormalized amplitude still features poles in ε that are of infrared and collinear origin. These singularities have a universal structure in that it can be expressed in a factorized fashion using lower orders in the amplitude [39,40]: where F is finite and M LO is the renormalized leading-order scalar amplitude, which in our case is trivially obtained by replacing m 0 t by m t in M 0 LO . Again, our splitting of the bare amplitude in eq. (2.3) is such that the IR subtraction term cancels M 0 IR to all orders in ε and to all relevant orders in α s so that where M fin,scale and M fin are obtained from M 0 fin,scale and M 0 fin by substituting m 0 t with m t . The complete renormalized and IR-subtracted NLO-contribution to gg → H in the above discussed schemes is simply given by The artificial splitting of the bare amplitude in (2.1) and the corresponding ancillary material is designed to make changes of renormalization or IR-subtraction schemes particularly simple. A change of renormalization schemes, e.g. to the on-shell scheme for the top mass renormalization with can straightforwardly be obtained by computing the corresponding finite piece 19) and adding it to the MS renormalized NLO piece in (2.17) obtaining

Kinematic limits
In the following we discuss the amplitude in the limits |s| m 2 t and |s| m 2 t . The second limit in particular can be used as an important check of the full result since it has a direct correspondence to the heavy top EFT, in which the inclusive cross-section of gg → H is known to N 3 LO [5, 42].

Small mass expansion
We perform the expansion around the limit m t → 0 (or |s| → ∞) with the code PolyLog-Tools [43] in the Euclidean regime and find for the leading order contribution The NLO pieces in (2.1) expanded in the limit of a small top-mass are

Large mass expansion
The large mass limit m t → ∞ (or |s| → 0) of the amplitude eq. (3.8) can easily be obtained as an all order expression in the dimensional regulator ε, by employing the method of regions. It therefore can be used as a non-trivial check of the higher order terms of M 0 NLO in eq. (2.3).
The LO-amplitude for the limit m t s is obtained by expanding on the integrand level. We find the all orders expression The NLO amplitude eq. (2.3) in the limit m t → ∞ factorizes as Where we separate the terms that correspond to different regions (in the sense of expansions by regions). The first term in (2.26) corresponds to the region hard soft region m t ∼ k 1 k 2 , p 1 , p 2 where k 1 denotes the top-loop momentum and the p i are external momenta. C LO ε is given in (2.25) and is the one-loop contribution to gg → H in the heavy top EFT (see e.g. eq (3.5) in [44]). The second term in (2.26) corresponds to the double hard region m t ∼ k 1 ∼ k 2 p 1 , p 2 . These all-order results in the dimensional regulator ε were obtained by expanding the momentum space representation of the loop integrals in each region and directly integrating the result.
On the other hand, we can expand our results for the bare amplitudes (2.1) in terms of harmonic polylogarithms in the limit m t → ∞ with the help of HyperInt and PolyLogTools. Including the normalization we have The result agrees with (2.25) expanded to O ε 4 , which is an important check of our computation.
For the two-loop pieces we find by directly expanding the result eq. (2.3) in terms of harmonic polylogarithms in the large top-mass limit and in particular 3 Amplitudes and results for H → bb

Higgs effective field theory
The Higgs Effective Field Theory (HEFT) is obtained by integrating out the top quark from the Standard Model [45]. In practice, as long as we do not describe electroweak corrections, it is equivalent but much simpler to describe our calculation in the context of QCD coupled to a singlet scalar H with mass m H 4 , yielding the following bare Lagrangian for the EFT: /v B and C 2 is the HEFT correction factor to the bottom Yukawa: when matching the HEFT to the SM We leave unspecified the details of the gauge-fixing and ghost Lagrangian of the gauge interaction L gf . The couplings mediated by C 1 and C 3,...,6 correspond to the next to leading power terms in the expansion of the exact Lagrangian in powers of the top mass, which we only show explicitly for C 1 since the other operators do not contribute to on shell amplitudes [45]. Note that due to the absence of a Higgs mechanism in our UV-complete theory, the top quark mass and Yukawa are not necessarily related so that we can consistently distinguish power counting in y t and m t . Consequently, C 1 is labelled as next-to-leading power despite being non-decoupling when matched to the full SM, where C 1 ∝ y t /m t .
We provide our results expressed in terms of HEFT parameters exclusively, leaving the matching to SM parameters to future applications. As a result, and for the sake of readability, we will drop explicit HEFT labels in the couplings and masses in the rest of this section, as we never refer to SM parameters.

Notation for bare amplitudes
The bare amplitude A B H→bb of the process H → b(p 1 )b(p 2 ) can be written to all orders as where M B is a Dirac matrix and i, j are color indices of the fundamental representation of SU (3). The external kinematics obey where both m 2 b and s are finite numbers (i.e. the p 2 i are not the bare masses). M B can be decomposed as 5 and we will therefore restrict further discussions to the scalar quantity M B 0 , which is obtained with the techniques discusses in section 2.1. We furthermore define the following decomposi- is generated diagrammatically by replacing the bottom quark propagator by its derivative, which we indicated with a red cross in eq. (3.10), i.e We leave the discussion of the renormalization of the mass and the precise definition of δm for the next section. In this section, we instead focus on defining the bare amplitudes in terms of simple components. We separate the mass-renormalized amplitude M B yt,1 (m b , m b ) according to and M B yt,2 as such that the poles are separated by IR-and UV-origin respectively. The UV-divergent contribution of the one-loop amplitude reads The two-loop UV-poles are separated according to The IR-divergences can be described using the factorization of next-to-leading order amplitudes with massive external colored particles [40] 6 as In the region where s < 0, or equivalently 0 < x < 1. We discuss the analytic continuation to the physical region s > 0 in section 4.

Renormalization and IR-subtraction
The renormalized amplitude of the process bb → H reads where the superscript B denotes bare quantities, Z b is the wave-function renormalization function of the massive b-quarks and µ the renormalization scale. The bare and the renormalized parameters are related by where the renormalization constants Z X are parametrized as The y b -renormalization is completely determined by the mass renormalization. The part of the renormalized amplitude that is proportional to C 1 reads where the scalar quantities M X (m b ) are defined in the previous section. We renormalize our amplitude in MS with the relevant renormalization constants [46][47][48] Comparing with eq. (3.13) and eq. (3.14) yields the MS renormalized amplitude We provide all contributions to the bare amplitudes M 0 yx,n as well as the renormalized and IR-subtracted amplitudes M f in yt, (1,2) in the ancillary material, such that results for a different choice of a renormalization scheme can be easily obtained (see section 2.2). We furthermore provide all necessary master integrals for this process up to weight 6, such that higher orders in the dimensional regulator are easily accessible for future computations.

Small mass expansion
In the limit where m 2 b |(p 1 · p 2 )| the renormalized and IR-subtracted amplitudes have the expansion and M fin.

Analytic continuation for gg → H and H → bb
Our results are provided in the Euclidean regime s < 0 and we discuss in the following how they can analytically be continued to the physical regime. The amplitude gg → H has the production threshold at s = 4m 2 t and the pseudo threshold s = 0. A parametrization of the amplitude in terms of the "natural" scaleless ratio s/m 2 t will yield undesirable roots of the form (4.1) The same happens for H → bb, where the scaleless variable s/m 2 b gives rise to the same roots. To make discussion valid for both amplitudes under consideration, we introduce the scaleless ratio To rationalize the roots we work with the scaleless complex variable x defined by and 0 < |x| < 1. Here Feynman's prescription is denoted by +iη, implicitly defining the branch on which to evaluate the roots in the definition of x eq. (4.3). More explicitly we have for all kinematic regions. Feynman's prescription is denoted by the explicit +iη. (4.5) The last line indicates that above threshold (s > 4m 2 q ) where −1 < x < 0, x has to be evaluated by approaching the negative real axis from the upper half plane. The variable x is shown in fig. 2.
The complete result of for gg → H at NLO as well as H → bb can be expressed in terms of harmonic polylogarithms [49] with argument x eq. (4.3). A harmonic polylogarithm of weight n is defined as the iterated integral H(a n , a n−1 , . . . , a 1 ; x) = x 0 H(a n−1 , . . . , a 1 ; t)f (a n , t)dt , (4.6) where a i ∈ {1, 0, −1} and For the case of all a n , . . . , a 1 being zero we define Harmonic polylogarithms form a shuffle algebra [49] and have a branch point at x = 0 if and only if a 1 = 0. If a 1 = 0 one can use the shuffle algebra and rewrite the HPL as a linear combination of products of HPLs such that every HPL of weight j ≤ n appearing in this linear combination which has a 1 = 0 has also a k = 0 for all k = 2, . . . , j, i.e. all discontinuities around x = 0 can be described by a polynomial in log(x). This method of explicitly extracting the logarithmic factors is implemented in several publicly available codes [38,[50][51][52], whereas we used PolyLogTools to perform this task. Once the logarithms are extracted explicitly, as in the provided ancillary material, the complete analytic continuation to the regime s > 4m 2 q is obtained by performing the limit lim η↓0 + log(x + iη) = log(−x) + iπ . (4.9) All other regimes have no subtleties and can be evaluated by using the explicit prescription in the right-hand side of eq. (4.5).

Computation of master integrals
We define a generic l-loop integral depending on the kinematic invariant s and the mass m q > 0 as where the D's denote the propagators, the ν i ∈ Z their respective powers and the normalization is chosen to render the integrals scaleless. For computing loop integrals analytically two main techniques are employed nowadays.
One is based on rewriting the integral in terms of Feynman parameters and attempting a direct integration. Powerful tools like programs HyperInt and PolyLogTools are dedicated towards performing these parametric integrals.
The second method is based on deriving a closed system of differential equations [53,54]. Instead of attempting a direct integration of the integrals, one tries solving a corresponding system of coupled first order differential equations obtained by taking derivatives with respect to all external and internal scales s i . In [55] it was conjectured that for a large class of Feynman integrals a particular basis choice of MIs can be found, such that the dependence of the dimensional regulator factors out completely. For such a canonical basis the total differential takes the particular simple form where I n denotes nth Laurent coefficient of the integrals and the matrix A depends on the external and internal scales s i only. A formal, general solution of the system of differential equations eq. (5.2) for every Laurent-coefficient in the ε-expansion of the canonical integrals can be written down in terms of Chen iterated integrals [56] directly. If the entries of the matrix A are Q-linear combinations of logarithms one can often find a solution in terms of multiple polylogarithms [57] defined for a k = 0 by the iterated integral For the special case where all a i ∈ {1, 0, −1} they reduce to harmonic polylogarithms defined in eq. (4.6), which appear in the amplitudes under consideration. The one-loop contribution M 0 LO in A 0 gg→H eq. (2.1) gives rise to one integral family shown in fig. 3 which has three MIs. We choose the following basis

Master integrals for
and we report all necessary coefficients to compute A 0 gg→H eq. (2.1) to O ε 2 in the ancillary files. The computation of the master integrals was done by performing the Feynman parameter integral with the help of the program HyperInt. For the computation of M 0 uv,m eq. (2.5), corresponding to the mass renormalization, we need the mass derivatives of the one-loop MIs. Since we defined a canonical basis, they take the particular simple form: where f k,n i denotes the nth Laurent coefficient of the kth MI.  All scalar integrals of the complete two-loop form factor M 0 N LO eq. (2.3) can be written in terms of integrals of the three auxiliary families a, b and c listed in table 1. In order to compute this NLO contribution to gg → H with full mass dependence to O ε 2 a total of 18 two-loop master integrals (MIs) have to be computed to higher orders in the dimensional regulator than known in the literature [8,[58][59][60]. The topologies corresponding to the MIs are depicted in fig. 4, where the dashed line corresponds to the Higgs, wavy lines denote massless and continuous straight lines massive propagators. The ν i denote the ith propagator and the superscripts A, B, and C denote the corresponding scalar family. As a canonical basis of integrals we chose the set eq. (5.8). The corresponding topologies appear already as a subset of integrals in [62] and we deviate from their choice of MIs only slightly.

Master integrals for
We derive the differential equation using LiteRed, [63] perform the necessary IBP-reduction with Kira and integrate the differential equation order-by-order in ε. As a boundary point we consider the point x = 1 corresponding to s/m 2 t = 0. The only non vanishing integrals in this limit are the basis integrals f a 1 and f a 5 f a 1 = e 2γ E ε ε 2 Γ(ε) 2 and f a 5 = e 2γε ε 3 .
We checked our results for the MIs numerically in every kinematic regime against the evaluation with the program FIESTA 4.1. The numeric evaluation of the HPL's is performed using the GiNaC [64,65]. We have complete agreement within the numerical uncertainties of  . As a basis of integrals we make the following choice: Family j: C 2 contribution Family k: C 1 contribution which we compute with the help of HyperInt. We provide the Laurent-coefficients up to weight six in the ancillary files.

Two-Loop master integrals in H → bb
The complete set of scalar integrals of the EFT process H → bb can be parametrized by the three auxiliary families l, m and n defined in tab. 3. As a set of MIs we take the 25 canonical integrals defined in eq. (5.11), eq. (5.12) and eq. (5.13). The topologies corresponding to the MIs are shown in fig. 6 and fig. 7, where the dashed line corresponds to the Higgs, wavy lines denote massless and continuous straight lines massive propagators and external lines of mass  In order to compute the MIs we use the method of differential equation as described in paragraph 5.2. The large mass expansion of the MIs (p 1 · p 2 ) m 2 b corresponding to the expansion x ≈ 1 (see eq. (4.3)) is not as straightforward as for gg → H. We therefore compute only the small subset f l 1 = e 2γ E ε ε 2 Γ(ε) 2 (5.14) x (x−1) 2 2ε Γ(1 − ε)Γ(−ε) 2 Γ(2ε + 1) to all orders. Furthermore we need f l We have derived canonical bases for the integral families relevant for both calculations, which will allow the systematic calculation of higher orders in ε should they be required. We have checked the results obtained for the integrals numerically against sector-decomposition programs and we have compared the pieces of our amplitude against existing results when available.
Although they do not consitute physical observables in themselves, our results can be combined with other components to improve the predictions on the production and decay rates of the Higgs boson and we hope to combine them to future results to further our understanding of its properties.