Higgs boson rapidity distribution in bottom annihilation at NNLL and beyond

We present precise resummed predictions for Higgs boson rapidity distribution through bottom quark annihilation at next-to-next-to-leading logarithmic (NNLL) accuracy matched to next-to-next-to-leading order (NNLO) and at next-to-next-to-next-to-leading logarithmic (N3LL) accuracy matched to next-to-next-to-next-to-leading order soft-virtual (N3LOsv) in the strong coupling. Exploiting the universal behavior of soft radiation near the threshold, we determine the analytic expressions for the process-dependent and universal perturbative ingredients for threshold resummation in double singular limits of partonic threshold variables $z_1,z_2$. Subsequently, the threshold resummation is performed in the double Mellin space within the standard QCD framework. The new third-order process-dependent non-logarithmic coefficients are determined using three-loop bottom quark form factor and third-order quark soft distribution function in rapidity distribution. The effect of these new resummed coefficients are studied at the $13$ TeV LHC. We observe a better perturbative convergence in the resummed predictions on the Higgs rapidity spectrum in bottom quark annihilation. We also find that the NNLL and N3LL corrections are sizeable which typically are of the order of $-2.5\%$ and $-1.5\%$ over the respective available fixed orders with the scale uncertainty remaining at the same level as the fixed order.


