Exact top-quark mass dependence in hadronic Higgs production

The impact of the finite top-quark mass on the inclusive Higgs production cross section at higher perturbative orders has been an open question for almost three decades. In this paper, we report on the computation of this effect at NNLO QCD. For the purely gluonic channel, it amounts to $+0.62\%$ relative to the result obtained in the HEFT approximation. The formally sub-leading partonic channels over-compensate this shift, leading to an overall effect of $-0.32\%$ at a $pp$ collider energy of $13\,\mathrm{TeV}$, and $-0.16\%$ at $8\,\mathrm{TeV}$. This result eliminates one of the main theoretical uncertainties to inclusive Higgs production cross section at the LHC.


INTRODUCTION
Gluon fusion is the dominant production process for a Standard Model (SM) Higgs boson at the Large Hadron Collider (LHC). Obviously, this makes a firm understanding of the pertinent cross section mandatory for current and future precision Higgs physics.
A comprehensive analysis of the theoretical prediction for the inclusive gluon fusion cross section has been performed in Ref. [1]. It arrives at an overall theoretical uncertainty of about ±5%, arising from six different sources, each of which contributes roughly 1% [49]. Thus, a significant reduction of the theoretical uncertainty cannot be achieved by eliminating a single source, but entails efforts on several of them.
Indeed, two of these sources have recently been addressed. The first one was due to the fact that the nextto-next-to-next-to-leading order (N 3 LO) QCD corrections were based on their expansion around z ≡ 1 − M 2 H /ŝ = 0 [2], giving rise to an uncertainty from the truncation of this expansion at finite order in z. Meanwhile, however, the exact dependence on z has become available [3]. The second one originated from the use of a factorization formula for the mixed QCD/electro-weak effects [4,5]. Recent progress indicates that this issue will be settled in the near future [6][7][8][9]. On the other hand, two other sources of uncertainty, namely the missing higher-order terms in both the partonic cross section and the parton density functions (PDFs), will require further technological advances before one can expect significant improvements.
The two remaining sources of uncertainty identified in Ref. [1] are related to quark mass effects. Both of them originate from the fact that the gluon-fusion process is induced by quark loops. The N n LO QCD corrections therefore involve an (n + 1)-loop calculation with at least two external mass scales (the Higgs and the quark mass, and possibly other quark masses from additional closed loops). While the next-to-leading order (NLO) result for arbitrary quark masses has been available for almost three decades [10], radiative corrections beyond this order were mostly restricted to top-loop induced terms [11][12][13], which make up around 95% of the total cross section in the SM. Their exact evaluation is the topic of the current paper. Progress in approximating bottom-and charm-quark effects beyond NLO has been made recently in Refs. [14][15][16][17][18].
It turns out that the dominant effect of the top-loop induced terms can be accounted for in the so-called Higgs effective field theory (HEFT) approximation, which is defined by multiplying the leading order (LO) cross section by the higher order (HO) correction factor in the limit of an infinite top-quark mass, M t → ∞ [19,20], which we take to be defined in the on-shell scheme throughout this paper: In this limit, the top-quark loop assumes the form of an effective Higgs-gluon vertex [21], thus reducing the number of associated loop integrations by one. At NLO, Eq. (1) approximates the full hadronic cross section for a SM Higgs boson to about 0.1%. This is remarkable for several reasons. On the one hand, the assumption that M t is the largest dimensional scale of the process is invalid over a large range of the partonic center-ofmass energy √ŝ , which reaches up to the collider energy √ s ≫ M t . On the other hand, less than 50% of the total cross section is due to the LO contribution, which means that the M t → ∞ approximation is applied to more than half of the total cross section. And finally, for the non-gg partonic channels such as qg or qq the HEFT approximation largely fails to capture the top-mass effects at NLO. It is only due to the strong numerical dominance of the gg channel that this hardly affects the total hadronic cross section.
Qualitatively, the high accuracy of the HEFT result can be explained by the suppression of the large-ŝ region by the PDFs. Also a dominance of the soft region in the total cross section could be made responsible for the small impact of the top-mass effects at higher orders. However, a solid quantitative understanding of this observation is still missing. The main reason for this is that higher order terms in 1/M t introduce positive powers ofŝ and thus enhance the large-ŝ region. Consequently, such terms cannot serve as uncertainty estimates of the heavy-top limit in a straightforward way.
So far the only estimate of top-mass effects beyond the HEFT approximation is therefore based on a combination of the 1/M t -expansion with the leading terms in the large-ŝ limit [22][23][24][25], from which an uncertainty of 1% due to top-quark mass effects was derived [1].
In this paper we eliminate this uncertainty by reporting on the exact calculation of the top-quark mass effects in hadronic Higgs production at next-to-next-to-leading order (NNLO) QCD.

