Analytic decay width of the Higgs boson to massive bottom quarks at next-to-next-to-leading order in QCD

The Higgs boson decay to a massive bottom quark pair provides the dominant contribution to the Higgs boson width. We present an exact result for such a decay induced by the bottom quark Yukawa coupling with next-to-next-to-leading order (NNLO) QCD corrections. We have adopted the canonical differential equations in the calculation and obtained the result in terms of multiple polylogarithms. We also compute the contribution from the decay to four bottom quarks which consists of complete elliptic integrals or their one-fold integrals. The result in the small bottom quark mass limit coincides with the previous calculation using the large momentum expansion. The threshold expansion exhibits power divergent terms in the bottom quark velocity, which has a structure different from that in $e^+e^-\to t\bar{t}$ but can be reproduced by computing the corresponding Coulomb Green function. The NNLO corrections significantly reduce the uncertainties from both the renormalization scale and the renormalization scheme of the bottom quark Yukawa coupling. Our result can be applied to a heavy scalar decay to a top quark pair.


Introduction
The Higgs boson was discovered more than ten years ago, and its couplings with other elementary particles in the standard model (SM) have been measured and found to be consistent with the SM predictions [1][2][3].The Yukawa couplings of the Higgs boson, along with the spontaneous electroweak symmetry breaking, play a crucial role in explaining the origin of the masses of fermions.The constraints on the Yukawa couplings with the top quark [4,5], the bottom quark [6,7], the charm quark [8,9], the tau lepton [10,11], and the muon [12,13] have been set by the ATLAS and CMS collaborations.Because the top quark is very heavy, the Higgs boson can not decay into top quarks.The dominant decay mode of the Higgs boson is thus H → b b.From the signal strength observed at the LHC [6,7], the bottom quark Yukawa coupling is derived with an error of ±10%.The accuracy will be significantly improved at future lepton colliders such as the Circular Electron Positron Collider (CEPC) [14][15][16], and the precision can reach the per mill level [15,16].
The corresponding theoretical predictions have been improved continuously in the past half-century.The next-to-leading order (NLO) QCD corrections [17][18][19][20] were computed more than 40 years ago.And the NLO electroweak (EW) corrections have also been known for a long time [21,22].The mixed QCD and EW corrections of order O(αα s ) have been obtained in Refs.[23,24].Later, the decay width of H → b b was calculated up to O(y 2 b α 4 s ) [25][26][27][28][29], when the final-state bottom quark is taken to be massless but the bottom Yukawa coupling is still finite.The corresponding differential results were investigated up to O(α 3  s ) [30][31][32].The inclusive and differential decay widths with finite m b were computed at NNLO numerically [33][34][35].The effects of parton shower in this decay process have been explored in [36,37].The differential decay rates of H → b bj at NNLO and H → b bb b at NLO were presented in [38] and [39], respectively.
The decay H → b b can also be induced by a top quark loop with the top quark Yukawa coupling y t , which is much larger than y b .However, this contribution starts at the two-loop level, i.e., O(y b y t α 2 s ).Because the helicity should be flipped along the bottom quark current, an additional bottom quark mass suppression emerges.As a result, the contribution is not significantly larger, or even smaller, than that of O(y 2 b α 2 s ).The analytic result of the decay width of H → b b at O(y b y t α 2 s ) with finite m b and m t has been calculated in [40].The calculation of the differential decay width that is sensitive to y t was presented up to O(α 3 s ) in [29,41].The decay rate of H → hadrons including contributions from both H → b b and H → gg has been computed up to O(α 4 s ) [42,43].In this paper, we provide an exact result for the decay H → b b at O(y 2 b α 2 s ) with full dependence on m b .Specifically, we consider the contributions from the decay processes H → b b(+X), X = g, gg, q q, b b, where q stands for the massless quark.We note that an analytical calculation has been performed in [44,45], but only the expansion result in m 2 b /m 2 H was obtained.Our analytic results enable an efficient computation of the NNLO decay width, and can be readily extended to other similar processes, such as a new heavy scalar decaying to top quarks.Moreover, given the analytic result, we can explore different aspects of the decay rate.The expansion of the decay width near the threshold shows power divergences, which violate the perturbative expansion in α s .To understand the origin of these divergences and to restore the power of perturbative expansion, we need to calculate a new Coulomb Green function.
This paper is organized as follows.Our calculation framework is introduced in section 2. The master integrals in the calculation of the decay width are computed in section 3.In section 4, the renormalization procedure is briefly described.The analytical and numerical results for the decay width are presented in section 5 and section 6, respectively.Our conclusion is given in section 7.