Introduction
The Standard Model (SM) Higgs boson is one of the important fundamental particles to study at colliders like the Large Hadron Collider (LHC).Understanding the Higgs properties is crucial to know the SM well and is also critical to the search of new physics beyond the SM (BSM) where new physics might couple to the Higgs sector.Testing the Higgs properties and understanding its interactions with other fundamental particles are indeed main tasks at the LHC in the upcoming runs.Precision calculations play a prominent role in these studies by calculating the higher-order contributions in the perturbation theory and improving the predictions of Higgs boson properties to a very high accuracy.
The dominant mode of Higgs production at the LHC is through gluon fusion.On the other hand, the Higgs production in bottom annihilation channel, despite being subdominant is also interesting to study.Firstly, Higgs dominantly decays to bottom quarks which can give direct access to the Higgs Yukawa coupling.This purely hadronic final state, however, is challenging to measure [1,2].This also requires production of Higgs boson at the first place, through bottom annihilation channel along with the dominant gluon fusion.Secondly, it also gives access to the Higgs Yukawa coupling to bottom quarks even when Higgs is decayed through cleaner channels like di-photons [3,4] or four-leptons [5,6] productions.Although in the SM the Higgs Yukawa coupling to bottom quark is suppressed by small bottom mass, in the extensions of SM e.g.two-Higgs doublet models, or minimal supersymmetric standard model the coupling could be enhanced and is crucial to the search for new physics.Therefore, a precise understanding of the SM contribution will be beneficial in BSM analyses.Thirdly, the Higgs production through bottom annihilation is also interesting on how bottom quark is treated: as whether it is a part of the proton, taking bottom quark as a massless parton except in the Yukawa coupling which is done in the 5 flavor scheme (5FS), or whether it is taken as a massive quark throughout and excluded from the proton structure as is done in 4 flavor scheme (4FS).While in 4FS a massive bottom quark is produced from gluon splitting from proton, in the 5FS scheme, massless bottom has its own parton distributions.
Due to its high importance, Higgs inclusive cross-section is now available theoretically to a very high accuracy to next-to-next-to-next-to leading order (N3LO) [7,8] in Higgs effective field theory (HEFT) in the gluon fusion channel providing a correction of around 2% with scale uncertainty around 3% reducing from 9% at NNLO [9][10][11].It is also known up to N3LO in the vector boson fusion [12][13][14][15] where the correction already stabilizes, and the scale uncertainty is found to be below 0.2% at N3LO.Beyond the fixed order, efforts were made to study the dominant threshold contributions by studying the soft-virtual (SV) corrections at fourth order as well as partial subleading logarithmic effects in gluon fusion [16,17] with a further enhancement of the cross-section ranging from 0.2%-2.7%depending on the scale choices.The fixed order cross-section is further improved by performing threshold resummation at the next-to-next-to-leading logarithmic (NNLL) accuracy [18][19][20][21] and to the third logarithmic accuracy (N3LL) [22,23] as well as resummation of π 2 terms [24] arising from time-like Sudakov form factor.The finite top mass effect is also known at NLO [25,26], partially at NNLO with top mass expansion [27,28,[28][29][30] and recently to exact NNLO [31] where an increment of 0.6% is observed compared to the HEFT approximation.The electroweak corrections are also known to NLO [32,33] which amount to a positive correction of about 5% compared to the NNLO QCD.Higgs production through bottom annihilation, despite being the subdominant channel, has also been studied extensively in the literature.Due to the availability of third order form factor [34], the soft distributions [35], and the relevant splitting functions [36,37], the inclusive cross-section is known up to N3LO [38][39][40][41][42] in the 5FS where the residual scale uncertainties are found to be around 5% and to NLO [43][44][45] in the 4FS.There have been several studies to combine the 5FS and 4FS prediction through different matching prescriptions [46][47][48][49][50][51][52][53][54].Further a complete N3LL resummation is performed in [55,56] where the scale uncertainty at this order reduces to about 4.9%.The pure QED and mixed QCD-QED effects are studied in [57] where the corrections are found to be below 0.03% of the LO.
Differential measurement like the rapidity distribution of the Higgs boson is important to understand Higgs interaction within the SM.Similar to the total production crosssection, rapidity is also inclusive to extra radiations.While this sheds light on the spin of the particle itself, it is also useful to constrain the parton distribution functions (PDFs).In particular, the region with large momentum fraction is not well-constrained where the resummed results could play an important role.Rapidity distribution has been known to NNLO for the Higgs production through gluon fusion [58,59], as well as in bottom quark annihilation [60,61].The accuracy is further extended beyond the NNLO level by studying the threshold contributions up to the third order [62,63] in gluon fusion as well in to bottom annihilation [64].It was observed there that the threshold contribu-tions play a prominent role in the rapidity distributions for these processes.Recently even complete N3LO corrections are also obtained for Higgs [65] and Drell-Yan (DY) [66] rapidities using q T subtractions [67][68][69][70] where the corrections are shown to be around 3% and −2% respectively over NNLO in the central rapidity range.However, unlike the DY case, the uncertainty band for the Higgs rapidity distribution at N3LO show a nice perturbative convergence, reducing the scale uncertainty below 5% and residing within the NNLO uncertainty band.Beyond NNLO, large threshold logarithms are also included at next-tonext-to-leading logarithmic (NNLL) accuracy [71] in the gluon fusion channel.Efforts [72] are also made to resum partial subleading logarithms in the gluon fusion.
In this article we focus on the bottom induced rapidity distribution which is known to NNLO [60,61] for quite some time.We aim to improve this by including the threshold effects at NNLL and beyond.Within the traditional QCD resummation framework, a formalism [71,[73][74][75][76][77] has been already developed to resum the large threshold logarithms in rapidity distribution for colorless particles.Originally the formalism was proposed for the x F distribution in the seminal work [73] by Catani and Trentadue.Later it was extended for rapidity1 following the framework developed in [62,78].The idea is to identify proper scaling variables z 1 , z 2 corresponding to partonic threshold (z) and rapidity (y p ).One then resums large rapidity logarithms by resumming these scaling variables simultaneously going to the threshold limit z 1 , z 2 → 1.This was termed [71,75,76] as the Mellin-Mellin (M-M) approach as the resummation is performed in the double Mellin space corresponding to z 1 , z 2 .Essentially in this approach one resums all the double singular terms arising from the delta function δ(z i ) and plus-distributions ln n zi zi + where zi = 1 − z i .This is also consistent with the generalized threshold resummation approach [79] employing softcollinear effective theory (SCET) in the double singular limit.A recent comparison for different approaches up to the next-to-leading power (beyond double soft) level can be found in [80].
In this article we follow the standard M-M approach at the leading power within the traditional QCD resummation framework and study the impact of these threshold logarithms for the Higgs boson rapidity distribution in the bottom quark annihilation channel.We organize the paper as follows: in sec.(2) we lay out the theoretical framework for the resummation of rapidity distributions in the M-M approach, in sec.(3) we present the results relevant at the 13 TeV LHC, and finally we conclude in sec.(4) and collect all the analytical results required upto N3LL in the appendices (A -D).

Theoretical Framework
The effective Lagrangian for the interaction of a scalar Higgs boson with the bottom quark is given as, where ψ b (x) and ϕ(x) are the bottom quark field and the Higgs field respectively.Here λ is the Yukawa interaction which is given by λ = m b /v, with v being the vacuum expectation value (VEV), and m b being mass of the bottom quark.We follow the 5FS where we use non-zero mass of the bottom quark only in the Yukawa coupling, elsewhere it is treated as a massless quark.The rapidity distribution of the Higgs boson at proton-proton collider takes the following form, dσ(τ, y) dy = i,j=q,q,g Here τ = m 2 H /S = x 0 1 x 0 2 with S being the hadronic center of mass energy.The hadronic rapidity is defined as y = 1 2 ln x 0 1 /x 0 2 .The rapidity-dependent partonic coefficient function ( σ d,ij ) can be decomposed in terms of singular SV piece consisting of plus-distributions and delta function in partonic threshold variables (z 1 , z 2 ) and non-singular or regular piece as, 3) The singular SV (∆ sv d,ij ) part gets contributions only from the diagonal channel i.e. in the present case i, j = b, b for the SV part [81].On the other hand, the regular terms (∆ reg d,ij ) are subdominant at the threshold and gets contributions from all partonic channels.The overall Born normalization factor σ 0 (µ r ) takes the following form, (2.4) Note that the µ r dependence in the Born factor above comes only through the running of bottom mass (m b (µ r )) or equivalently the Yukawa (λ(µ r )).The Yukawa running is performed using the mass anomalous dimensions γ m which we collect in the appendix (B) up to the third order.
Resummation is conveniently performed in the Mellin space where the double Mellin transformation is performed taking Mellin transformation on both partonic threshold variables z 1 , z 2 in the following way (suppressing the µ r , µ f dependence), Here N i is the Mellin variable corresponding to partonic threshold z i .The singular SV contribution can be resummed through the integral form in terms of universal cusp anomalous dimensions A b and rapidity-dependent threshold non-cusp anomalous dimension D b d as well as process-dependent coefficients g ′ b 0 .In the double-Mellin space this can be written in the following integral form, where z 12 = M2 H z 1 z 2 and zi = 1 − z i .Expansion of the above expression at a fixed order in strong coupling will reproduce the Mellin space fixed order results corresponding to eq. (2.5).Performing the Mellin integration above will produce some non-logarithmic constants in double Mellin space which can be combined with the process-dependent prefactor (g ′ b d,0 ) and one can finally organize the large threshold logarithms in the Mellin space in the following form, The non-logarithmic coefficient (g b d,0 ) has the following perturbative expansion To achieve N3LL accuracy, they are needed up to the third order in strong coupling.
We have obtained these using the third order b bH form factor [34] and third order soft distribution function [35], and we present these in appendix (C).The exponent is universal and resums the large logarithms (which now appear in N i → ∞ limit) to all orders.It can be expanded in the strong coupling and the inclusion of successive terms defines the resummed order, 2 where ω = a s β 0 ln(N 1 N 2 ), with N i = e γ E N i , i = 1, 2. These process-independent resummed exponent are same as the quark-initiated Drell-Yan process and up to N3LL accuracy these can be found e.g. in [82].The resummed expression in eq.(2.7)only resums the leading singular terms which appear through the large logarithms of N 1 and N 2 in Mellin space.This lacks the subleading regular pieces which can be included through the available fixed order results to improve the resummed predictions.However, one can not simply add them as the fixed order also contains the same logarithmic contributions up to a certain order which are already taken into account in the resummed expression.Therefore, a matching procedure has to be invoked removing these logarithms which also appear in the fixed order.This is done through the following all order matched expression, The first term on the right-hand side of the equality is the fixed order contribution containing singular and regular contribution up to a fixed order in strong coupling.The first term inside the square bracket organizes the resummed series up to a certain logarithmic accuracy provided by the knowledge of cusp and rapidity anomalous dimensions as well as the process-dependent g b d,0 coefficients.The symbol 'f.o.' in the second term inside the square bracket means the function is truncated to a fixed order in order to avoid double counting of singular terms already present in the fixed order ( dσ f.o.dy ).Here ) are the PDFs in the Mellin space.In practice, we use the x-space PDF through the LHAPDF6 [83] interface using the derivatives of PDF as described in [84].
The Mellin inversion in eq.(2.10) is not straightforward as the resummed expression diverges when ω = 1.This corresponds to the Landau pole where strong coupling diverges.The perturbative formalism thus break down in this region.One way to proceed is by choosing the contour of the Mellin inversion according to the minimal prescription (MP) [85].The basic idea is to choose the contour in such a way that all the poles remain at the left of the contour except for the Landau pole which remains far right of the contour.In double Mellin space this is little involved as the Landau pole is now a function of two Mellin variables.Typically, one needs to project the complex integration on real variables r i and chose the contour accordingly.One can still fix the contour of one of the Mellin variable (N 1 = c 1 + r 1 exp(iϕ 1 )) according to MP.Once one fixes the contour for N 1 through c 1 and ϕ 1 , the Landau pole is not constrained anymore on the real axis of the second Mellin variable (N 2 = c 2 + r 2 exp(iϕ 2 )) and in fact, it now depends on the first Mellin variable (N 1 ).In order to satisfy the MP, a reasonable choice [74,82] of the contour for the second Mellin variable as