CALCULATION
The calculation requires the combination of the purely virtual three-loop corrections to the cross section with the contributions from the real emission both of a single parton (quark q, anti-quarkq, or gluon g) at two-loop level, and of two partons at one-loop level. Factorization scheme dependence demands to take all possible partonic initial states into account. This is also important in the light of the failure of Eq. (1) for the non-gg channels as mentioned before, combined with the fact that they increase by roughly 100% from NLO to NNLO within HEFT.
Calculations of all the relevant amplitudes, including their full top-mass dependence, have already been reported on in the literature. In fact, the double-real emission amplitudes have been known for two decades [26]. Today, they can be obtained with public automated tools, and we use OpenLoops [27] for this purpose.
Complete results for the three-loop virtual amplitude are very recent. Its full top-mass dependence at NNLO has been first obtained with the help of Padé approximants constructed from the heavy-top expansion and the non-analytic terms at the thresholdŝ = 4M 2 t [28]. For the present study, we use a numerical result that was derived by expressing the amplitude in terms of master integrals, and subsequently evaluating them numerically [29]. Note that a fully analytic result is only available for the part which involves light (massless) fermion loops [30].
The main obstacle when calculating the total cross section with full top-mass dependence are the two-loop single-emission amplitudes. Unfortunately, existing results are not suitable for our purpose. After a number of approximate results [31,32], the amplitudes have been evaluated including their full top-mass dependence in the context of the Higgs-plus-jet production [33,34], but their numerical accuracy is insufficient, in particular since we need them also in the soft and collinear regions. Semianalytic results for the master integrals in the form of one-dimensional generalised power series have been presented as well [35]. In lack of a public code for these results, it would be necessary to implement the algorithm of Ref. [35] from scratch, which is a very demanding task.
In order to arrive at the required numerical precision, we have calculated the single-emission contribution by following the strategy of Ref. [29], which itself is based on Ref. [36]. In short, the amplitudes have been reduced to a set of master integrals with the help of the public software Kira⊕FireFly [37][38][39][40][41]. Algebraic manipulations have been simplified by setting the ratio of the top-quark and Higgs-boson mass to a fixed value of M 2 t /M 2 H = 23/12, corresponding to M t ≈ 173.055 GeV for M H = 125 GeV. The same software has also been used to derive a system of first-order homogeneous linear differential equations in M t satisfied by the master integrals. Using initial conditions in the heavy-top limit, obtained with a diagrammatic large mass expansion, the system of differential equations has been solved numerically at a very large number of phase-space points. As a result, the amplitudes have been obtained on a dense grid that could, in principle, be used for interpolation. However, since the grid does not extend to the boundaries of the phase-space where the amplitudes diverge, direct inclusive phase-space integration requires a non-trivial extrapolation to the singular soft and collinear regions.
In order to evaluate the phase-space integrals, the amplitudes have been regulated in the soft and collinear limits by subtracting their counterparts in HEFT. For instance, for the gg → gH process assuming |t| < |û|, witĥ t,û the partonic Mandelstam variables, the contraction of the two-loop, |M (2) ⟩, and one-loop, |M (1) ⟩, amplitudes, both treated as vectors in color and spin space, has been replaced with: with |F 1,2 ⟩ the respective amplitudes for the gg → H process, and ⟨P  (1) exact |M (2) exact ⟩ regulated , defined in Eq. (2), for the processes gg → gH (first row), qg → qH (second row) and qq → gH (third row), separated into the region below (left column) and above (right column) threshold for intermediate top-quark pair production,ŝ = 4M 2 t . A factor of α 4 s /(4π) 2 · 1/v 2 ·ŝ, with v the Higgs-field vacuum expectation value, has been factored out. The renormalization scale has been set at µR = MH/2. The kinematics is parameterized with z ≡ 1 − M 2 H /ŝ and λ ≡t/(t +û).ŝ,t andû are the standard Mandelstam variables.
It is not difficult to convince oneself that ⟨M (1) exact |M (2) exact ⟩ regulated is devoid of soft and collinear singularities. We use an analogous expression for the qg → qH process with: instead of ⟨P (0) gg (z)⟩, while there is no splitting-function contribution in the qq → gH case. The regulated amplitudes for the three processes are illustrated in Fig. 1. Notice that the amplitudes for the first two processes are still singular in the collinear limit, but these singularities are integrable and occur very close to the edge of the phase space. In order to obtain a reliable inclusive phase-space integral, we approximate the regulated amplitudes in the limitt → 0 with the ansatz a ln |t| + b, for each value of z, and fit the coefficients to the available numerical values of the amplitudes at |t| > |t 0 |. The ansatz is subsequently integrated analytically in the region 0 > |t| > |t 0 |. The uncertainty of the procedure is estimated by using the more general ansatz, a ln |t| + b + ct ln |t| + dt.
Having the amplitudes for all contributions at hand, it is necessary to actually integrate them in order to obtain the cross section contributions. Since the effect of the top-quark mass beyond the heavy-top limit is expected to be small, we directly evaluate the difference of the cross sections at each phase-space point: rather than the cross sections themselves separately. This has the additional advantage that ultraviolet and infrared divergences in the form of 1/ϵ k poles in the dimensional regularization parameter, ϵ, as well as soft and collinear singularities first appear at the NNLO level. Hence, for example, the NLO contributions to the difference are welldefined separately for the virtual and real corrections. This delay of the appearance of divergences and singularities is one of the reasons for the smallness of the topquark mass effects beyond HEFT. Ultimately, Eq. (5) is evaluated with Monte Carlo methods using the sector-improved residue subtraction scheme [42][43][44] implementation in the C++ library Stripper. Note that it suffices to use the subtraction term in the square brackets of Eq. (2) in order to cancel the IR divergences with the double real emission. Since this subtraction term is given in terms of compact analytic formulae [45], it allows for a fast and numerically stable Monte Carlo integration. The phase space integration and PDF convolution of the r.h.s. of Eq. (2) is done separately. Adding it to the output from Stripper cancels the subtraction term contribution and leads to the final result.

