Resummed prediction for Higgs boson production through $b\bar{b}$ annihilation at N$^3$LL

We present an accurate theoretical prediction for the production of Higgs boson through bottom quark annihilation at the LHC up to next-to-next-to-next-to leading order (N$^3$LO) plus next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy. We determine the third order perturbative Quantum Chromodynamics (QCD) correction to the process dependent constant in the resummed expression using the three loop bottom quark form factor and third order quark soft distribution function. Thanks to the recent computation of N$^3$LO corrections to this production cross-section from all the partonic channels, an accurate matching can be obtained for a consistent predictions at N$^3$LO+N$^3$LL accuracy in QCD. We have studied in detail the impact of resummed threshold contributions to inclusive cross-sections at various centre-of-mass energies and also discussed their sensitivity to renormalization and factorization scales at next-to-next-to leading order (NNLO) matched with next-to-next-to leading logarithm (NNLL). At N$^3$LO+N$^3$LL, we predict the cross-section for different centre-of-mass energies using the recently available results in \cite{Duhr:2019kwi} as well as study the renormalization scale dependence at the same order.

Contents coupling in the dominant decay channel of the Higgs boson to a pair of bottom quarks is a challenging task. Associated production of Higgs boson with vector bosons or with top quarks, its subsequent decay to bottom quarks is the promising one while there are also other proposals [8].
While the Higgs boson dominantly decays to bottom quarks in SM, its production from bottom quark annihilation is much smaller than the gluon initiated subprocess. However, as the precise measurement of the Higgs cross-sections is underway at the LHC, inclusion of bottom quark initiated channels in the theoretical predictions is unavoidable. Unlike in the SM, bottom quarks in the Minimal Supersymmetric SM (MSSM) [9] couple to neutral Higgs boson with the coupling proportional to 1/ cos β which in some parameter region can increase the production rate. Here, the angle β is related to the ratio, denoted by tan β, of the vacuum expectation values of two Higgs doublets. Hence, there is a considerable interest in studying the production mechanism of a single and a pair of Higgs bosons through bottom quark annihilation. The production of Higgs boson(s) in this mode is often studied using two approaches namely four flavor and five flavor schemes often called 4FS and 5FS respectively. In 4FS, the bottom quark distribution in the proton is set equal to zero, however they are radiatively generated through gluons in the proton allowing the possibility of their annihilation to produce the Higgs boson. Such contributions are enhanced by large logarithms that are proportional to bottom quark mass and hence they need to be resummed to get reliable predictions. The alternate approach, 5FS, avoids this enhancement through the introduction of non-zero bottom quark distributions in the proton. The origin of these distributions can be traced back to the bottom quarks resulting from gluon distributions inside the proton, thanks to the DGLAP evolution equation of parton distribution functions, which resums collinear enhanced logarithms to all orders in perturbation theory. While both these schemes should give the same result at the observable level, care is needed while comparing their predictions. The leading order contribution in 4FS comes from two to three scattering reaction g + g → b + b + H and the corresponding cross-section is of order α 2 s λ 2 . Here, the strong coupling constant α s is given by α s = g 2 s /4π and the Yukawa coupling by λ = m b /v where m b is the mass of the bottom quark and v, the vacuum expectation value. On the other hand, in 5FS, the leading order reaction is b + b → H and the corresponding cross-section is proportional to only λ 2 and is two orders less in the strong coupling constant. Hence, only at NNLO, the partonic cross-sections in the 5FS will have same order coupling strength, O(α 2 s λ 2 ), as compared to the leading order (LO) cross-section in the 4FS. Because of this, the higher order computation in perturbative QCD in the 5FS is relatively easier. Note that only next to leading order QCD effects [10][11][12] are known in 4FS while in 5FS, the state of the art N 3 LO prediction [1] is already available. The later computation provides an opportunity to compare N 3 LO predictions against those at next-to-leading order (NLO) computed in 4FS in a consistent manner. In 5FS, the complete NLO [13,14] and NNLO [15] as well as dominant threshold effects at N 3 LO [16,17] are known for quiet sometime. In this article we improve the 5FS crosssection by resumming the threshold logarithms up to N 3 LL accuracy. It was observed in [18,19] that the 5FS cross-section provides dominant cross-section in a matched prediction. Thus threshold improved result justifies in the context of precision study for this process. Recently the 5FS prediction has been also improved [20] by resumming time-like logarithms in SCET framework.
Fixed order (FO) predictions in perturbative QCD are often plagued with large logarithms resulting from certain boundaries of the phase space and hence their reliability in those regions are questionable. In the inclusive production rates, when partonic scaling variable z = m 2 h /ŝ → 1, which corresponds to the emission of soft gluons, large logarithms are generated at every order in perturbation theory. Here, m h is the Higgs boson mass andŝ is the partonic centre-of-mass energy. One finds a similar problem in the transverse momentum and rapidity distributions of Higgs boson when there are soft gluon emissions. This is resolved by resumming these large logarithms to all orders in perturbation theory through a systematic resummation approach. For inclusive rates, we refer [21][22][23] to the earliest approach. Catani and Trentadue, in their seminal work [22], demonstrated the resummation of leading large logarithms for the inclusive rates in Mellin space.
Needless to say that the inclusion of higher order QCD effects is of utmost importance to achieve precise prediction. In addition, these terms reduce the dependence on the unphysical scales such as renormalization and factorization scales at the observable level. Note that for the production of scalar Higgs boson through gluon fusion, the N 3 LO contribution is now known [24,25], which was further improved by the resummation of threshold logarithms, arising from soft gluon emissions, to N 3 LL accuracy [26][27][28] (the prime denotes that in addition to the N 3 LL terms, the resummed result includes higher logarithmic order terms coming from the matching to N 3 LO). Such an analysis leads to a precise determination of the SM Higgs cross-section at the LHC with small uncertainty.
The goal of this paper is to present the prediction for the inclusive production of Higgs boson in bottom quark annihilation at the LHC taking into account the resummed threshold corrections at next-to-next-to-next to leading logarithmic accuracy. We can achieve this using the recent N 3 LO prediction [1] and N independent threshold constant g b,0 computed to third order in QCD in this paper. Here N denotes the Mellin moment. The g (3) b,0 is obtained using the three loop form factor [29] of the Higgs-bottom-anti-bottom operator and the third order soft distribution function computed in [17].

