Higgs Rapidity Distribution in $b {\bar b}$ Annihilation at Threshold in N$^{3}$LO QCD

We present the rapidity distribution of the Higgs boson produced through bottom quark annihilation at third order in QCD using the threshold approximation. We provide a framework, based on the factorization properties of the QCD amplitudes along with Sudakov resummation and the renormalization group invariance, that allows one to perform the computation of the threshold corrections in a consistent, systematic and accurate way. The recent results on threshold N$^3$LO correction in QCD for the Drell-Yan production and on three loop QCD correction to Higgs form factor with bottom anti-bottom quark are used to achieve this task. We also demonstrate the numerical impact of these corrections at the LHC.


Introduction
With the spectacular discovery of the Higgs boson at CERN LHC [1,2] , the full spectrum of matter particles and force carriers of the Standard Model (SM) has been established very successfully. Though the mass of the newly discovered boson is already pinned down with an impressive experimental uncertainty of just a few hundred MeV in the range of 125 -126 GeV, to fully validate the mechanism of electroweak symmetry breaking and to shed light on the possible potential deviations from its SM apprehension, it is indispensable to study the inclusive as well as exclusive observables associated with the production and decay channels of the Higgs boson to a very high accuracy. Within the framework of SM, the production mechanism of the Higgs boson is dominated by gluon fusion, whereas one of the alternative channels, namely, bottom quark annihilation is severely suppressed by the small Yukawa coupling of bottom quark to the Higgs boson. However, in extensions of the SM with an enlarged spectrum of Higgs sector, as in the case of two-Higgs doublet model, the Yukawa coupling of bottom quark to some of the Higgs bosons can be enhanced significantly, such that the production channel of bottom quark annihilation could be the dominant one. Moreover, the contribution from gluon fusion channel decreases due to enhanced negative top-bottom interference diagrams. Furthermore, the bottom quark initiated processes at hadron colliders are of much theoretical interest on account of the freedom in treating the initial state bottom-quarks. In the four flavor scheme (4FS), alternatively known as the fixed flavor number scheme (FFS), the mass of the bottom quarks is considered to be non-zero throughout and they are excluded from the proton constituents, whereas, in the framework of five flavor scheme (5FS), also known as the variable flavor number scheme (VFS), the bottom quarks are considered as massless partons, except in the Yukawa coupling, with their own parton distribution functions (PDF).
While the theoretical predictions at NNLO and next-to-next-to leading log (NNLL) [26] QCD corrections and of two loop electroweak effects [27][28][29][30][31][32] played an important role in the discovery of the Higgs boson, the theoretical uncertainties resulting from the unphysical factorization and renormalization scales are not fully under control. In addition, the interpretation of the experimental data with higher accuracy from the upcoming run at the LHC demands the inclusion of higher order terms in QCD in the theoretical computation. Hence, the efforts to go beyond NNLO are going on intensively in past few decades. The computation of N 3 LO corrections is underway and some of the crucial ingredients, like the quark and gluon form factors [33][34][35][36][37], the mass factorization kernels [38] and the renormalization constant [39] for the effective operator describing the coupling between the Higgs boson and the SM fields in the infinite top quark mass limit are available up to three loop level in dimensional regularization. In addition, NNLO soft contributions are also known [40] in n dimensions. These results were already used to compute the partial threshold contributions at N 3 LO to the production cross-section of di-leptons in Drell-Yan (DY) and of the Higgs boson in gluon fusion as well as in bb annihilation, see [41][42][43][44][45]. Since then, there have been several advances [46][47][48][49] towards obtaining the complete N 3 LO result for the inclusive Higgs production. The milestone in this direction was achieved by Anastasiou et al. in [50] to obtain the complete threshold N 3 LO corrections. This result provided a crucial input in [51] to obtain the corresponding N 3 LO threshold corrections to DY production. Independently, in [52], using light-like Wilson lines threshold corrections to the Higgs boson as well as Drell-Yan productions up to N 3 LO were obtained. Catani et al. in [53] used the universality of soft gluon contributions near threshold and the results of [50] to obtain general expression of the hard-virtual coefficient relevant for N 3 LO threshold as well as threshold resummation at next-to-next-to-next-to-leading-logarithmic (N 3 LL) accuracy for the production cross section of a colourless heavy particle at hadron colliders. There have been several attempts to go beyond threshold corrections [54,55] for the inclusive Higgs production at N 3 LO. Recently, [56], the full next to soft as well as the exact results for the coefficients of the first three leading logarithms at this order have been obtained for the first time. For the Higgs boson production through bb annihilation, the recent results of the Higgs form factor with bottom-antibottom by Gehrmann and Kara [57] and the universal soft distribution obtained for the Drell-Yan production [51] enabled us to obtain the missing δ(1 − z) contribution (see [44,45,58] for the partial results to this order) to the production cross-section at threshold at N 3 LO [59].
Like the inclusive one, the differential rapidity distributions are computed for the dilepton pair in DY [60] and the Higgs boson produced through gluon fusion in [61,62], the Higgs boson through bb annihilation in [63] and associated production of the Higgs with vector boson in [64,65] to NNLO in QCD. Using the formalism developed in [44,45], the partial N 3 LO threshold correction to the rapidity distributions of the dileptons in DY and the Higgs boson in gluon fusion as well as bottom quark annihilation were computed in [66]. Following the same technique, we obtained the complete N 3 LO threshold correction to the rapidity distributions of both dilepton pair in DY and the Higgs boson in gluon fusion [67]. We had seen the dominance of the threshold contribution to the rapidity distribution in these processes. A significant amount of reduction in the dependence on the unphysical renormalization and factorization scale of the rapidity distribution takes place upon inclusion of the N 3 LO threshold corrections. In addition, these computations provide first results beyond NNLO level and will serve as a non-trivial check for a complete N 3 LO results. Keeping these motivations in mind, we intend to extend the existing result of the rapidity distribution of the Higgs boson produced through bb annihilation to higher accuracy, namely the inclusion of complete N 3 LO threshold correction.
In Sec. 2.1, we perform an explicit calculation of threshold correction to the rapidity distribution of the Higgs boson in bb annihilation at NLO, using the factorization properties of QCD amplitude, Sudakov resummation of soft gluons and renormalization group invariance. This helps us to build an elegant framework to calculate the rapidity distribution at threshold, of a colorless state produced at hadron colliders, to all orders in QCD perturbation theory. In Sec. 2.2, we use that general framework to achieve the goal of computing the complete analytic expression for the threshold corrections beyond NLO and provide the result up to N 3 LO. Sec. 2.3 contains the discussion on the numerical impacts of our results. Finally, we conclude with our findings in Sec. 3.

