Exact Top Yukawa corrections to Higgs boson decay into bottom quarks

In this letter we present the results of the exact computation of contributions to the Higgs boson decay into bottom quarks that are proportional to the top Yukawa coupling. Our computation demonstrates that approximate results already available in the literature turn out to be particularly accurate for the three physical mass values of the Higgs boson, the bottom and top quarks. Furthermore, contrary to expectations, the impact of these corrections on differential distributions relevant for the searches of the Higgs boson decaying into bottom quarks at the Large Hadron Collider is rather small.


INTRODUCTION
The discovery of the Higgs boson [1,2] by the AT-LAS [3] and CMS [4] experiments at CERN has ushered in a new era in particle physics phenomenology. The Standard Model (SM) of elementary particles is now complete and there is no decay or scattering phenomenon at low energies that significantly deviates from what is predicted by the SM. Still, we know that the SM cannot be the ultimate theory, if not for the lack of consistency with mainly cosmological observations, but for the fact that it contains quite a large number of parameters, which makes it unreasonable to think of it as a fundamental theory. The LHC is guiding the experimental community towards the study of the Higgs potential and the Higgs boson direct couplings with all the other particles, an experimentally previously completely unexplored sector of the SM Lagrangian. The SM makes very precise predictions for all the vertices and couplings of the Higgs boson and the verification of these predictions is among the fundamental questions addressed at the LHC.
In particular, the gluon fusion production mechanism has given direct access to Higgs boson decay into vector bosons and an indirect access the top Yukawa coupling. More recently, also the direct coupling of the Higgs boson to the top quark has been observed [5,6]. Measurements of the Higgs coupling to the tau lepton have been extracted by combining all production modes [7,8], while the direct coupling to the bottom quark has been observed by exploiting the features of the VH (V = W or Z) associated production mechanism [9,10]. The decay to bottom quarks is quite special, because it is the one with the largest branching ratio. The decay width of the Higgs boson into bottom quarks has been computed at up to four loops in QCD [11][12][13][14][15][16][17][18] using an approximated treatment of the bottom quark mass, up to one loop including electroweak corrections [19,20] and also including mixed QCD-electroweak effects [21]. The exact bottom mass corrections have been computed up to two loops [22]. A relatively large component of the two loop computation is represented by the diagrams in which the Higgs boson couples to a top quark loop, see Figs. 1 and 2. These two sets of diagrams are both UV and IR finite separately, and their contributions have been computed only approximately in [16], finding a very compact formula that should be considered valid for values of the masses such that m b m H m t . Comparing this formula with the rest of the two loop contributions, it turns out that these pieces, proportional to the top Yukawa coupling y t , account for about 30% of the total two loop result.
The aim of this letter is twofold. First, we want to assess the impact of the neglected terms in the expansion of [16], that in principle could be of the order of (m H /m t ) 4 ∼ 20% (see Eq. (3) in [16]  quarks up to two loops for the massless case, and merged this computation to the two loop corrections to the associated production in hadronic collisions [23,24]. The corrections to key distributions like the transverse momentum and the mass spectra of the Higgs boson (reconstructed using the two hardest b-jets in the final state) are found to be very large. y t contributions to the Higgs decay into bottom quarks are not included in the differential results mentioned above and so it is natual to ask what the impact of these corrections is (see for example [24]). We answer this second question by presenting differential results which include the y t contributions to the decay and retain the full mass dependence on the top and bottom quark masses.

Double virtual
The decay H(q) → b(p 1 ) +b(p 2 ) receives O(α 2 s y t ) contributions from the interference |M yt,bb | 2 of the Born amplitude with two loop virtual corrections that involve a closed top quark loop, where M (2) yt,bb is given by the two diagrams shown in Fig. 1. By evaluating the Feynman diagrams, we decomposed |M yt,bb | 2 as with a j ∈ Z and i = t, b, H. In Eq. (2), the c a1···a7 are rational coefficients and I a1···a7 are two loop integrals of the type defined by the set of inverse propagators: We computed the loop integrals through the consolidated differential equations (DEs) method [25][26][27][28]. First, we used integration-by-parts identities (IBPs) [29][30][31], generated with the help of Reduze2 [32], in order to reduce the integrals that appear in |M yt,bb | 2 to a set of 20 independent master integrals (MIs) I = (I I , . . . , I 20 ). Subsequently, we derived the analytic expression of the MIs by solving the system of coupled first-order DEs in the kinematic ratios m 2 H /m 2 t and m 2 b /m 2 t . The structure of the DEs, and hence of their solutions, is simplified by parametrizing such ratios in terms of the variables x and y, defined by and by using the Magnus exponential method [33][34][35][36][37][38] in order to identify a basis of MIs that fulfil a system of canonical DEs [39], In Eq. (6), the coefficient matrix dA is a dlog-form that contains 12 distinct letters, with M i ∈ M 20×20 (Q). Since all letters are algebraicallyrooted polynomials, we could derive the -expansion of the general solution of Eq. (6) in terms of twodimensional generalized polylogarithms (GPLs) [40][41][42][43][44], by iterative integration of the dlog-form, which we performed up O( 4 ), i.e. to GPLs of weight four. In order to fully specify the analytic expression of the MIs, we complemented the general solution of DEs with a suitable set of boundary conditions. The latter were obtained by demanding the regularity of the MIs at the pseudothresholds m 2 H = 0 and m 2 H = 4m 2 b that appear as unphysical singularities of the DEs.
The expression of the MIs obtained in this way is valid in the Euclidean region m 2 where the logarithms of Eq.  Upon inserting the expressions of the MIs into Eq. (2), we observed the expected analytic cancellation of allpoles and obtained a finite result for |M yt,bb | 2 , with C(x, y) being a polynomial combination, with algebraic coefficients, of 256 distinct GPLs. The explicit expression of C(x, y), as well as of the newly computed MIs, can be released by the authors upon request.

Real-virtual
The real-virtual part of the computation involves the interference |M yt,bbg | 2 of the tree level amplitude for H(q) → b(p 1 ) +b(p 2 ) + g(p 3 ) with the loop diagrams in Fig. 2 containing a closed top quark loop, We used standard techniques to evaluate the one loop amplitude. As expected, |M yt,bbg | 2 is finite (in ) and can be written as plus terms that vanish in four dimensions. In Eq. (10), s ij denotes twice the dot-product of momenta, s ij ≡ 2p i · p j . We integrated the real-virtual contribution over the whole phase space both analytically and numerically using Monte Carlo integration, finding perfect agreement. The analytic computation was performed by direct integration of |M yt,bbg | 2 over the three-particle phase space. The phase space measure for the decay where ∆ 3 is given by The integral is finite in four dimensions and was evaluated in terms of GPLs after suitable transformations of the integration variables. In particular, square roots involving the integration variables appear at intermediate stages of the calculation (both from the one loop matrix element and from resolving the phase space constraint implied by the positivity of ∆ 3 ) and must be linearized, e.g. by using the techniques of [47]. The full result is represented in terms of a formula with 1841 distinct GPLs and can be released by the authors upon request. The numerical integration of the real-virtual contribution is straightforward and has been used to validate the analytic computation. It also allows to build Monte Carlo simulations with acceptance cuts and has been used to obtain the differential result of the next section.

RESULTS
We begin the presentation of our results by discussing the inclusive decay rate. In Table I, we compare our exact formula, obtained from the sum of the double virtual and real-virtual contributions described in the previous section, to the approximated one of Ref. [16]. The numbers in the table are obtained with the following formula for the relative discrepancy among exact and approximated 4 results: The agreement is excellent for the physical mass values, proving for the first time and in a completely independent way the validity of the approximated formula, and the fact that it works much better then expected. We now turn to the second question regarding the impact of the y t contribution at differential level. To this extent, we present results for Higgs boson associated production and to avoid the contamination from initial state radiation we consider pp → W + (l + ν l )H(bb) at leading order and add the corrections to the decay process at the next-to-leading order. Then, we compare this result with the one obtained by adding also the y t contribution. Note that in both cases we normalize the cross section to the total Higgs boson decay width into bottom quarks reported in the Yellow Report of the Higgs Cross Section Working Group [48] (HXSWG), that includes higher order corrections. So, effectively, we are comparing the shapes of distributions. To obtain our results, we use the SM parameters recommended by the HXSWG and the NNPDF3.0 [49] parton distribution functions. Furthermore, we impose the following typical lepton acceptance cuts: selected events must have a missing transverse momentum greater than 30 GeV, the charged lepton is required to have a transverse momentum greater than 15 GeV and an absolute rapidity smaller than 2.5 and, finally, the W boson transverse momentum is required to be larger than 150 GeV. We reconstruct jets using the anti-kt algorithm with the resolution variable set to 0.5 and require two b-jets with transverse momentum greater than 25 GeV and absolute rapidity smaller than 2.5. the impact of the y t contribution on both the transverse momentum distribution and the mass distribution of the putative Higgs boson is extremely small with at most a 5% effect in the low energy tail of the transverse momen-5 tum distribution. These corrections are much smaller than the scale variation uncertainty of the computation.

CONCLUSIONS
In this letter we presented the results of the full analytic computation of the top Yukawa contribution to the Higgs boson decay width into bottom quarks. First, we demonstrated that the approximate formula used so far in the literature works exceedingly well for physical values of the masses. This nice behaviour was not predictable a priori and, with respect to a possible estimate of about a 20% error, we have instead found smaller than per mill deviations of this formula from the exact result. Then, we showed that the impact of this contribution at the differential level is very small and Monte Carlo simulations performed so far are not affected by an additional significant source of uncertainty due to the ne-