Numerical Results
With the setup introduced in the previous section, we now focus on the numerical impact of the threshold logarithms in the rapidity distribution.We focus on 13 TeV LHC with CT14 as our default PDF choice which is used through the LHAPDF6 [83] interface. 3The fixed order results are obtained from [61] using the N-jettiness slicing method [86,87] as implemented in MCFM [88][89][90].The strong coupling is taken through the corresponding PDF sets.At the third order we evolve the strong coupling using four-loop QCD beta function [91][92][93][94][95][96][97][98][99][100][101][102] for which we set the initial condition as α S (m Z ) = 0.118 where m Z = 91.1876GeV is the Z boson mass.We set the Higgs mass to be m H = 125 GeV.The central scale choice for this process is taken as (µ c r , µ c f ) = (1, 1/4)m H GeV. The choice of the low µ c f scale is done following the observation in [40] to minimize the effect of large collinear logarithms which appear at µ c f ∼ m H /4. The Yukawa coupling is also evolved through renormalisation group equation using 4-loop mass anomalous dimensions with bottom mass m b (m b ) = 4.18 GeV.The required mass anomalous dimensions are collected in the appendix (B).To have an estimation of the residual scale uncertainties we follow the standard seven-point scale variations around the central scale choice stated above, with the restriction 1/2 ≤ (µ r /µ c r )/(µ f /µ c f ) ≤ 2. This amounts to seven configurations for the scale (µ r , µ f ).For each bin, the uncertainty envelope is obtained by considering the maximum and minimum deviations from the central scale choice.Note that the scale variation also includes the scale dependence as arising from the Yukawa running in the MS scheme.The double Mellin inversion in eq.(2.10) is performed with an in house code which we also interface to LHAPDF6 as well as to Cuba [103,104] for the final integration.Accordingly we chose the contour in eq.(2.10) as c 1 = c 2 = 1.9 and ϕ 1 = 3π/4 and ϕ 2 given in the previous section.
We further define the following perturbative quantities in order to assess the higher order effects.We define the ratios K-factor and R-factor [105][106][107] corresponding to the fixed order and resummed order respectively as, Further we define the RF ij rations to estimate the resummed contributions over the fixed orders as, For such ratios, we calculate the correlated error by taking both the numerator and the denominator at the same scale and obtain the seven-point uncertainties from seven such ratios.In order to shorten the notation, we denote +LL, +NLL, +NNLL to indicate the complete matched resummed results, meaning the resummed corrections are matched to the FO results according to eq. (2.10).At the third order similar notation is adopted where the complete resummed N3LL resuls are matched to the N3LOsv results and are denoted as +N3LLsv.
In fig.
(1), we present the rapidity distribution of the Higgs boson in bottom quark annihilation up to NNLO (left panel) and to +NNLL (right panel).The asymmetric band is obtained by taking the envelope of maximum(minimum) deviation from the central scale according to the seven-point scale variation.On the fixed order side, the NLO gets a correction of similar size to LO, which gets further increased at NNLO by 9.3% compared to NLO in the central rapidity region (at y = 0 − 1.6).In the higher rapidity region (at y = 2−4), the behavior is different where NLO gets negative correction up to 50% compared to LO, whereas on the other hand, the NNLO gets positive correction up to 27% of NLO.The corresponding scale uncertainties are respectively +62.9 −53.9 % at LO, +22.0 −31.1 % at NLO, and −20.0 % at NNLO in the central rapidity region.In the resummed case, the convergence is faster with +NLL getting a correction 73% compared to +LL whereas +NNLL gets a further increment of 10.1% compared to +NLL in the central region.Corresponding corrections at the higher rapidity (y = 3.2) region are −31.7% and 8.6% respectively.The scale uncertainty does not improve compared to the FO, a behavior which is also seen in the neutral and charged DY productions [82].In the central region, the asymmetric scale uncertainties are +63.at +NNLL respectively.We notice that while +LL gets a positive contribution ranging from 13.7% to 22.8% over LO, the +NLL gets −3.2% to −2.3% corrections over NLO and +NNLL gets relatively flat correction ranging −2.4% to −3.1% going from central rapidity to higher rapidity.On the right panel of fig.
(1), we present the new third order SV results which is further matched with the third order resummed results.While SV results at the third order gets a increment up to 0.5% compared to the NNLO, the scale uncertainty is not improved.The third order matched results however gets a flat correction of about 1.8% throughout the rapidity region compared to +NNLL with the scale uncertainties are at the similar level as the fixed order.It is also possible to match the resummed result in a multiplicative way instead of the additive matching in eq.(2.10).In order to estimate the effect of such a procedure, we followed the prescription4 as presented in [108].We find that the +NLL gets an increment of around 84% at y = 0 compared to +LL whereas at +NNLL the correction is less that 1% of +NLL.Similarly, at y = 3.2, the corrections are −31% and 5% respectively.Thus, we observe a faster convergence for multiplicative matching compared to the additive matching.The corresponding asymmetric scale uncertainties at y = 3.2 are +53.9−50.8 % at +LL, +25.0 −33.8 % at +NLL, +6.4 −2.9 % at +NNLL respectively.Up to +NLL, the uncertainties remains similar compared to the additive matching, at +NNLL, the uncertainties gets better at higher rapidities.On the other hand we see similar scale uncertainties as the additive case at y = 0 at +NNLL level.In the rest of the article we follow the additive matching as provided by eq.(2.10).A better way to visualize the higher order effects is through the ratios viz. the K, R, and RF factors as defined in eq.(3.1)-eq.(3.2) which are presented in fig.
(2).In the left panel we show these ratios up to the second order.It is clear that NLO (or +NLL) corrections shapes the rapidity distributions very well, whereas the corrections from NNLO (or +NNLL) are rather flat over a large rapidity region.Compared to the K 21 factor we observe slight increment of the central scale in the case of R 21 in the higher rapidity region.On the right panel of fig.
(3), we present the these ratios at the third order.Again we observe a relatively flat QCD correction over the large rapidity range at N3LOsv amounting to about 8% uncertainty at the higher rapidities.The matched result becomes almost flat even in the higher rapidity region and the correlated scale uncertainty reduces to 4%.
To estimate the intrinsic PDF uncertainty, we define the quantity δPDF = 1 ± δ[ dσ dy ]/[ dσ dy ] 0 × 100%, where the numerator δ[ dσ dy ] is the intrinsic PDF uncertainty and the denominator [ dσ dy ] 0 is the central prediction.We study only the SV part of the fixed order which might be improved by resummed results.This is shown in fig.(4).Notice that the central predictions are different for FO and resummed cases.Indeed, we observe a −2.5% to −2.1% change in the central predictions for resummed case compared to FO SV while going from y = 0 to y = 4.The PDF uncertainties at the second order (SV) results are about 8% in the central rapidity (y = 0) which gets reduced to about 4.5% at y = −2.At the higher rapidity region the uncertainties further increases.This is a known behavior which is also observed in the case of DY rapidities [82].The matched resummed results show a slight improvement ( below 0.1%) over the fixed order PDF uncertainty in all rapidity region.