Calculation framework
We are going to calculate the decay width of H → b b, denoted by Γ Hb b, via the optical theorem, where Σ represents the amplitude of the process H → b b → H.In this method, we do not need to calculate the virtual and real corrections separately, and the complicated multi-body phase space integration can be avoided.
We (2.2) At O(y 2 b α 2 s ), we have to compute the imaginary part of three-loop self-energy diagrams, as shown in figure 1.We generate the diagrams with the package FeynArts [46], and perform the Dirac algebra using the package FeynCalc [47].The amplitudes are expressed as linear combinations of three-loop scalar integrals, which are reduced to a minimal set of integrals called master integrals (MIs) due to the identities from integration by parts (IBP) [48,49].The solution of these identities is found by the Laporta algorithm [50] that is implemented in the packages FIRE [51] and Kira [52].These MIs can be represented by two integral families, i.e., NP1 and P1 in figure 2. In this procedure the package CalcLoop1 has been used to shift the loop momenta.In our calculation, the Higgs boson and the bottom quark are taken to be massive while the light quarks are massless.We do not consider the top quark contribution.According to the cutting rules in the optical theorem, the imaginary parts of Σ come from the cut diagrams in which some propagators are put on-shell simultaneously.We find that most cut diagrams in Σ have at least two on-shell bottom-quark propagators.And   Notice that some cut diagrams contribute to the imaginary part of Σ, but do not have two bottom-quark cut propagators, such as the diagram of a cut through two gluon propagators as depicted in figure 5.This cut diagram is related to the process H → gg, which should not be included in the corrections of H → b b.To subtract this contribution from the imaginary part of the forward scattering amplitude, we calculate the squared amplitude of H → gg with one triangle bottom-quark loop, and perform the two-body phase space integration.

Calculation of master integrals
As mentioned in eq.(2.2), we have to calculate the contributions from both two and four bottom quark final states.The corresponding MIs can be computed separately.
The MIs contributing to the decay to four bottom quark final state contain elliptic integrals, such as those shown in figure 6.They have been calculated in [53] where the authors consider the total cross section of e + e − → Q QQ Q.After choosing a regular basis, only the O(ϵ 0 ) parts of the MIs are needed, which in turn can be expressed either as complete elliptic integrals of the first kind or one-fold integrals of them.Following this method, we present the decay width Γ y b y b Hb bb b in terms of one-fold integrals over expressions depending on complete elliptic integrals and multiple polylogarithms.The explicit form can be found in eq.(5.8) below.
Then we turn to the MIs that are needed in the calculation of Γy b y b Hb b .We adopt the method of differential equations [54,55].Once a set of canonical basis F is found, the dimensional regulator ϵ is factorized out from the kinematic variables [56], where {x n } denote the kinematic variables.When dA(x n ) is expressed in d log forms, the differential equations can be solved recursively in terms of multiple polylogarithms (MPLs) [57], which are defined by G(x) ≡ 1 and 3) The number of elements in the set {l 1 , l 2 , . . ., l n } is referred to as the transcendental weight of the MPLs.Note that we have dropped MIs that possess an imaginary part only from cuts on four bottom quarks.And all the boundary values are evaluated in the phase space where it is not possible to apply a cut on four bottom quarks.For instance, the left diagram in figure 6 is dropped in the basis.The right diagram in figure 6 is included but only the cut through two bottom quarks, as shown in the right diagram of figure 7, is considered.