RESULTS
Table I collects our main results. It shows the hadronic cross section σ NNLO HEFT in the HEFT approximation through NNLO QCD, including only top-loop induced contributions and without electro-weak effects, and separately for the partonic sub-channels (qq denotes the sum over all quark initial states). The absolute numbers are split into the contributions from the individual orders in α s . The uncertainties indicate the Monte Carlo integration errors.
While the finite-mass effects are small and positive for the gg channel (and largely independent of the collider energy), the relative effect on the other channels is negative and much larger. For the pure quark channels, the HEFT approximation is off by more than 100% at each perturbative order. Taken individually, this would already exhaust the uncertainty estimate associated with the missing mass effects of Ref. [1], despite the fact that these channels contribute to the total cross section only Table I: Effects of a finite top-quark mass on the total hadronic Higgs-boson production cross section for the LHC @ 13 TeV and 8 TeV, separately for the partonic channels and including Monte Carlo integration error estimates. Results obtained with the PDF set NNPDF31 nnlo as 0118 [46], renormalization and factorization scales µR = µF = MH/2, Higgs-boson mass MH = 125 GeV, and top-quark mass Mt = 23/12 × MH ≈ 173.055 GeV. The NNLO cross section within HEFT (σ NNLO HEFT ) has been obtained with SusHi [47,48] and is split into contributions from the individual orders in αs.

CONCLUSIONS AND OUTLOOK
A calculation of the hadronic Higgs production cross section including the full top-mass dependence at NNLO was reported. It results in a slight decrease relative to the result in the HEFT approximation of −0.32% at 13 TeV, and −0.16% at 8 TeV collider energy. This result confirms and at the same time eliminates the commonly accepted uncertainty estimate arising from the lack of knowledge of these effects.
Our calculational techniques are also applicable to the bottom-and charm-loop induced terms and the associated interference with the top-loop terms. This is deferred to future work. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grants 396021762 -TRR 257 and 400140256 -GRK 2497: The physics of the heaviest particles at the Large Hardon Collider.
Simulations were performed with computing resources granted by RWTH Aachen University under projects rwth0414 and rwth0643.