Conclusions
At the LHC, soft gluons play an important role in predicting observables correctly in different phase space corners.Particularly in the threshold region their contribution dominates and hence to have a reliable prediction one needs to resum them and match them to the available fixed order to better predict an observable.For rapidity distribution in the threshold region both threshold variables corresponding to partonic threshold and rapidity become equally important and one needs to resum both of them consistently order by order in the perturbation theory.In QCD, this is achieved by the M-M approach, and we exploit this to resum large threshold logarithms to NNLL accuracy matched to available NNLO result.The threshold effects at fixed order amount to −3.2% enhancement over the NLO distribution over a large range of rapidity.On the other hand, the resummed corrections at NNLO+NNLL amount to −2.5% enhancement over the fixed order NNLO.
In general, we find a better perturbative convergence in the resummed spectrum which combines merit of both resummed logarithms and non-singular contributions from fixed order which are not captured in purely resummed predictions.While the third order analytic ingredients presented in this article will be useful to match to the third order fixed order results once it will be available, we present the first predictions at the third order coming from the dominant threshold logarithms and also match them to N3LL.We observe a relatively flat correction in the third order matched results for all rapidities with reduced scale uncertainties to below 4% particularly in the higher rapidities, a region dominated by the large threshold logarithms.Our resummed results can be useful in constraining the bottom quark PDFs at the large momentum fraction.
G.D. is grateful to C. Williams for providing the NNLO results.G.D. also thanks V.
Ravindran for encouraging to work on this project.Part of the analytical computation has been performed using symbolic manipulation system Form [109,110].Simulations were performed with computing resources granted by RWTH Aachen University under project rwth1298.The research of G.D. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 -TRR 257 (Particle Physics Phenomenology after Higgs discovery.).