Canonical basis of NP1
The master integrals of the NP1 integral family in figure 2 are defined by with where all n i ≥ 0, i = 1, • • • , 8, and n 9 ≤ 0. The denominators D i read where the Feynman prescription +iε has been suppressed.The external momentum k satisfies k 2 = m 2 H .We focus only on the imaginary part of the integral, which gets contributions from all possible cuts except for the ones on four bottom quarks.The results of the master integrals for specific cuts have been obtained in [58].
In the NP1 integral family, there are 29 MIs contributing to ΓHb b.To construct the canonical basis, we first select the MIs Here all M NP1 i are dimensionless, and the corresponding topology diagrams are displayed in figure 11 in appendix A.
Then the canonical basis integrals F NP1 i , i = 1, . . ., 29, can be constructed as linear combinations of M NP1 i using a method similar to that in ref. [59].Examples can be found in [60].For simplicity, we define a dimensionless variable z ≡ m 2 H /m 2 b + i0 + .Then we obtain the following canonical basis of the NP1 family: where r = z(z − 4). (3.9) To rationalize the square root r, we write The corresponding differential equations for the canonical basis F NP1 become and N NP1 i being constant rational matrices.

Canonical basis of P1
The P1 integral family is defined by where the denominators D i read The topology diagrams are displayed in figure 12 in appendix A. We find the corresponding canonical basis integrals, (3.17) The corresponding differential equations are and N P1 i being constant rational matrices.

Boundary conditions
Given the above canonical differential equations, we obtain the results of the basis integrals at each order of ϵ up to some integration constants.The MPLs in the results can be evaluated efficiently using the PolyLogTools package [61], which depends on the algorithm of ref. [62] in the GiNaC [63] framework.The integration constants are calculated numerically by comparing the results of basis integrals with those obtained using the AMFlow package [64][65][66][67].Their analytical form is reconstructed using the PSLQ algorithm [68,69].We find that the transcendental constants relevant to our calculation consist of up to weight four, where ζ(n) is the Riemann zeta function.All the master integrals are then checked against the numerical computation with AMFlow, and over 200-digit agreement is achieved.We provide the analytical results of the canonical bases in an ancillary file along with the arXiv submission.

Renormalization
In the framework of the optical theorem, the imaginary part of each three-loop self energy diagram arises from all possible cuts, corresponding to the sum of both virtual and real corrections.The infrared divergences cancel out in the sum, and only ultraviolet divergences remain.
We perform renormalization following the standard procedure, i.e., by including the contribution from diagrams with counter-terms.The bottom quark mass is renormalized in the on-shell scheme, while the strong coupling α s is renormalized in the MS scheme, where C ϵ = (4πe −γ E ) ϵ .The renormalization constants up to α 2 s are given by [70][71][72][73], where and color factors C F = 4/3, C A = 3, T R = 1/2 have been substituted.The bottom quark Yukawa coupling is related to the bottom quark mass.However, one has the freedom to choose a different renormalization scheme.Though it is expected that the final result for the decay width should be the same if all order corrections are included, the fixed-order results in different schemes may exhibit obvious deviations, which should be taken as one kind of theoretical uncertainty.The results with Yukawa couplings in the on-shell and MS scheme can be converted to each other by applying the relations between the two schemes, which have been computed up to O(α 4 s ) [74][75][76][77].We have performed two separate calculations using different renormalization schemes.The ultraviolet divergences cancel in either scheme, and the final results agree with each other after the conversion of the renormalization scheme.