Differential Distribution with Respect to Rapidity
The interaction of bottom quarks and the Higgs boson is encapsulated in the following action where, ψ b (x) and φ(x) denote the bottom quark and scalar field, respectively. The Yukawa coupling λ is given by √ 2m b /v, with the bottom quark mass m b and the vacuum expectation value v ≈ 246 GeV. Throughout our calculation, we consider five active flavours (VFS scheme), hence except in the Yukawa coupling, m b is taken to be zero like other light quarks in the theory.
We study infrared safe differential distribution, namely rapidity distribution of the Higgs boson at hadron colliders, in particular those produced through bottom anti-bottom annihilation. Our findings are very well suited for similar observables where the rapidity distribution is for any colorless state produced at hadron colliders. We will set up a framework that can provide threshold corrections to rapidity distribution of the Higgs boson to all orders in perturbation theory. It is then straightforward to obtain fixed order perturbative results in the threshold limit.
The general frame work that we set up for the computation of threshold corrections beyond leading order in the perturbation theory for such observables is based on the factorization property of the QCD amplitudes. Sudakov resummation of soft gluons, renormalization group equations and most importantly the infrared safety of the observable play important role in achieving this task. QCD amplitudes that contribute to hard scattering cross sections exhibit rich infra-red structure through cusp and collinear anomalous dimensions due to the factorization property of soft and collinear configurations. Massless gluons and light quarks are responsible for soft and collinear singularities in these amplitudes and also in partonic subprocesses. Singularities resulting from soft gluons cancel between virtual and real emission diagrams in infrared safe observables. While the final state collinear singularities cancel among themselves if the summation over degenerate states are appropriately carried out in such observables, the initial state collinear singular configurations remain until they are absorbed into bare parton distribution functions. In the upcoming section, we present one loop computation for the rapidity distribution in order to demonstrate how the various soft singularities cancel and also to give a pedagogical derivation of how the most general resummed threshold correction to the rapidity distribution can be obtained.

