Threshold enhanced cross sections for colorless productions

We study the threshold effect for neutral and charged Drell-Yan productions, associated production of Higgs boson with a massive vector boson and Higgs production in bottom quark annihilation at LHC to the third order in QCD. Using the third order soft-virtual results for these processes and exploiting the universality of the threshold logarithms, we extract the process-dependent coefficients for these processes and resum large threshold logarithms to next-to-next-to-next-to leading logarithmic (N$^3$LL) accuracy. By matching our results to the recently available N$^3$LO results, we provide the most precise theoretical predictions for these processes. We present numerical results for invariant mass distribution and total production cross-sections. We find the conventional scale uncertainties of about $0.4\%$ at N$^3$LO level in the fixed order results get reduced to as small as less than $0.1\%$ at N$^3$LO+N$^3$LL level in the high invariant mass region.

capture the dominant contributions from the large threshold logarithms. On the other hand, the universality of these logarithms allows us to resum them to all orders and capture significant contribution which is essential to compare to the experimental predictions.
The singular part of the fixed order contribution viz. the soft-virtual (SV) corrections to these processes are known for a long time through N 3 LO and beyond [39][40][41][42][43][44][45][46][47][48][49][50]. Some of these results have been used to perform threshold resummation up to N 3 LL. Threshold resummation has been well studied from the early days of QCD for a wide range of processes [51,52]. This is possible due to refactorization [53][54][55][56][57][58][59][60][61][62][63][64][65] of partonic cross sections in terms of soft and hard functions in the threshold region. For color-singlet productions, they are extensively studied to N 3 LL accuracies [22,42,47,[66][67][68][69][70][71][72] and even to N 4 LL [48][49][50] for some processes in the SM. In general the inclusion of threshold resummation results into better perturbative convergence with an improved scale uncertainty. In [69], it was pointed out that the threshold resummation is important in the DY invariant mass distribution in the moderate and high invariant regions and results 2% correction to the fixed order. The scale uncertainty reduces to below 1% in the high invariant mass region Q > 1500 GeV where a matching at third order was limited only to the SV part. A general framework has been provided for arbitrary color singlet production in [73,74]. Beyond the SM (BSM), these are also studied for a wide range of models and stringent bounds have been obtained with models parameters with precise N 3 LL results [75][76][77][78][79][80][81][82][83][84][85]. However, one also needs to match the threshold logarithmic contributions with the FO results to capture the subleading terms which become important at higher orders as well. The recent results on the FO frontiers allow us to study this with the availability of a fixed order code n3loxs [1] shipped with a range of color-singlet processes at N 3 LO.
The purpose of the present article is thus to perform a complete study with properly matching N 3 LL resummed results with the newly available N 3 LO results. We extract all the process-dependent and universal coefficients needed for N 3 LL resummation following the formalism developed in [59][60][61][62][63]. It is worth mentioning that the Higgs boson production in gluon fusion at the LHC has been studied extensively in the literature. The mass of the Higgs boson being around 125 GeV, the underlying parton fluxes at high energy hadron colliders are large. In addition to the fact that the lowest order dominant contribution comes from the gluon fusion channel, the QCD corrections cannot be neglected even beyond NNLO in QCD. This led to the computation of higher order QCD corrections to N 3 LO QCD in the fixed order and to N 3 LL accuracy in the context of resummation, in order to achieve the robust precision studies for this process. The details of these precision studies can be found in [22,39,40,67,75,[86][87][88][89][90][91]. Hence, we will not repeat them here but rather focus on the other color-singlet production processes at hadron colliders. The article is organized as follows: in Sec. II, we briefly lay out the theoretical framework providing the essential formulas to perform threshold resummation and the matching. In Sec. III we present the phenomenological results for different color-singlet processes in the context of LHC, and finally we conclude in Sec. IV.

II. THEORETICAL FRAMEWORK
The hadronic cross section for colorless production at the hadron collider is given by, σ(Q 2 ) = a,b=q,q,g where σ(Q 2 ) ≡ Q 2 dσ/dQ 2 for the DY-type processes and σ(Q 2 ) ≡ σ(M 2 H ) for bbH process. Note that for the total production cross section for V H, we integrate this over the invariant mass Q of the final state V H. The hadronic and partonic threshold variables τ and z are defined as where S andŝ are the hadronic and partonic center of mass energies respectively. τ and z are thus related by τ = x 1 x 2 z. The partonic coefficientσ ab can be further decomposed as follows, σ ab (z, Q 2 , µ 2 F ) = σ (0) (Q 2 ) ∆ sv ab z, µ 2 F + ∆ reg ab z, µ 2 The term ∆ sv ab is known as the soft-virtual (SV) partonic coefficient and captures all the singular terms in the z → 1 limit. The ∆ reg ab term contains regular contributions in the variable z. The overall normalization factor σ (0) depends on the process under study. In particular, the processes under consideration take the following forms, where, Here, V qq are the CKM matrix coefficients with Q q + Q q = ±1 and where Q a is the electric charge and T 3 a is the weak isospin of the fermions. Here, M V and M H are the masses of the weak gauge boson and Higgs boson respectively, m b is the mass of the bottom quark and v being the vacuum expectation value. The function λ which appears in the V H case is defined as The singular part of the partonic coefficient has a universal structure which gets contributions from the underlying hard form factor [92][93][94][95][96], mass factorization kernels [97,98], and soft radiations [62,63,[99][100][101][102]. All of these are infrared divergent which, when regularized and combined give finite contributions. The finite singular part of these has the universal structure in terms of δ(1 − z) and plus-distributions These large distributions can be resummed to all orders in the threshold limit (z → 1). Threshold resummation is conveniently performed in the Mellin (N ) space where the convolution structures become a simple product. The partonic coefficient in the Mellin space is organized as follows: The factor g 0 is independent of the Mellin variable, whereas the threshold enhanced large logarithms are resummed through the exponent G N . The resummed accuracy is determined through the successive terms from the exponent G N which up to N 3 LL takes the form, where N = N exp(γ E ). These coefficients are universal and only depend on the partonic flavors being either quark or gluon. Their explicit form can be found in e.g. [56,59] and are also given in the Appendix A.
In order to achieve complete resummed accuracy one also needs to know the N -independent coefficient g 0 up to sufficient accuracies. In particular, up to N 3 LL, it takes the form, All these coefficients g 0i are given in Appendix A. It is also possible to resum part (or all) of the g 0 by including them in the exponent [22,49,58,67,69], which however have subleading effect as these contributions are not dominated in the threshold region.
To obtain the results in z space, one needs to do the Mellin inversion as While doing this complex integral, one encounters the Landau pole at N = exp 1/(2a S β 0 ) − γ E , and hence the choice of the contour becomes crucial. We choose the value of c in the minimal prescription [53] so that the Landau pole will be on the right of the integration contour and all other singularities will be on the left side of the contour. The Mellin inversion can then be performed [103] along the contour N = c + x exp(iφ), where x is real variable, and c and φ determine the contour. We choose c = 1.9 and φ = 3π/4 to obtain a stable result. The final matched result can be written as, The f a,N are the Mellin transformed PDF similar to the partonic coefficient in Eq. (8) and can be evolved e.g. using the publicly available code QCD-PEGASUS [103]. However, for practical purposes, it can be also approximated by directly using z-space PDF following [52,56,104]. The last term in the bracket denotes that the resummed partonic coefficient in Eq. (8) has been truncated to the fixed order to avoid double counting the singular terms already present in the fixed order through σ N n LO .

III. NUMERICAL RESULTS
In this section, we present numerical results for various color singlet production processes discussed in the previous section in the context of the LHC. Unless specified otherwise, in our numerical analysis we use MMHT2014 [105] parton distribution functions (PDFs) throughout taken from the LHAPDF [106]. The LO and NLO cross sections are obtained by convolving the respective coefficient functions with MMHT2014lo68cl and MMHT2014nlo68cl PDFs, while the NNLO and N 3 LO cross sections are obtained with MMHT2014nnlo68cl PDF sets. In all these cases, the central set (iset=0) is the standard choice. Our default choice of the a S is the same as the one used in the n3loxs code, and it varies order by order in the perturbation theory. The fine structure constant is α 1/132. The unphysical renormalization and factorization scales are chosen to be µ R = µ F = Q, where Q is the invariant mass of the dilepton or the invariant mass of the vector and Higgs bosons in the final state. For the case of Higgs production in bottom annihilation, however, the default scale choice is µ R = m H and µ F = m H /4 following the suggestion from Ref. [18]. For all the processes we have considered here, the scale uncertainties are estimated by varying the unphysical scales in the range so that | ln(µ R /µ F ) | < ln 4. The symmetric scale uncertainty is calculated from the maximum absolute deviation of the cross section from that obtained with the central/default scale choice. To estimate the impact of the higher order corrections from FO and resummation, we define the following ratios of the cross sections which are useful in the experimental analysis:

A. Neutral DY production
For the neutral DY case, the fixed order results and the associated uncertainties have been discussed in Ref. [1] in greater detail. Hence, we will not repeat them here, instead we focus on the resummed results. In the left panel of Fig. 1, we present the invariant mass distribution Q 2 dσ/dQ 2 up to N 3 LO+N 3 LL by varying Q from 250 GeV to 3000 GeV. The corresponding R i0 -factors defined above are given in the right panel. It can be seen that R 20 is larger than R 30 here up to about Q = 2000 GeV, and then they slowly converge to each other, while R 10 being smaller than these two for the entire Q region considered. We also present the R ii -factors which estimate the contribution of higher order threshold logarithms over the respective FO results. The effect of threshold resummation is prominent at NLO. However, its contribution at N 3 LO level is very small. This is expected as the FO results for DY case are converging already at NNLO level onwards, unlike the case of Higgs production in gluon fusion channel [14,86,87,91]. Here, the resummed results at N 3 LO+N 3 LL demonstrate excellent convergence of the perturbation theory.
At this level of precision, it is also important to consider effects from power corrections. As a first step, one can consider the next-to-soft-virtual (NSV) corrections by taking into account the logarithms of the kind ln i (1 − z). In the context of the nDY process, a detailed phenomenology at LHC has been studied in Refs. [71,107]. In Ref. [71], these NSV logarithms have been resummed to next-to-next-to-leading logarithmic accuracy and a comparison of the same has been carried out with the SV resummed result. It is found that the differences between the SV and NSV resummations are about 4% at LL, 2.84% at NLL and 0.85% at NNLL accuracy for Q = 1000 GeV. At higher invariant mass, the differences do not change much, e.g. at Q = 2000 GeV, the differences are about 4%, 2.8% and 0.89% respectively. A general observation is that the difference between SV and NSV resummations is less sensitive to the value of Q and that this difference decreases with higher logarithmic accuracy. If this trend were to continue, at N 3 LL the difference between SV and NSV resummed results could be about 0.1%. This translates to a rough estimate of the order of 0.1% theory uncertainty due to the missing power corrections. It is also worth noting that the NSV logarithms from other partonic channels will also contribute to the cross section. However, a detailed study of such power corrections is beyond the scope of the present work. the different underlying parton fluxes that are different starting from the Born level, the behavior of these R i0 -factors is expected to be different from those of both neutral DY and W − * .
Besides, the resummed results at N 3 LO+N 3 LL play a significant role in reducing the conventional 7-point scale uncertainties. The scale uncertainties in the resummed predictions up to N 3 LO+N 3 LL accuracy are given in Fig. 4 for neutral DY (left panel), for W − * (middle panel) and for W + * (right panel). While the FO scale uncertainties for neutral DY case at Q = 3000 GeV are found to get reduced from about 1.5% at NNLO level to about 0.4% at N 3 LO level, the same in the resummed results are found to get significantly reduced from about 0.2% at NNLO+NNLL level to almost about 0.1% at N 3 LO+N 3 LL level. Thus, DY process is one of the processes for which the theoretical predictions available to-date are the most precise, any uncertainties that can still be present in these cross sections are only due to the PDFs. We will discuss the uncertainty due to PDFs later. It should be noted here that the scale uncertainties in the resummation context will not show improvement over the FO results in the low Q-region, say below 1000 GeV, where the threshold logarithms are not the sole dominant contributions to the cross sections. However, the scale uncertainties in the resummed results will get reduced in the high Q region. To elaborate on this, we compare and contrast the 7-point scale uncertainties in FO and resummed results at third order in QCD (see   To see the effects of resummation over FO results at a given order, it is quite useful to use the factors R ii defined in Eq. (13). We plot in Fig. 6 these factors R 11 , R 22 and R 33 for all three different DY processes as a function of the invariant mass of the dileptons. In all these plots, we can see the R 11 contribution is dominant particularly in the high Q region, while the R 33 is almost unity except for a small contribution in the high invariant mass distribution.

C. V H production
In this section, we present the numerical results for the Higgs production in association with a massive vector boson V = Z, W − , W + . We present results for both the invariant mass distribution and the total production cross sections at hadron colliders for different center of mass energies. However, for invariant mass distribution our results are confined to only DY process i.e. the production of an off-shell gauge boson followed by its decay to on-shell V and H.  It is worth noting that the total production cross sections for these processes up to NNLO are available through the public code vh@nnlo [108]. The updated version vh@nnlo 2.1 [109] of the code can handle both the SM Higgs and other BSM scenarios where Higgs can be produced in association with a gauge boson. In this version, the code is also interfaced with MCFM [110] to produce invariant mass distributions up to NNLO level.
In the present context, we provide the invariant mass distribution up to third order (N 3 LO+N 3 LL). To achieve this we use the invariant mass distributions of DY processes available at N 3 LO level through n3loxs [1] code, and  we numerically incorporate the decay of the off-shell vector boson to V and H instead of to leptons, using the Eq.(2) and Eq.(3) of Ref. [109]. For consistency, we reproduce the results obtained from vh@nnlo 2.1 up to NNLO for V H invariant mass distribution. We then extend our analysis to the fixed order N 3 LO level. For the resummation case, we use our in-house numerical code, similar to the one used for dilepton production case discussed in previous sections. In Fig. 7, we plot the invariant mass distribution of V H at FO (left panel) and at resum level (right panel) up to third order (N 3 LO+N 3 LL) for ZH production process by varying Q from 250 GeV to 3000 GeV. Because the branching of off-shell V * to V H is different from that to dileptons, the production cross sections will certainly be different from that of the neutral DY production of dileptons. However, the corresponding K-factors and R-factors are defined in Eq. (13), expected to be almost same as those of neutral DY case due to the cancellation of the branching part in these ratios.
In Fig. 8, we present these K-factors (left panel) and R-factors (right panel) up to third order (N 3 LO+N 3 LL). The corresponding scale uncertainties for the ZH production case are shown for FO in Fig. 9 and for the resum case in Fig. 10. We notice that the behavior of the scale uncertainties for FO is the same as that for the neutral DY case. For completeness, in Fig. 9, we show separately the uncertainties due to the 7-point scale variations (left panel), those due to only µ R for fixed µ F = Q (middle panel) and those due to only µ F for fixed µ R = Q (right panel) for FO results. We find that the factorization scale uncertainties at higher orders, namely, at NNLO and N 3 LO level are smaller than those due to the renormalization scale. In Fig. 10, we show similar plots as those in Fig. 9 but for resummed results. However, here we find that the factorization scale uncertainties are larger than the uncertainties due to the renormalization scale variations. This is expected because in the resummation case, the factorization scale dependence is included only from the threshold regions and that too in the qq channel. However, at the fixed order at N 3 LO level, this is not the case as the full µ F scale dependence is included in the coefficient functions at this order. We also notice that in general the scale uncertainties in FO results in the high Q region increase but very slowly. On the contrary, the scale uncertainties in the resummed predictions decrease with increasing Q. This is because in the  high Q region the bulk of the cross sections are dominated by threshold logarithms and are resummed to all orders in the perturbation theory.
In Fig. 11 and Fig. 13, we show the results for invariant mass distribution of W − H and W + H production processes respectively. The corresponding K-factors and R i0 -factors are shown in Fig. 12 and Fig. 14. The results for the invariant mass distribution differ from those of the respective dilepton production processes through off-shell W − and W + gauge bosons. However, the ratios K and R-factors are almost the same. Due to the underlying parton fluxes, these K-factors and R-factors for W H production case, on the other hand, will certainly be different from those of ZH case. In Fig. 15, we show the 7-point scale uncertainties in the invariant mass distribution of W − H and W + H processes.
For these new results i.e. the invariant mass distribution for V H process at N 3 LO and N 3 LO+N 3 LL, we compare the 7-point scale uncertainties in FO and in resum results Fig. 16. The behavior of the scale uncertainties is almost identical to that of the neutral DY case (see Fig. 5), namely, the scale uncertainties are smaller in the low Q region for FO case, while the same is true for resum case in the high Q region. There is a very small and negligible difference between ZH case and neutral DY case, which is mostly due to the presence of photon contribution in the latter case.
We also present in Fig. 17, the ratios R ii for i = 1, 2, 3 for ZH, W − H and W + H cases. It is also worth noting that with the automation of most of the NLO calculations [111,112], the NLO results are readily available for these processes and hence it will be particularly useful to estimate the R i1 for i = 1, 2, 3. We plot these R i1 factors for all these V H processes in Fig. 18   For the V H production process, we also give the total production cross sections obtained by integrating the invariant mass distributions over the entire Q region that is kinematically accessible. First, we give these total cross sections for DY type ZH production processes from LO to N 3 LO+N 3 LL for different center of mass energies √ S = 7, 8, 13, 13.6 and 100 TeV in Table I. The corresponding 7-point scale uncertainties are also provided for each case. For the total production cross sections, the bulk of the contribution comes essentially from the low Q region and hence the corresponding scale uncertainties in the resummed results are predominantly due to those coming from the FO results that enter through matching procedure. Hence, the scale uncertainties are smaller for FO case than those for the resum case. However, one can clearly see these scale uncertainties decrease for any center of mass energy as we go from LO to N 3 LO in FO, or as we go from LO+LL to N 3 LO+N 3 LL in the resummation series. For example, for the 13.6 TeV case, the scale uncertainties at LO are 4.06% and get reduced to about 0.33%, while those at LO+LL are 4.44%, and they get reduced to about 0.58% at N 3 LO+N 3 LL. In Table III, we present similar results for W − H case and in Table V for W + H case. In all these results, the general observation is that the scale uncertainties do increase with the center of mass energy √ S of the incoming hadrons. Finally, in the context of V H production, we note that the DY type contributions do not fully give the complete FO predictions for V H case. Starting from O(a 2 S ), particularly the ZH process receives contributions from the gluon fusion channel [113], bottom annihilation channel as well as the heavy top-loop contributions [114]. For the gluon fusion channel, the LO order contribution formally comes at the same order as NNLO of DY type ZH production process and the results for this channel are available already. The higher order NLO corrections to this gluon fusion  channel which contribute at O(a 3 S ) level appear at the same order as the N 3 LO of DY type corrections, and these NLO corrections to the gluon fusion channel in the effective theory [115] have been computed already. It is worth noting that in the recent past, there has been a significant amount of progress toward the computation of the NLO corrections to this gluon fusion process considering the finite top quark mass effects [116][117][118][119][120][121][122]. The top-loop contributions are available at the a 2 S level [114] and are implemented in the vh@nnlo 2.1 code. In the present context, we consider all these three contributions to the ZH production process and define the total production cross section defined as where the power of a S in the parenthesis denotes the highest order up to which the respective contribution has been taken into account. We present these results for the ZH case in Table II. For the case W H production, there will be no contribution from gluon fusion as well as bottom annihilation processes, however, there will be contribution from top-loops from NNLO onwards. Hence, we define the total production cross sections for W H case as These production cross sections up to N 3 LO+N 3 LL are given in Table IV for W − H and in Table VI for W + H cases. We also note that at this level, the electroweak corrections [123,124] are also competitive. For the current LHC energies, they are about −5.28% for ZH and −6.88% for W H total production cross sections. This can easily be implemented in the analysis, however, we did not consider them for simplicity and focus only on the QCD corrections.  Finally, we estimate the uncertainties in our predictions for V H invariant mass distributions due to the choice of the PDFs used in our analysis. For this, we compute the invariant mass distributions at N 3 LO+N 3 LL using different choice of PDFs : ABMP16 [125], CT18 [126], NNPDF23 [127] and PDF4LHC15 [128]. All these sets are considered at NNLO level and the invariant mass distributions are obtained for central set (iset=0). We obtain these results normalized with respect to those obtained from our default choice of MMHT2014nnlo PDFs (iset=0) and for the central scale choice and present them in Fig. 19 for ZH (left panel), W − H (middle panel) and for W + H (right panel). The bands for each of the PDF sets in these figures represent the corresponding 7-point scale uncertainties, which are found to get reduced with Q for all the PDF sets considered here except those for NNPDF23 sets. For the latter choice of PDF sets, the scale uncertainties are found to decrease first with Q up to about 1500 GeV, and then they slowly increase with Q. The uncertainty in these predictions due to the choice of different PDF sets is smaller in the low Q-region and is at the most 4%, but these uncertainties tend to increase with Q. In the high Q region It is worth studying the intrinsic PDF uncertainties in our resummed predictions due to the parametrization of the PDFs themselves. For this, we use the MMHT2014nnlo68cl PDF sets and compute the cross sections due to all the 51 different sets and estimate the asymmetric uncertainties as provided by the LHAPDF routines. We present these uncertainties in our resummed results for V H productin cross sections at N 3 LO+N 3 LL accuracy normalized with respect to those obtained with the central set. These results are shown in Fig. 20 for the ZH, W − H and W + H processes, with uncertainties reaching about 5% in the high invariant mass region.

D. bbH production
As a final process, we consider the Higgs boson production in bottom quark annihilation process at the LHC. This process has already been studied up to NNLO+NNLL by some of the authors in [42] including a detailed uncertainty analysis. In [42], the work also has been extended to N 3 LO+N 3 LL including the estimation of uncertainties only due to unphysical renormalization scale. However, in this work, we complement these studies by considering the bb → H production cross sections along with a detailed 7-point scale uncertainties for the choice of parameters given before. We present in Fig. 21 these production cross sections from LO+LL to N 3 LO+N 3 LL along with the respective scale uncertainties at hadron colliders for different center of mass energy √ S values from 7 TeV to 100 TeV. For √ S = 13.6 TeV, we find that these 7-point scale uncertainties is 5.26% at N 3 LO and is 4.98% at N 3 LO+N 3 LL. The scale uncertainties are found to get reduced with higher logarithmic accuracy in resummation while at any given order in perturbation theory they are in general found to increase with √ S.

IV. CONCLUSIONS
We have studied the threshold effect on the color singlet processes DY, associated V H productions and Higgs production through bottom quark annihilation. To achieve this, we have used the universal and process dependent resum coefficients that are known by some of us up to N 3 LL. Thanks to the recent publicly available n3loxs code that we are able to perform the complete matching at the third order in QCD thus incorporating missing regular contributions at this order. We performed a detailed phenomenology for neutral DY, charged DY production of dileptons, Higgs production in association with massive vector boson (V = Z/W ) in the DY-type process as well as for the Higgs production in bottom quark annihilation process. We presented our results for the dilepton and V H invariant mass distributions as well as the total production cross sections for V H and bbH processes. For the case of V H production process, the resummed corrections over the fixed order N 3 LO results are found to be very small, thus demonstrating a very good convergence of the perturbation series. We observe that the FO results have smaller scale uncertainties in the low Q region. However, the resummation gives significant reduction of the conventional 7-point scale uncertainties to as small as 0.1% in the high invariant mass region Q ≥ 1500 GeV of V H production processes. For the associated Higgs production with a Z boson, the gluon fusion channel and top-loop contributions from NNLO onward are also important [116][117][118][119][120][121][122] at the accuracy that we considered here. We have included these contributions in our inclusive total production cross sections for V H processes. To tame down the uncertainties further, we note that the resummation in gluon fusion channel is equally important. For the DY and DY-type V H production processes we have considered here, the theory uncertainties after the resummation to N 3 LL accuracy has been achieved are very much under control. However, the PDF uncertainties are still large which can be reduced with the availability of N 3 LO level PDFs and more experimental data. Further, at this precision level the electroweak corrections are also important to bring the total theory uncertainties under control. Here N 4 = (n 2 c − 4)/n c and n f v is proportional to the charge weighted sum of quark flavors [95].
The process-independent universal resum exponent defined in Eq. (9) which are used for DY-type processes are given as, Here A i are the universal cusp anomalous dimensions, D i are the threshold noncusp anomalous dimensions, and ω = 2a S β 0 ln N . Note that all the perturbative quantities are expanded in powers of a S . The cusp anomalous