Full results
The decay width of H → b b up to O(α 2 s ) can be expressed as ) . ( The LO result is where N c = 3 is the number of colors and y b is the MS Yukawa coupling with m b the MS mass of the bottom quark and .The latter denotes the contribution from the interference of the bottom and top quark Yukawa coupling induced amplitudes and starts at O(α 2 s ).The corresponding analytical result has been calculated in ref. [40].The former receives contributions from both two and four bottom quark final states, as discussed in section 2, + 6G(0, x 2 , 0, −1, w) − 4G(0, x 2 , 0, 0, w) + 3G(0, x 2 , 0, 1, w) − G(0, x 2 , 1, 0, w) − 2G(0, 0, w) + 2iπG(0, w)   with where β = 1 − 4/z measures the velocity of the bottom quark.The master integrals F 4b 1 (z), F 4b 2 (z) and F 4b 3 (z) are linear combinations of complete elliptic integrals or their derivatives, and the other F 4b i (z) can be expressed by one-fold integrals of MPLs and complete elliptic integrals.The complete expressions of all F 4b i (z) can be found in [53].We also provide these expressions in electronic form in the ancillary file attached to this paper.

Asymptotic expansion
It is interesting to explore the decay width in various limits.First, we study the small bottom quark mass limit, i.e., m 2 b → 0 or z → ∞, while the Yukawa coupling is still kept finite.The NLO coefficient becomes Therefore, it is safe to take the continuous limit of m b = 0.This property is in accordance with the Kinoshita-Lee-Nauenberg theorem [78,79].Here the parameter m b can be considered as a regulator for infrared divergences.Though the virtual and real corrections contain log(z) terms induced by collinear gluons separately, the dependence on the regulator cancels in the decay width at O(z 0 ).Notice that there are large log(z) terms at O(z 0 ) if the Yukawa coupling is renormalized in the on-shell scheme.
In the limit of z → ∞, the NNLO coefficient where we have shown explicitly the first two orders in the expansion.It is obvious that there are large logarithms log i (z) with i ≤ 3 at O(z 0 ).This means that the contribution of this part with two bottom quarks in the final state to the decay width is divergent in the massless limit. ( Figure 8: Typical cut diagrams that have large logarithms log i (z).The thick black and red lines stand for the massive bottom quark and the Higgs boson, respectively.And blue line represents the particle that is soft.
These logarithms log i (z) arise from the diagrams with a bottom quark loop, e.g., see the diagram (a) in figure 8.These diagrams contain up to 1/ϵ 3 infrared divergences in the massless case, which are converted to log i (z), i ≤ 3, in the massive bottom quark case.
It is interesting to note the log 4 (z) term at order O(z −1 ).It appears due to a different reason.It comes from the contribution of diagrams with two soft quarks exchanged between two separate directions, e.g., see the diagram (b) in figure 8.
It can be seen that X y b y b 2,b bb b also have large logarithms when z → ∞.The logarithms log i (z), i ≤ 3, at O(z 0 ) are obtained by calculating the diagrams with a soft gluon splitting to a soft bottom quark pair; see the diagram (c) in figure 8.In contrast, the logarithm log 4 (z) at O(z −1 ) arises from the diagrams with a collinear quark splitting to a soft quark and a collinear gluon which splits into a collinear quark and a soft quark further, such as the diagram (d) in figure 8.
Combining the results in eq.( 5.10) and eq.( 5.11), we obtain (5.12) The large logarithms of z at O(z 0 ) cancel, as expected.However, the large logarithms at O(z −1 ) survive.The structure of this kind of large logarithm at higher orders in α s is not well understood.It is promising to explore this issue using the effective field theory, which we leave for future study.After adopting the same renormalization scheme, the above expression is consistent with the result in ref. [45], which was obtained with the method of large momentum expansion.
It is also interesting to investigate the threshold limit of bottom quark pair production, where the coefficients X y b y b 1 and X y b y b 2 are expanded in the limit of m H → 2m b or β → 0. Expanding the full analytic results, we have In obtaining the above compact expressions, we have extensively used the relations among MPLs in [80,81].We check that the coefficients of 1/β2 and 1/β are consistent with the results of two-loop form factors [82], which means that the power divergences at β → 0 in X y b y b 1 and X y b y b 2 come only from virtual corrections.Actually, these velocity-enhanced terms are due to the Coulomb potential interaction between the bottom quark pair.They violate the convergence of perturbative expansion 2 and need to be resummed to all orders in α s by calculating the corresponding Coulomb Green function [83].Notice that the coefficients of π 2 /β in X y b y b 1 and π 4 /β 2 in X y b y b 2 agree with those in e + e − → t t production near threshold [84], but the π 2 /β 2 term is different.They can be reproduced by calculating the imaginary part of the first a few terms of the zero-distance Coulomb Green function, The all-order result of this Coulomb Green function can be found in refs.[85,86].
The LO decay width in the on-shell scheme does not depend on the scale.The scale uncertainty is about 11% and 13% at NLO and NNLO, respectively.At a typical scale µ = m H , the NLO correction decreases the LO result by 35%, and the NNLO correction reduces the decay width further by 23%.This behavior indicates that the perturbative expansion with the on-shell Yukawa coupling converges slower than that with the MS Yukawa coupling.
Comparing the two schemes, one can find that the results in the on-shell scheme are larger than those in the MS scheme.The main reason is that the pole mass m b = 5.07 GeV is much larger than the MS mass m b (µ = m H ) = 2.78 GeV.However, higher-order corrections reduce the difference between the two schemes prominently, signifying its importance in providing reliable theoretical predictions.
Given that the scheme dependence constitutes the dominant theoretical uncertainty now, it is necessary to see how it changes as higher orders are taken into account.We adopt the O(y b y t α 2 s ) correction in the expansion of m ) corrections in the massless bottom quark limit [26][27][28].One can see that the central value for the decay width in the MS scheme becomes stable quickly and that the scale uncertainty is decreased to ±0.3% after the N 4 LO correction is included.In contrast, the N 3 LO correction in the on-shell scheme decreases the NNLO result by 12%, and the N 4 LO correction reduces the decay width by 5.5% further.The scale uncertainty is still 2% at N 4 LO.All these features manifest that the results in the on-shell scheme have a bad convergence behavior.It is more appropriate to take the MS scheme in this high-energy process.
Our analytic results can also be applied to the decay of heavy scalar particles.If there is a new scalar Higgs with a mass of more than 350 GeV, it can decay to a top quark pair.In this process, we have five massless light quarks and one massive top quark.The decay width Γ Ht t of H → t t with m t = 172.69 GeV and different values of m H are shown in figure 10.The NNLO correction increases the NLO result by 3% − 17% at µ = m H as m H varies from 400 GeV to 1 TeV, and the scale uncertainty reduces from 7% − 13% at NLO to 2% − 8% at NNLO.

Conclusion
We present an exact result of the decay width of  diagrams of the Higgs boson which have two bottom quarks in the propagators.The master integrals are calculated using the method of differential equations combined with the choice of a canonical basis.The result has been expressed in terms of multiple polylogarithms.The contribution from the four bottom quarks is also calculated analytically and written as a linear combination of complete elliptic integrals and one-fold integrals of them.The expansion of our analytical result in the small quark mass m b limit agrees with the previous computation in a different method.The threshold expansion shows velocity enhancement which is caused by Coulomb potential interactions.Though the structure is not the same as that in the top quark pair production near the threshold, we checked that the power divergent terms can be reproduced by calculating a new Coulomb Green function.We discuss the numerical results in both the MS and on-shell renormalization schemes for the Yukawa coupling.Higher-order QCD corrections reduce both the scale uncertainty and the difference between the two schemes.The analytical results presented in this paper enable accurate and fast evaluation of the decay width of heavy scalar particles.
define the contribution from the final states of two bottom quarks as Γy b y b Hb b , and that from final states of four bottom quarks as Γ y b y b Hb bb b.The partial decay width of H → b b is thus given by Γ y b y b Hb b ≡ Γy b y b Hb b + Γ y b y b Hb bb b .

Figure 1 :Figure 2 :
Figure 1: Typical three-loop Feynman diagrams of H → b b → H.The thick black and red lines stand for the massive bottom quark and the Higgs boson, respectively.
these cut diagrams correspond to the processes H → b b, H → b bg, H → b bgg, H → b bq q and H → b bb b.All of them contribute to Γ Hbb at O(y 2 b α 2 s ).Two sample cut Feynman diagrams are shown in figure 3. The requirement of cuts on two bottom quarks simplifies the calculations since some MIs do not have imaginary parts, e.g., the three-loop massive vacuum bubble diagram shown in figure 4.

Figure 3 :
Figure 3: Typical cut diagrams.The black dashed line stands for the cut line, and the propagators crossing the cut line are on-shell.The cut line passes through two (left graph) and four (right graph) bottom quarks.

Figure 4 :
Figure 4: A three-loop massive vacuum bubble diagram.

Figure 5 :
Figure 5: A cut diagram contributing to the H → gg decay.

Figure 6 :
Figure 6: Two master integrals contributing to the decay to four bottom-quark final states.The thick black and red lines stand for the massive bottom quark and the Higgs boson, respectively.

Figure 7 :
Figure 7: The four (left) and two (right) bottom quark cut diagrams.The thick black and red lines stand for the massive bottom quark and the Higgs boson, respectively.

. 15 )
There are 19 MIs in this integral family contributing to ΓHb b.Most of them have appeared in the NP1 integral family discussed above.Here we only show the new MIs,