Threshold Correction at NLO
The process under consideration is the production of the Higgs boson through bottom quark annihilation in hadron colliders. The leading order process is where, k i 's are the momenta of the incoming bottom and anti-bottom quarks involved in partonic reaction and q is the momentum of the Higgs boson. The hadronic center of mass energy squared is defined by S ≡ (p 1 + p 2 ) 2 , where p i 's are the hadronic momenta and the corresponding one for the incoming partons is given asŝ = (k 1 + k 2 ) 2 . The fraction of the initial state hadron momentum carried by the parton is denoted by x i i.e. k i = x i p i . The rapidity of the Higgs boson is defined through The differential distribution with respect to rapidity of the Higgs boson can be expressed as with τ ≡ q 2 /S, q 2 = m 2 H , m H -the mass of the Higgs boson. λ(µ 2 R ) is the Yukawa coupling defined at the renormalization scale µ R , N = 3 is the number of QCD colors and σ b,(0) is the leading order cross-section. Defining z ≡ q 2 /ŝ, we find In this expression, X is the remnants other than the Higgs boson, Z b (µ 2 R ) is the ultraviolet (UV) renormalization constant for the Yukawa coupling λ and dP S 1+X is the phase space element for the H + X system. M ac→H+X denotes the scattering amplitude at partonic level. The functionĤ ac (x 1 , x 2 ) is the product of unrenormalized parton distribution functions (PDF)f a (x 1 ) andf c (x 2 ), (2.6) The PDF f a (x 1 , µ 2 F ), renormalized at the factorization scale µ F , is related to the unrenormalized ones through Altarelli-Parisi (AP) kernel Γ ad as follows: where, the scale µ is introduced to keep the strong coupling constantĝ s dimensionless in space-time dimensions n = 4 + , regulating the theory andâ s ≡ĝ 2 s /16π 2 . Expanding the AP kernel in powers ofâ s , we get where, P ad (z) is the leading order AP splitting function. S = exp (γ E − ln 4π) 2 where γ E is the Euler-Mascheroni constant. Using Γ ad , W b can be written in terms of renormalized H, given by (2.9) The LO contribution arises from the Born process b + b → H and the NLO ones are from one loop virtual contributions to born process and from the real emission processes, namely . For LO and virtual contributions, dP S 1+X = dP S 1 and for real emission processes we have two body phase space element dP S 1+X = dP S 2 . In order to define the threshold limit at the partonic level and to express the hadronic cross-section in terms of the partonic one through convolution integrals, we choose to work with the symmetric scaling variables x 0 1 and x 0 2 instead of y and τ which are related through In terms of these new variables, the partonic subprocess contributions can be shown to depend on the ratios z j = x 0 j x j which take the role of scaling variables at the partonic level. The dimensionless partonic differential cross-section denoted by∆ b d,ac through 1 is UV finite. Here subscript d stands for differential distribution. The collinear singularities that arise due to the initial state light partons are removed through the AP kernels resulting in the following finite ∆ b d,ac (2.12) Therefore, expressing W b in terms of renormalized H ac and finite ∆ b d,ac , we get Since, W b involves convolutions of various functions, it becomes normal multiplication in the Mellin space of the Mellin moments of renormalized PDFs, AP kernels and bare differential partonic cross-section. The double Mellin moment of The threshold limit is defined by N i → ∞, which in z j variables corresponds to z j → 1. In this limit, only diagonal terms in the AP kernel Γ −1 and ∆ b d contribute to the differential cross-section. Hence, ln∆ b d is simply a sum of the contributions from 1) diagonal terms of the AP kernels and 2) bare differential partonic cross-section. Due to the born kinematics, the form factor contribution can be further factored out from the differential partonic cross sections to all orders in perturbation theory. Hence, the remaining part of the differential partonic cross-sections contains contributions from only real emission processes, namely those involving only soft gluons. Taking into account the renormalization constant of the Yukawa couplingλ, we find whereF b and S b (N 1 , N 2 ) are bare form factor and real emission contributions of partonic subprocesses, respectively. The inverse Mellin transform will bring back the expressions in terms of the variables z j and they will contain besides regular functions, the distributions namely δ(1 − z j ), D i and D i , defined as (2.17) The subscript '+' denotes the customary 'plus-distribution' f + (z) which acts on functions regular in z → 1 limit as where, g(z) is any well behaved function in the region 0 ≤ z ≤ 1. In the threshold limit, we drop all the regular terms and keep only these distributions.
In the following, we perform NLO computation in the threshold limit. The overall renormalization constant (Z b ) 2 is found to be The form factor contribution |F b | 2 at one loop level gives The contribution from Γ bb in the threshold limit is found to be 1 z 1 Note that the regular terms in the limit z j → 1 in Γ bb do not contribute in the threshold limit and hence dropped. The inverse Mellin transform ofS b (N 1 , N 2 ), namely S b (z 1 , z 2 ) can be obtained directly from the real gluon emission processes in bottom anti-bottom annihilation processes: b + b → H + g. The two body phase space is given by The phase space in the limit z j → 1 becomes The spin and color averaged matrix element square in threshold limit is found to be where terms that are regular in z j as z j → 1 have been dropped. It is then straightforward to obtain the threshold contribution resulting from the real gluon emission process: Using the identity it can be shown that ∆ b d (z 1 , z 2 ) in the threshold limit contains only the distributions such as δ(1 − z j ), D i and D i . Decomposing ∆ b d,ac into hard and soft parts, At the hadronic level, decomposing W b as where all the parton densities are defined at µ F = m H . In general, The Spence function (li 2 (x)) is defined as The exact result computed at NLO level confirms our expectations, see for example [68] where the rapidity distribution of di-leptons in the Drell-Yan production for a physics beyond the SM (BSM) involving a generic Yukawa type interaction was obtained to NLO level. After the suitable replacement of the BSM coupling in [68], we obtain where W 's can be expanded in the strong coupling constant a s (µ 2 F ) as and the corresponding coefficients are given by and (2.37) In the threshold limit, after setting µ F = m H , we find that the above result reduces to given in Eq. 2.30.