A Anomalous Dimensions
The cusp anomalous dimensions have the following perturbative expansion in strong coupling, We collect here the coefficients up to fourth order [111][112][113] needed at N3LL, The quartic casimirs appearing above are defined as The coefficients are the same as the quark ones (see e.g.[82]).Up to the N3LL accuracy, these are needed up to third order [62] which we collect below,

B Yukawa running
The Higgs Yukawa coupling with bottom quark is given as λ = m b /v.Here m b (µ r ) is the MS running mass of the bottom quark.Thus, the running of Yukawa goes through the running of bottom mass as, The mass anomalous dimension (γ m ) has the following perturbative expansion: where the coefficients are known up to four loops [99,[114][115][116][117].We collect all the coefficients up to four-loop order, Below we present the new process dependent coefficients (g b d,0i ) up to N3LL accuracy,

D Soft-virtual coefficients in double Mellin space
For completeness, here we collect all the singular SV coefficients [81] up to third order in the double Mellin space as defined in eq.(2.5).The perturbative expansion in Mellin space takes the following form, Defining L ≡ ln(N 1 N 2 ), the coefficients up to third order take the form, ∆

Figure 1 .
Figure 1.Higgs rapidity distribution in bottom annihilation for 13 TeV LHC.The left figure compares different fixed order and matched resummed orders.The notation +NkLL indicates that the matched resummed results are obtained by matching with the corresponding FO results.The right figure compares the SV results at the third order against the matched result at the same order where the matching is done with corresponding SV results.

Figure 2 .
Figure 2. RF factor (as defined in eq.(3.2)) along with correlated errors up to the third order.For the third order, resummed N3LL results are matched to the N3LOsv results.

Figure 3 .
Figure 3. K, R factors (as defined in eq.(3.1)) along with correlated errors upto the second order (left) and the same at the third order (right).The third order results are up to SV accuracy and same for matched.

Figure 4 .
Figure 4. Effects of matched NNLL results to NNLOsv are presented and compared against the NNLOsv for CT14 PDF at 13 TeV LHC.