. 7 )
Here the color factors C F = 4/3, C A = 3, and T R = 1/2 have been substituted by their numeric values, and n l denotes the number of light quark flavors.The result of X y b y b The other NNLO coefficient X y b y b 2,b bb b consists of master integrals that depend on elliptic integrals.We present the asymptotic expansion results of the master integrals in appendix B. Based on these results, we obtain the asymptotic expansion, X y b y b 2,b bb b| z→∞ = log

= m b β 2 +
O(β 4 ).The p i • p j numerator arises because the scalar current ψb ψ b with ψ b the bottom quark field matches to the current j (s) = ψ † σ • iDχ/m b in non-relativistic QCD where ψ and χ represent the two-component quark and anti-quark field, respectively.

[ 3 ]X y b y b 1 (X y b y b 1 (X y b y b 2 (X y b y b 2 (Table 1 : 1 ,b y b 2 and X y b y b 2 coefficients. The numerical results of X y b y b 1 and X y b y b 2 at 2 .
m b (m b ) = 4.18 GeV, m b = 5.07 GeV, m H = 125.09GeV, m Z = 91.1876GeV, (6.1) α s (m Z ) = 0.1181, G F = 1.166378 × 10 −5 GeV −2 , where the on-shell value of m b is obtained from m b (m b ) in the MS scheme with the fourloop mass relation.The package RunDec [87, 88] was used to obtain m b at other scales, e.g., m b (m H /2) = 2.956 GeV, m b (m H ) = 2.784 GeV and m b (2m H ) = 2.639 GeV.our res.)+3.037513882 +5.810102605 +8.582691327 Ref. [45]) +3.037513882 +5.810102605 +8.582691327 Xy b y b 2,b b (our res.) −4.537634223 +29.22467780 +78.04118427 our res.)−3.181291730 +30.58102030 +79.39752676Ref. [45]) −3.181291730 +30.58102030 +79.39752676The numerical results for the X y b y b Xy the scale µ = m H /2, m H , 2m H are shown in table 1.We also provide the results from the asymptotic expansion expression 3 in [45], which have been calculated up to O(m 8 b /m 8 H ). Perfect agreements were found for both the NLO X y b y b 1 and NNLO X y b y b 2 coefficients.The reason is that the expansion parameter m 2 b /m 2 H ≈ 5 × 10 −4 is rather small.Actually, the difference between the expansion up to O(m 2 b /m 2 H ) and the exact result at the scale µ = m H is already 0.002% for X y b y b