Threshold Corrections Beyond NLO
Following the factorization approach that we used in the previous section to obtain the threshold correction to NLO rapidity distribution, we now set up a framework to compute threshold correections to rapidity distribution to all orders in strong coupling constant. Our approach is based on the fact that the rapidity distribution in the threshold limit can be systematically factorized into 1) the exact form factor, 2) overall UV renormalization constant, 3) soft gluon contributions from real emission partonic subprocesses and 4) the diagonal collinear subtraction terms involving only δ(1−z) and D 0 (z) terms of AP splitting functions. We call such a combination soft-virtual (SV) part of the rapidity distribution and the remaining part as hard. Hence, we propose that The symbol 'C' means convolution with the following definition where, ⊗ indicates double Mellin convolution with respect to the variables z 1 and z 2 and the function f (z 1 , z 2 ) is a distribution of the kind δ(1−z j ) and/or D i (z j ). The finite distribution ) and the mass factorization kernels Γ bb (â s , µ 2 , µ 2 F , z j , ): We have expressed all the quantities in the above equation in terms of unrenormalized strong coupling constantâ s related to the standardα s throughâ s =α s /4 π and the dimensional regularization scale µ. The UV renormalization ofâ s is done at the renormalization determines the structure of the Z(µ 2 R ), up to O(a 3 s ), we find (2.44) The first three coefficients of the QCD β function, β 0 , β 1 and β 2 are given by [69] with the SU (N ) color factors and n f is the number of active flavours. The overall operator renormalization constant Z b renormalizes the bare Yukawa couplingλ resulting λ(µ 2 R ) through the relation with the anomalous dimensions γ b i given by [70][71][72] Upon solving the above RGE in 4 + space-time dimensions, we obtain up to O(a 3 s ). The bare form factorF b (â s , Q 2 , µ 2 , ) satisfies the following differential equation which follows from the gauge as well as renormalization group invariances [73][74][75][76] where, all the poles in are encapsulated within K b and G b contains the terms finite in .
where, A q i 's are the cusp anomalous dimensions, found to be [38,[77][78][79][80] A q 1 = 4C F , Being flavor independent, A b i 's are same as A q i . Solving the RGE 2.52 satisfied by K b we get Similarly upon solving the RGE 2.52 for G b , we obtain Expanding the finite function G b (a s (Q 2 ), 1, ) in powers of a s (Q 2 ) as Note that the single pole term of the form factor depends on three different anomalous dimensions, namely the collinear anomalous dimension B q i , anomalous dimension of the coupling constant γ b i and the soft anomalous dimension f q i . B q i can be obtained from the δ(1 − z) part of the diagonal splitting function known up to three loop level [38,77] which are (2.59) The f q i for i = 1, 2 can be found in [81] and in [38] for i = 3. We list them below: (2.60) Since B q i and f q i are flavour independent, we have used i are controlled by the beta function of the strong coupling constant through renormalization group invariance of the bare form factor: The coefficients g b,k i can be extracted from the finite part of the form factor. Up to two loop level, we use [19,44,45] and at three loop level the recent computation by Gehrmann and Kara [57] enable us to compute the relevant g b, 1 3 in [59] where g b,1 3 was already used to obtain threshold correction to inclusive Higgs production in bottom anti-bottom annihilation process: Using the expressions for K b and G b given in Eq. 2.54 and Eq. 2.58, respectively, we obtain the renormalized form factor up to order O(a 3 s ) as (2.63) Note that the poles of ln |F b | 2 are fully controlled by the universal anomalous dimensions A q , γ b , B q and f q while the constant terms require vertex dependent constants g b,k i . In MS scheme, the mass factorization kernels Γ bb (â s , µ 2 , µ 2 F , z j , ) remove the collinear singularities which arise due to massless partons. These kernels satisfy the following RG equation : where P bc z j , µ 2 F are AP splitting functions. We can expand the P bc z j , µ 2 F in powers of a s as The off diagonal splitting functions are regular as z j → 1. The diagonal ones contain in addition distributions such as δ(1 − z j ) and D 0 multiplied by the universal anomalous dimensions B q i and A q i , respectively: As we are interested in results from the threshold region, we can ignore all the non-diagonal splitting functions and also the regular part P (i) reg,bb arising from the diagonal terms. Hence, the solution to Eq. 2.64 takes the following form: (2.67) Finally, we need to determine the soft distribution function Φ b d (â s , q 2 , µ 2 , z 1 , z 2 , ) in ∆ SV d,b . Its most general form can be systematically constructed if Φ b d also satisfies a differential equation similar to the form factor. It is indeed the case because the q 2 dependence and pole structure of Φ b d have to be similar to those of ln |F b | 2 in order to obtain finite distribution Ψ in the limit → 0 [44,45]. Hence, we propose that Φ b d satisfies and consequently The right hand side of the above equation is proportion to δ(1 − z 1 )δ(1 − z 2 ) as the most singular terms resulting from K b d should cancel with those from the form factor contribution which is proportional to only pure delta functions. To make the ∆ SV d,b finite, the poles from Φ b d (â s , q 2 , µ 2 , z 1 , z 2 , ) have to cancel those coming fromF b and Γ bb . Hence the constants A q should satisfy The RGE 2.70 for G b d can be solved using the above relation to get With these solutions, it is now straightforward to solve the above differential equations 2.68 for Φ b d to get where the constants G b d,i ( ) are flavour independent and they satisfy the following structure similar to G b i ( ) of the form factor, i.e., (2.79) In the above expression, we have used G d,1 are needed to obtain Φ b d . We achieve this using the following identity: where σ b is known to NNLO level [19] exactly and to N 3 LO level in the threshold limit [59]. In large N limit i.e. N → ∞ the above Eq. 2.80 relatesφ q,(i) d ( ) toφ q,(i) ( ) that appears in inclusive threshold corrections to Drell-Yan process (see [44,45,51,59]) as followŝ We present below our results for ∆ At the stage, we can demonstrate that integration over the rapidity correctly reproduces inclusive threshold contribution to the Higgs production in bottom anti-bottom annihilation reported in [59] : The integration over the rapidity y leads to the following relation between ∆ SV d,b (z 1 , z 2 ) obtained in this paper and ∆ SV b (z) in [59]: We have explicitly checked that the results presented here for ∆ SV d,b and those for ∆ SV b in the [59] up to N 3 LO level satisfy the above relation confirming the consistency of the formalism used. For completeness, we present the results for ∆ SV,(i) d,b up to N 3 LO after substituting all the constants that are required to this order:

Numerical Results
In this section, we present the numerical impact of the rapidity distribution of the Higgs boson, produced via bottom anti-bottom annihilation subprocess at the LHC. The rapidity distribution can be expanded in powers of the strong coupling constant a s as (2.92) Beyond LO, the distribution is split into hard and SV parts as In the following, for our numerical study we will use the exact results up to NLO level but at NNLO, we use exact NLO and only threshold contribution at O(a 2 s ) as we do not have access to the hard part at O(a 2 s ) computed in [63] 1 . We call it NNLO(SV). Similarly at N 3 LO level, we will use NNLO(SV) and threshold contribution at O(a 3 s ), denoted by N 3 LO(SV) hereafter. We present results for the center of mass energies 8 and 13 TeV at the LHC. The standard model parameters which enter into our computation are the Z boson mass m Z = 91.1876 GeV, top quark mass m t = 173.4 GeV and mass of the Higgs boson m H = 125 GeV. The strong coupling constant is evolved using the 4-loop RG equations with α N 3 LO s (m Z ) = 0.117. Following the Ref. [84], the solution to RGE 2.47 for λ(µ 2 R ) is given by, The K i are given by and where µ 0 is some reference scale at which λ is known. We have numerically evaluated λ(µ 2 R ) to relevant order namely LO, NLO, NNLO and N 3 LO by truncating the terms in the RHS of Eq. 2.47. We have used λ(µ 2 0 ) = √ 2m b (µ 0 )/v and m b (µ 0 ) = 3.63 GeV with the choice µ 0 = 10 GeV. We use the MSTW2008 [85] parton density sets with errors estimated at 68% confidence level with five active flavours. Parton densities and α s are evaluated at each corresponding perturbative order. Specifically, we use (n + 1)-loop α s at N n LO, with n = 0, 1, 2, 3. However, we use MSTW2008NNLO PDFs at N 3 LO, the N 3 LO kernels not being available at the moment. We set the renormalization scale µ R = m H and factorization scale µ F = m H /4 [16] as their central values.
Several checks have been performed on our numerical code. We have found complete agreement with the literature on the inclusive Higgs production rate [19,59] after performing an additional numerical integration over the rapidity Y of our distribution. The check was also performed at the analytical level. However, we were not able to reproduce the   plot given in [63], after using the same set of values of the input parameters. We begin our discussion with the results at NLO level. In Sec. 2.1, we presented the contributions coming from the exact results, containing the regular as well as pure threshold ones to the rapidity distribution at O(a s ). In Fig. 1, we plot both the NLO(SV) and exact NLO rapidity distri- butions to exhibit the dominance of threshold over the entire rapidity range after setting the values of the renormalization and factorization scales to their central values. From now onward, we adopt a consistent representation to display the figures corresponding to our results. In every figure, the left panel shows the result for 8 TeV whereas the right panel corresponds to 13 TeV at the LHC. We observe that the exact NLO contribution is well approximated by the NLO(SV), thanks to the intrinsic property of the matrix element, where the phase-space points corresponding to the born kinematics contribute towards the largest radiative corrections for the low τ (m 2 H /s ≈ 10 −4 ) values. So, we expect that the trend of approximating the exact results by threshold corrections at that order to remain same after the inclusion of higher order terms also.
With this in mind, we present the results at LO, NLO, NNLO(SV), N 3 LO(SV) for different values of the rapidity Y after setting the central values for renormalization and factorization scales for 8 TeV in Table 1 and for 13 TeV in Table 2 at LHC. The hadronic cross-section, obtained by the convolution of the partonic cross section with the parton densities, suffers from the theoretical uncertainties, arising from the missing higher order corrections, through the renormalization (µ R ) and factorization (µ F ) scales. These can be estimated through the variation of the differential hadronic cross section with µ R and µ F , thereby exhibiting the size of the higher order effects.
In Fig. 2, we plot two curves for each order for the predictions at NLO, NNLO(SV), N 3 LO(SV) corresponding to two different choices of the renormalization scale, µ R = 0.1m H and µ R = 10m H , keeping the factorization scale fixed at µ F = m H /4, whereas in Fig. 3, we plot the predictions at each order corresponding to two different choices of the factorization scale, µ F = m H /8 and µ F = m H /2, keeping the renormalization scale fixed at µ R = m H . We observe a consistent improvement in the accuracy of the predictions with the inclusion of the higher order terms, the width of the bands being an clear indicator of the theoretical uncertainties. Moreover, we can see that the dependence on the renormalization scale for this process is very mild. Another way to assess the reliability of the prediction is to study the rate of convergence of the perturbation series, represented by the K-factor.
In the Fig. 4, we plot the K-factors defined as K 1 = dσ N LO /dσ LO and K (SV) i = dσ N i LO(SV ) /dσ LO , i = 2, 3 as a function of Y . For 8 TeV LHC, we see that the K 1 varies from 1.57 to 0.63 over the entire rapidity range, while the value of K 1 for the inclusive rate is 1.37. Similarly, for K (SV) 2 ,the variation is from 1.67 to 0.66, while for the inclusive rate it is 1.35. It shows, particularly, that the shape at higher orders can not be rescaled from lower orders as the differential K-factor varies significantly over the full rapidity range. In the Fig. 5 we plot K factors defined by K

Conclusions
To summarize, we present threshold enhanced N 3 LO QCD correction to rapidity distribution of the Higgs boson produced through bottom quark annihilation at the LHC. We show in detail the infra-red structure of the QCD amplitudes at NLO level as well as the cancellation of the various soft and collinear singularities through the summation of all possible degenerate states and the renormalization of the PDFs in order to demonstrate a general framework to obtain threshold corrections to rapidity distributions to all orders in perturbation theory. We have used factorization properties, along with Sudakov resumma-tion of soft gluons and renormalization group invariance to achieve this. The recent result on three loop form factor by Gehrmann and Kara [57] and the universal soft distribution obtained in [51] provide the last missing information to obtain threshold correction to N 3 LO for the rapidity distribution of Higgs boson in bottom quark annihilation. We find the dominance of the threshold contribution over the entire rapidity range at NLO. We extend this approximation beyond NLO to make predictions for center of mass energies 8 and 13 TeV. We observe that the inclusion of N 3 LO contributions reduces the scale dependency further, as expected, through the variation of the renormalization and factorization scales around their central values and that K-factors show stability at higher orders.