Theoretical Framework
The Lagrangian that describes the Yukawa interaction of the SM Higgs boson with bottom quarks is given by is the bottom quark field and φ(x) is the SM Higgs field. λ is the Yukawa coupling given by λ = m b v . Note that we will use non-zero mass of the bottom quark only in the Yukawa coupling, elsewhere it is treated as massless quarks i.e. we use the VFS scheme throughout our analysis. The inclusive cross-section for the production of Higgs boson in proton collision is given by where f c (x, µ 2 f ) is the non-perturbative parton distribution function (pdf) with c denoting the parton type and x its momentum fraction. τ = m 2 h /S is the scaling variable, where S is hadronic centre-of-mass energy. The renormalization and factorization scales are denoted by µ r and µ f respectively. The born cross-section σ (0) bb (µ 2 r ) is given by (2. 3) The mass factorized parton level cross-section (∆ ab ) is calculable order by order in strong coupling constant, a s (µ 2 r ) = g 2 s (µ 2 r )/16π 2 in perturbative QCD, At each order in perturbation theory one can write with z = m 2 h /ŝ. In the above equation ∆ SV collects all those contributions that result from soft and collinear configurations of partons in the scattering events. They are proportional to distributions of the kind δ(1 − z) and D j (z) where The superscript SV is the short form of soft plus virtual. The remaining contribution is called the regular part of the cross-section denoted by ∆ reg ab . In QCD, ∆ ab was computed up to NNLO level in [15], at N 3 LO level in the threshold limit it was obtained in [16,17] and recently the complete N 3 LO result has been reported in [1].
The threshold contributions to inclusive cross-section ∆ ab (z) originate from soft and collinear partons in the virtual and real emission subprocesses. These contributions demonstrate remarkable factorization property through process independent cusp, soft and collinear anomalous dimensions. Consequently, the leading contributions resulting from the large threshold logarithms of the form D j (z) can be summed to all orders in a systematic fashion. Following [21][22][23], the resummation of these logarithms can be efficiently achieved in Mellin N -space and the resulting resummed threshold contribution takes the following form: where we have set µ r = µ f . The cusp anomalous dimension A (b) (a s ) and the constant D (b) (a s ) are process independent and hence can be obtained from the resummation result of Drell-Yan process. They are known up to third order in QCD and are listed in the Appendix C. The N independent constantg b,0 (a s (µ 2 f )) on the other hand is process dependent. It gets contribution from the process dependent virtual Higgs-bottom-anti-bottom quark form factor. The pre-factorg b,0 (a s (µ 2 f )), however gets modified due to the finite delta contribution from the soft function in Eq. (2.7) after integration. This leads to a new pre-factor g b,0 (a s (µ 2 f )). Following the method described in [16] and using the three loop form factor for the Higgs-bottom-anti-bottom quark [29] and the third order soft distribution function [30], we determine g b,0 (a s ) up to third order in strong coupling constant.
where C A = N and C F = (N 2 − 1)/2N are the Casimirs of SU (N ) and L f r = log(µ 2 f /µ 2 r ) and L qr = log(q 2 /µ 2 r ) with q = m h . In addition ζ i are the Riemann zeta functions. These constants g (i) b,0 can also be obtained from the soft-virtual results ∆ SV,(i) bb computed in [17] after taking the Mellin moment in N space ( collected in App. B ) and setting all the logarithms log(N ) and Euler Gamma γ E to zero.
To compute G N , we use the method described in [31] up to N 3 LL accuracy. Defining λ = 2β 0 a s (µ 2 r ) ln(N ), we obtain Thus, the resummed contribution to the Higgs production in bottom quark annihilation in Mellin N -space takes the simple form (2.14) The coefficients g i for i = 1, 2, 3, 4 are listed in the App. A. Note that except A (b) 4 both the cusp A (b) and collinear D (b) anomalous dimensions are known to third order in a s . Since the approximate A (b) 4 is available in the literature [33][34][35][36][37][38][39][40][41][42] we can readily predict N 3 LL contributions to inclusive Higgs production in bottom quark annihilation. In the next section, using the recently available predictions [1] for the fixed order N 3 LO contributions, we present the resummed prediction to N 3 LO + N 3 LL accuracy.