Figure 9 :
Figure 9: Γ y b y b Hb b at different values of µ.The left (right) plot corresponds to the result with MS (on-shell) Yukawa couplings.The LO, NLO and NNLO predictions are represented by the blue, red, and green lines, respectively.
H → b b at O(y 2 b α 2 s ) with full dependence on m b .The result is obtained by calculating the imaginary part of the three-loop self-energy

Figure 10 :
Figure 10: Γ Ht t at different values of m H .The LO, NLO and NNLO predictions with QCD scale uncertainties are represented by the blue, red, and green bands, respectively.

Figure 11 : 1 z 1 z 1 z 1 z 10 + 1 z 15 − 1 z 3 60 − 1 z
Figure 11: Master integrals in the NP1 topology.The thick black and red lines stand for the massive bottom quark and the Higgs boson, respectively.One black dot indicates one additional power of the corresponding propagator.The numerators can be found in the text.

Table 2 :
H , 2m H ]. The scale uncertainty is significantly reduced down to about 9% and 3% at NLO and NNLO, respectively.The NLO correction increases the LO decay width by a factor of 12% − 28% depending on the scale.The NNLO correction is very small at µ = m H /2 but can be as large as 7% at µ = 2m H . Decay width of Γ Hb b (in MeV) at different orders in QCD.The second row shows the results with the MS Yukawa coupling while the third row shows the results with the on-shell Yukawa coupling.The renormalization scale uncertainties around µ = m H are also shown.The NNLO(y b y t ), N 3 LO(y 2 b ), and N 4 LO(y 2 b