Phenomenology
In this section, we present a detailed discussion on the numerical impact of the resummed threshold contribution to the inclusive production of the Higgs boson in bottom quark annihilation at the LHC up to N 3 LO+N 3 LL accuracy in perturbative QCD. The resummed part of the inclusive cross-section for the Higgs production in N space can be obtained by taking N -th Mellin moment of Eq. (2.2) as r , µ 2 f ) whose N -th moment ∆ res N is given in Eq. (2.7). The resummed cross-section is then added to the fixed order one after subtracting the Mellin moment of ∆ SV,(n) (z) in the large N limit. This is done because they are already present in the fixed order result and hence this will avoid double counting. This is achieved through a matching procedure at every order. Finally, the matched result takes the following form: In the above equation the superscript N n LL in ∆ res implies that we retain up to g (n) b,0 terms in the g b,0 and up to g n+1 (λ) in the exponent G N given in Eq. (2.13). Similarly, N n LO implies that we retain the fixed order result up to order a n s . The subscript tr in the last term of the above equation means truncation of the series in a s to desired accuracy.
The fixed order analytical results [15] up to NNLO order have been implemented in a fortran code. We have validated our predictions to a very good accuracy with the available public code SuShi [43]. For N 3 LO+N 3 LL analysis, we have used the numbers presented in Tab. II in the arXiv version of the paper [1] for N 3 LO. We perform the inverse Mellin transformation of the resummed N -space result using an in-house fortran code. Minimal prescription [44] has been used to deal with the Landau pole in the Mellin inversion routine. Since we work in the 5FS, we take n f = 5 throughout. We use the MMHT2014(68cl) parton distribution set [45] and renormalization group (RG) running for a s at each order. One could in principal use the strong coupling as provided through LHAPDF [46] interface, which at the NNLO level changes the cross-section by 0.09%. The Yukawa coupling is evolved using 4-loop RG with bottom mass m b (m b ) = 4.3 GeV. It is well known that the optimal choice for the central scale to study the Higgs production in bottom quark annihilation is µ r = m h and µ f = m h /4. This choice mimics most of the higher order contributions. Hence, we have made this choice throughout. In addition, we have predictions for other central scale choice namely µ r = µ f = m h .
[pb] In the left panel of Fig. 1, we present the resummed cross-section up to NNLO+NNLL level against the hadronic centre-of-mass energy for 7 TeV to 100 TeV. The bands in the plot correspond to the scale variations obtained by varying the unphysical renormalization and factorization scales in the range (1/2, 2)(µ defined as the ratio of the cross-section at a particular order (NLO+NLL, NNLO+NNLL) over the same at the LO+LL order. At NLL the k-factor increases by 14% and at NNLL by 21% compared to LL for 13 TeV LHC. We find that the uncertainty due to µ r and µ f scales increases with the energy of the collider, however at the current energy of the LHC it is within 11% at NNLL. The reason for the large uncertainty at high E CM could be due to the lack of knowledge of the pdf sets at these energies. In the right panel of Fig. 1, we compare the NNLO and NNLO+NNLL predictions for various centre-of-mass energies. We find that the resummed cross-section decreases by 3% over the fixed order one for the entire range of centre-of-mass energies. The scale uncertainty however does not improve much compared to NNLO.  In Fig. 2, we present our findings on the sensitivity of our predictions to the renormalization and factorization scales at the LHC energy 13 TeV. We plot the cross-sections at various orders as a function of µ r and µ f in the range (0.1, 10)m h . In the later part of our analysis we will choose the range (1/2, 2)m h . In the first panel, we vary µ f keeping µ r = m h . We find that predictions for the cross-sections from fixed order as well as resummed ones are very sensitive to µ f at first two orders and become better at third order. This could be due to lack of precise knowledge on the bottom quark distribution function. At NNLO level, the uncertainty due to µ f is about −21 −16 % in the range µ f = (0.1, 10)m h and at NNLO+NNLL this goes down to −21 −8 %. In the second panel we vary µ r keeping the factorization scale µ f at m h /4. Interestingly, predictions from both the fixed order and the resummed one are less sensitive to µ r compared to µ f . In addition, as we increase the order of perturbation, the sensitivity to µ r goes down significantly. At LO, the cross-section de-creases by 26% when µ r is varied ten times its central value. On the other hand it increases rapidly as we decrease µ r . We find 59% increase in the cross-section when the scale is taken one tenth of the central value. In the resummed case, the corresponding changes are −28% +67% . At NLO level the fixed order cross-section changes by −16% +21% whereas for NLO+NLL case it improves and changes by −15% +19% . The corresponding numbers at NNLO are −6% −11% for FO and −6.7% −4.9% for resummed one. Finally in the third panel, we set µ r = µ f = µ and vary µ to see the combined effect. We observe that the resultant uncertainty from the renormalization and factorization scales for resummed result compared to the one coming from fixed order result does not show reduction. This could be due to the fact that the resummation takes into account only the bottom quark initiated channels beyond N 3 LO accuracy and in addition the factorization scale dependence dominates over the renormalisation scale. The inclusion of the other channels and better determination of the bottom quark distribution can lead to the reduction of scale uncertainty.
[pb] In Fig. 3, we present the resummed cross-sections against the fixed order ones for 13 TeV LHC. The uncertainty at each order has been obtained by using the 7-point variation i.e., by varying µ r and µ f scales around their central values m h and m h /4 respectively in the range (1/2, 2). At NNLL level the uncertainty is comparable to the NNLO prediction. However the central value of the resummed result shows a better perturbative convergence compared to the FO result. We find that the cross-section is larger by 38% compared to that at LO whereas at NNLO, there is an increase of 5.2% over the NLO prediction. In the resummed case, the NLO+NLL cross-section is 14.8% larger than the LO+LL cross-section, while NNLO+NNLL is 5.8% larger than NLO+NLL one. The perturbative convergence is seen to improve in the resummed case whereas the scale uncertainty does not improve compared to the fixed order. This is mainly due to strong dependence of the cross-sections on the factorization scale at every order in perturbation theory. However, we expect that the uncertainty from µ r will be mild.
There are several pdfs that exist in the literature and the numerical predictions on the cross-section do depend on the choice of pdf. A detailed analysis with respect to the choice of pdfs is thus required. In particular, it will help us to understand how the underlying theoretical assumptions and models in the pdf fits affect the b-quark pdf determination [47]. We use some of the widely used pdf sets; e.g. the ABMP16 5 nnlo, CT14nnlo, MMHT2014nnlo68cl, NNPDF31 nnlo as 0118, PDF4LHC15 nnlo 100 pdf sets and corresponding strong coupling through LHAPDF interface. The central values for all pdf sets are given by the zeroth subset of the pdfs except from the NNPDF where the central value is given by the mean corresponding to all subsets. The results are tabulated in Table-1 with the central values and corresponding pdf uncertainties for different pdf groups. We see that the central value as well as the uncertainties are small at NNLO+NNLL compared to the NNLO order for all the pdf sets. The difference in the central prediction among different pdf groups can be traced back to the different treatment of bottom quark pdf by different groups. We notice that the largest pdf uncertainty of 7.3% comes from PDF4LHC sets whereas NNPDF sets give least uncertainty of 1.3%. We also observe that the pdf uncertainty is marginally improved in the case of NNLO+NNLL compared to that of NNLO. Recently, the fixed order predictions for ∆ ab to third order, namely ∆ (3) ab have become available [1] and this allows us to predict the matched cross-section to N 3 LO+N 3 LL accuracy using the approximate A (b) 4 . Note that SV contribution to N 3 LO cross-section is already available, see [17]. We have used same set of input parameters as in [1] for our study. For example, we use the PDF4LHC2015nnlo set with strong coupling evolved at four loops as in [1]. The bottom quark mass is taken to be m b (m b ) = 4.18 GeV. At the LHC with centre-of-mass energy 13 TeV, we find that the threshold enhanced SV part of N 3 LO cross-section differs from the complete N 3 LO result by around 2%. Our prediction for the N 3 LO+N 3 LL cross-section at the central scale is 0.537 pb which differs from the fixed order contribution by around −1%. In Table- A detailed study on the dependence of the cross-section at N 3 LO+N 3 LL level on the renormalization and factorization scales is possible provided the numerical code at N 3 LO [1] is publicly available. In addition, for the µ f variation, one also needs the pdfs at the N 3 LO. However, using the renormalization group equation and the numerical predictions given in table-II of [1] for a selected set of values of centre-of-mass energies, it is possible to study the dependence on renormalization scale. In Fig. 4, we show how predictions at N 3 LO and N 3 LO+N 3 LL depend on the renormalization scale. In fact, entire dependence from µ r at N 3 LO level is controlled by the contributions from previous orders, namely LO, NLO and NNLO. We observe that the sensitivity to renormalization scale goes down at N 3 LO level compared to the previous order. At NNLO level the µ r variation changes the cross-section by +0.8 −1.3 % compared to the central value for the variation in the range µ r ∈ {1/2, 2}m h . At N 3 LO level this decreases to −0.5 +0.6 %. For the resummed prediction, we find that the dependence on µ r goes further down at N 3 LO+N 3 LL level. The corresponding uncertainty is +0. 3 −0.2 % compared to the central value. This shows the improvement in the stability against µ r at N 3 LO+ N 3 LL level compared to fixed order.

Conclusions
In this article we have studied the role of resummed threshold corrections to the inclusive cross-section for producing the Higgs boson through bottom quark annihilation at the LHC. While this is a sub-dominant channel compared to gluon fusion subprocess, the precise measurements that are carried out at the LHC in the context of processes involving Higgs bosons demand the inclusion of this channel for precision studies. Complete NNLO QCD corrections [15] and soft plus virtual corrections [17] to N 3 LO level to this observable are known for a while. More recently, a complete N 3 LO contributions resulting from all the partonic channels where the Higgs boson couples to bottom quarks became available [1]. In hadronic colliders, soft gluons play an important role in most of the observables. In the fixed order perturbative computations, in certain kinematic regions the soft gluons dominate. The large logarithms resulting from these soft gluons often spoil the reliability of the fixed order predictions. The resolution to this problem is to resum these logarithms to all orders in perturbation theory. The framework to resum such logarithms in a systematic fashion to all orders in perturbation theory for inclusive observables is well established and results based on this demonstrate better and reliable predictions compared to the fixed order ones. In the present article we have done a detailed study to understand the role of the soft gluons on the inclusive cross-section for producing Higgs boson in bottom quark annihilation. This is achieved within the frame work of threshold resummation in Mellin-N space. We have done this to N 3 LO+N 3 LL accuracy. We have used the recent prediction at N 3 LO level from [1] for the fixed order contribution and for the resummed part at N 3 LL level, except the process dependent constant g (3) b,0 , rest of the ingredients are already known. We have computed this constant for the first time using the three loop form factor [29] and the quark soft distribution function [17] known to third order in QCD. Our numerical result at N 3 LO+N 3 LL in QCD is the most precise prediction for the inclusive cross-section for the production of Higgs boson through bottom quarks at the LHC. We have predicted the inclusive rates at various centre-of-mass energies with the corresponding k-factors. In addition, we studied in detail the sensitivity of the resummed predictions to the renormalization and factorization scales and choice of parton distribution functions in order to estimate the theoretical errors precisely.
A Resummation constants g i (λ, µ 2 r , µ 2 f ) The resumamtion constants g i (see 2.13) in the N -space are found to be B Soft-Virtual coefficients in N -space Here we have collected all the large N coefficients for this process. DefiningL = ln N + γ E they are given below:

C The Cusp and the soft anomalous dimensions
The quark cusp anomalous dimensions are given as [48],