Resummed Higgs boson cross section at next-to SV to $ \rm NNLO + \rm \overline {NNLL}$

We present the resummed predictions for inclusive cross section for the production of Higgs boson at next-to-next-to leading logarithmic ($\rm \overline {NNLL}$) accuracy taking into account both soft-virtual ($\rm SV$) and next-to SV ($\rm NSV$) threshold logarithms. We derive the $N$-dependent coefficients and the $N$-independent constants in Mellin-$N$ space for our study. Using the minimal prescription we perform the inverse Mellin transformation and match it with the corresponding fixed order results. We report in detail the numerical impact of $N$-independent part of resummed result and explore the ambiguity involved in exponentiating them. By studying the K factors at different logarithmic accuracy, we find that the perturbative expansion shows better convergence improving the reliability of the prediction at $\rm NNLO + \overline{NNLL}$ accuracy. For instance, the cross-section at $\rm NNLO + \overline{NNLL}$ accuracy reduces by $3.15\%$ as compared to the $\rm NNLO$ result for the central scale $\mu_R = \mu_F = m_H/2$ at 13 TeV LHC. We also observe that the resummed $\rm SV + NSV$ result improves the renormalisation scale uncertainty at every order in perturbation theory. The uncertainty from the renormalisation scale $\mu_R$ ranges between $(+8.85\% ,-10.12\%)$ at $\rm NNLO$ whereas it goes down to $(+6.54\% , - 8.32\%)$ at $\rm NNLO + \overline{NNLL}$ accuracy. However, the factorisation scale uncertainty is worsened by the inclusion of these NSV logarithms hinting the importance of resummation beyond $\rm NSV$ terms. We also present our predictions for $\rm SV + NSV$ resummed result at different collider energies.


Introduction
The milestone discovery of the Higgs boson at the Large Hadron Collider (LHC) [1,2] puts the Standard Model (SM) on a firm position to describe the dynamics of all known elementary particles. There have been tremendous efforts in constructing models that address the shortcomings of SM and, at the same time, search for new physics signatures. These studies demonstrate rich phenomenology that can be explored at present and future colliders. In addition, model-independent studies using the sophisticated effective field theory approach have also aided to a better understanding of the underlying theory. To improve the search strategies of new physics signatures at the LHC and thereby to arrive at sensible conclusions, we require precise predictions from SM as well as beyond SM (BSM) scenarios with minimum theoretical uncertainties. Such studies within the SM by both theorists and experimentalists have set stringent constraints on the parameters present in various BSM scenarios [3]. At the LHC, Higgs bosons are produced dominantly in the gluon fusion channel, and hence most of the theoretical studies centers around this. In addition, the radiative corrections from perturbative Quantum Chromodynamics (QCD) beyond the leading order (LO) were found to be sizeable. Hence efforts to include higher order effects became inevitable. As the accuracy in Higgs boson measurements improved, the level of precision in the predictions was also increased by including contributions from other partonic channels, the radiative effects from the electroweak (EW) sector of the SM, and other tiny effects resulting from finite quark masses.
In hadronic collisions, the Higgs bosons are dominantly produced in the gluon fusion subprocess through a top quark loop. Beyond next to LO (NLO), QCD corrections were obtained in the Higgs effective field theory (HEFT) where top quark degrees of freedom are integrated out [4][5][6] and this approach turned out to be the most advantageous technically. By multiplying higher order corrections computed in HEFT with the exact leading order cross section, one can retrieve most of the top-quark mass effects. See [5,[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] for a complete list of results for the Higgs boson production in gluon fusion both in exact and in HEFT up to NNLO level and also for threshold effects beyond NNLO in QCD. The first step towards complete N 3 LO QCD corrections was obtained by expanding the partonic contributions near the threshold [24]. The electroweak corrections to this channel were computed in [26][27][28][29][30][31]. A complete study taking into account all these effects was done in [32][33][34]. In [35,36] the full z-dependence for all the partonic channels was recovered at N 3 LO. To estimate the error resulting from finite top quark mass in HEFT, the top quark mass effects were also systematically studied. In particular, they were included at large partonic center of mass energy as well as from 1/m t terms [37][38][39][40] and they amount to 1% correction. Recently, in [41], the exact top quark effects were included at NNLO for all the channels and this results in a decrease in cross section compared to HEFT to −0.26% at 13 TeV and −0.1% at 8 TeV.
While these fixed order predictions often provide reliable estimates of the radiative corrections, there are certain kinematical regions in phase space, where these estimates can go wrong due to the presence of very large logarithms at every order in perturbative series. In the inclusive reactions, this happens in the threshold regions where most of the partonic center of mass energy goes in producing the Higgs boson and at the same time the Parton distribution functions (PDFs) that are driving these partonic reactions dominates. These large logarithms can be systematically resummed using standard threshold resummation [42,43]. For the Higgs boson production, in [44] the results up to N 3 LL accuracy [45][46][47] were used to predict matched results at N 3 LO+N 3 LL for µ R = µ F = m H which resulted in 4% error from the missing higher effects at 13 TeV LHC.
The successive terms in the threshold expansion of various partonic channels are found to be numerically important and in many cases, or they are even larger than the leading contributions. In particular, the subleading threshold logarithms, called the next-to soft virtual (NSV) terms in the inclusive production of Higgs boson dominate over the leading ones as pointed out in [48]. There have already been several studies [20,25,[49][50][51][52][53][54][55][56][57][58][59] to understand the structure of NSV terms. Recently, in [60], we set up a formalism that can resum these NSV logarithms in Mellin space and in [61], we studied their numerical significance for the production of a pair of leptons in Drell-Yan process. In this article, we extend the study to the production of Higgs boson in gluon fusion at the LHC.
The plan of the paper is as follows. Sec.2 is devoted to provide a theoretical overview of our formalism for completeness. Here, we first describe the importance of threshold expansion as well as the soft gluon exponentiation and then review the formalism given in [60] to include NSV terms for the production of Higgs boson at the LHC. We also discuss different resummation prescriptions to explore the ambiguity involved in exponentiating the N-independent constants. In Sec. 3, we study the numerical impact of our NSV resummed result in detail and present our findings. Finally we conclude in Sec.4.

Theoretical overview
In QCD improved parton model, the inclusive cross section for producing Higgs boson in the hadronic collisions can be written as a convolution of perturbatively calculable coefficient functions (CFs), ∆ ab , and non-perturbative fluxΦ ab : . We work in the HEFT framework, where the CFs are multiplied with the factor G B which is the product of born cross section computed in full theory with finite quark masses and the square of Wilson coefficient in HEFT. The factor G B explicitly depends on both m H and m t and is renormalised at the renormalisation scale µ R (see [15] for the details). The flux, renormalised at the factorisation scale µ F , is the convolution of parton distribution functions f a , f b of incoming partons a, b respectively and is given bỹ The CFs in HEFT are perturbatively calculable in powers of a s = g 2 s /16π 2 where g s is the QCD strong coupling constant: with the normalisation ∆ (0) ab = δ(1 − z). Beyond leading order, the general structure of CFs can be decomposed into, ∆ ab (µ 2 R , µ 2 F , z) = δ bā ∆ SV aā (µ 2 R , µ 2 F , z) + ∆ reg ab (µ 2 R , µ 2 F , z) , (2.4) where ∆ SV aa denotes the soft-virtual (SV) corrections while ∆ reg ab constitutes to the regular part of CFs. Here, the SV corrections resulting from pure virtual contributions in the born process g + g → H and also the leading contributions from gluon gluon initiated partonic channels with at least one emission of on-shell parton. While the pure virtual contributions are proportional to only δ(1 − z), the latter ones contain additionally the distributions of the form D k (z): (2.5) These distributions satisfy for any regular function f (z). Further decomposing the SV CFs order by order, we get with ∆ (i) δ refers to the coefficients proportional to δ(1 − z) and ∆ (i) The regular part of CFs, denoted by ∆ reg ab , contain functions of the form (1−z) m ln k (1− z), m, k = 0, 1, · · · ∞. It can be seen that, the computation of the SV part of the CFs are relatively easy compared to the regular part. For instance, the Feynman integrals that contribute to the pure virtual corrections depend only on the hard scale m H and those corresponds to the remaining soft part of SV contributions proportional to the soft scale m H (1 − z). However, the regular part depends on more than one scale through various powers of m H (1 − z) α , with α ≥ 1 leading to their rich analytical structure. Already at NNLO level, the Feynman integrals in the perturbative expansions contain rich analytical structure through the appearance of iterated integrals [62], which are often called generalised polylogarithms. Further the N 3 LO [35,36] computation demonstrates the appearance of a new class of functions called elliptic functions (see for example [63]) in the iterated integrals for the first time.

Threshold expansion and soft gluon exponentiation
Owing to the complexity involved in computing CFs exactly in the variable z at higher orders, one resorts to the method of expansion near threshold region, namely around z = 1.
Note that the scaling variable z quantifies the fraction of partonic center of mass energy that goes into producing on-shell Higgs boson in hadronic collisions. In the threshold region, the partons that accompany the Higgs boson are mostly soft and the Feynman integrals from the real emission subprocesses in this region can be evaluated as a power series in 1 − z. These expansions show remarkably simple structure in terms of logarithms of 1 − z: The above threshold expansion of CFs in partonic scaling variable z will be useful provided the partonic fluxΦ ab (τ /z) that multiplies them to give hadronic cross sections also dominates in the same region for a given hadronic scaling variable τ . It is interesting to note that the first NNLO result [13] for the inclusive Higgs boson production was achieved through the expansion of CFs around z = 1. Similarly, at N 3 LO level, in [48], next-to soft contributions from all the partonic channels were obtained for the Higgs production. Later on, in [24], beyond next-to soft terms at N 3 LO were obtained using threshold expansion.
In [35], the analytical results for all the channels were reported for the first time and numerically they were found to be consistent with those obtained using threshold expansion [24]. While the state-of-the-art predictions at N 3 LO level is good enough for most of the phenomenological analysis, the results in perturbation theory demonstrates rich analytical structure which can be exploited to study their all-order behaviour. One such study was pioneered for the inclusive processes by Sterman [42] and Catani and Trentedue [43] who showed that the leading threshold contributions can be systematically resummed to all orders in a process independent way. This goes by the name of threshold resummation and the resummed results are organised in terms of certain exponents that are controlled by universal anomalous dimensions. It is convenient to perform resummation in Mellin N space and the leading threshold logarithms in z space translate to large ln N terms in N space. These large ln N s accompanied by small a s at every order give order one terms, ω = 2a s (µ 2 R )β 0 ln N , spoiling the reliability of the perturbative results. The resummation of these order one terms provide sensible expansion. The successive terms in the expansion are arranged as leading logarithm (LL), next-to LL (NLL) and so on. For the Higgs boson production, the results up to N 3 LL level is available [45][46][47]64] and in [44], analysis at N 3 LO+N 3 LL at µ R = µ F = m H provided the estimate of 4% error from the missing higher order effects at 13 TeV LHC. In addition to standard QCD approach, the threshold effects in the large order perturbation theory can be studied using soft collinear effective theory (SCET) which captures the soft and collinear modes at the Lagrangian level. Using SCET [65][66][67] one can conveniently compute the SV results order by order in perturbation theory. The intrinsic scales that separates various modes in theory can be used to set up renormalisation group equations [68] whose solutions sum up large logarithms from threshold regions to all orders. The numerical impact of these results provide reliable estimates of missing higher order effects and they are found to be very close to standard QCD ones.

Going beyond soft gluon exponentiation
Inclusive results that are available to third order for DY and Higgs productions suggest that the subleading logarithms present in the regular part of CFs, ∆ reg ab from various partonic channels also play important role in the precise predictions. In particular, the leading logarithmic terms in ∆ reg ab (z) near z = 1, namely where ∆ N SV,(i) ab (z) is defined by setting l = 0 in (2.8),i.e., were shown to spoil the approximations based on threshold expansions as the order of perturbation increases [48]. One finds that these NSV contributions are larger than SV ones at the LHC energies. Hence, these corrections as well as those from O(1 − z) can not be ignored [24]. Theoretically, NSV terms also show rich perturbative structure like the SV ones. The NSV terms in the inclusive processes have been a topic of interest for more than a decade, see [49,[52][53][54][55][56][57][58][59]. In [50], using physical evolution equations for the CFs of deep inelastic scattering and exploiting renormalisation group invariance, the all order structure of the logarithms were studied. In particular, by suitably modifying the kernels of these equations, certain leading logarithms in NSV terms were correctly reproduced. In [51] (and [20,25]), a remarkable development was made by Moch and Vogt by studying the logarithmic structure for the physical evolution kernels (PEK). They observed that there is an enhancement of single logarithms near threshold up to third order obtained from the fixed order results of DIS, semi-inclusive e + e − annihilation and Drell-Yan production of a pair of leptons in hadron collisions. Based on this finding, they conjectured that it will hold true to all orders in 1 − z, consequently they could predict certain NSV logarithms to all orders in a s like the way the SV terms can be predicted. Very recently, in [87] the resummation of NSV terms has been addressed in the context of semi-inclusive DIS(SIDIS) process. The NSV corrections to various inclusive processes were also studied in a series of papers [49,[52][53][54][55][56][57] and lot of progress have been made leading to a better understanding of the underlying physics. Recently, in [60], we investigated the structure of NSV terms present in the quark anti quark initiated channels in the inclusive production of a pair of leptons in Drell-Yan process and gluon/bottom anti bottom initiated ones for production of Higgs boson exploiting the factorisation properties and renormalisation group invariance along with the certain universal structure of real and virtual contributions obtained through Sudakov K+G equation. We also studied the allorder perturbative structure of the NSV logarithms in the CFs of deep inelastic scattering (DIS) and semi-inclusive e + e − annihilation (SIA) processes in [88]. The formalism was later extended for studying the all-order behaviour of NSV terms in rapidity distributions of Drell-Yan and Higgs production through gluon fusion and bottom quark annihilation in [89]. The goal of [60] was to set up a theoretical formalism to systematically include the NSV contributions in the diagonal channels of inclusive cross sections in general. We used collinear factorisation and studied the structure of various building blocks that constitute the CFs near threshold. The building blocks include the form factors from 2 → 1 processes, the real emission contributions and the Altarelli-Parisi kernels. For the diagonal channels, for SV and NSV terms, the collinear factorisation of bare partonic channels remarkably simplifies. For example, one observes that only diagonal part of collinear singular Altarelli-Parisi kernels containing only diagonal splitting functions contribute. We then use a set of differential equations that govern these building blocks. For instance, the form factors and real emission contributions satisfy K+G equation while AP kernels satisfy AP evolution equations. In addition, each of these terms is renormalisation group invariant. Note that these differential equations are governed by set of universal anomalous dimensions and a few process dependent constants. The solutions to these differential equations can be used to obtain an all order result that encapsulates both SV and NSV terms. The all order result can be displayed in a compact form through an integral representation in z space. The integral representation can be used to resum order one terms in Mellin N space to obtain reliable theoretical predictions at colliders.
For the Higgs boson production in gluon initiated channels, the sum of SV and NSV contributions is denoted by ∆ g , Following [60] in dimensional regularisation (n = 4 + ε),we obtain where Ψ g is finite in the limit ε → 0 and is given by where Z U V,g is the overall renormalisation constant. The constantâ s =ĝ 2 s /16 is the bare strong coupling constant of QCD. The scale µ results from dimensional regularisation. The form factorF g is from virtual contributions to g + g → H while the soft-collinear function Φ g is obtained from real emission subprocesses normalised by square of the form factor. The diagonal AP kernels are denoted by Γ gg which contain only diagonal splitting functions P gg .
The symbol "C" refers to convolution which acts on any exponential of a function f (z) takes the following expansion: We keep only SV distributions namely δ(1 − z), D i (z) and NSV terms ln i (1 − z) with i = 0, 1, · · · resulting from the convolutions. The solution to the K+G equation of the form factor is expressed in terms of cusp (A g ), soft (f g ), collinear (B g ) anomalous dimensions and process dependent constants, while from the AP evolutions equations one finds the solution in terms of cusp (A g ), collinear (B g , C g , D g ) anomalous dimensions. The K+G equation of Φ g contains kernels that contain δ(1 − z), + distributions and ln(1 − z) terms. Due to the explicit z dependence in Φ g , one needs the knowledge of logarithmic structure of CFs, ∆ g , from fixed order perturbative results. Thanks to CFs that are known to third order, we have successfully parametrised Φ g in terms of a set of specific functions of z in dimensional regularisation. The function Φ g is also found to be dependent on the cusp (A g ), soft and the collinear (C g and D g ) anomalous dimensions. Except the overall renormalisation constant the remaining terms are independently IR singular. However, the sum of these contributions is finite and it takes the following integral representation in z space: The z independent coefficient C g 0 is expanded in powers of a s (µ 2 R ) as The AP splitting function P gg is The constants C g and D g are known to three loops in QCD [90,91] ********************************* BORN ************************************* (see [90][91][92][93][94][95][96][97][98][99] for the lower order ones). The cusp, soft and the collinear anomalous dimensions are all expanded in powers of a s (µ 2 F ) as: where X g i to third order are available in [90,91] and are listed in appendix A. The SV function G g SV a s (q 2 (1−z) 2 ), is related to the threshold exponent D g a s (q 2 (1− z) 2 ) via Eq.(46) of [17]. The NSV function ϕ f,g is given by where q 2 = m 2 H and the coefficients ϕ g,1 = 4C A , * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * BORN * * * * * * * * * * * * * * * * * * * * For SU(n c ) gauge group, the color factor C F = (n 2 c − 1)/2n c and C A = n c . Here n f is the number of active flavors.

Resumming next-to SV terms
The resummation of threshold logarithms can be conveniently done in Mellin space. The threshold limit z → 1 in z space translates to N → ∞ in the Mellin space. To include both SV and NSV terms we need to keep ln N as well as O(1/N ) terms in the large N limit. In the large N limit one encounters order one terms from ω = 2β 0 a s (µ 2 R ) ln N at every order in a s in the exponent which spoils the reliability of fixed order truncation of the perturbative series. This can be solved by reorganising the series by using the resummed strong coupling constant. Up to O(1/N ) in N space, the N space resumed will reproduce the fixed order result. The Mellin moment of ∆ g was computed in [60] and is given by We split all the 1/N 0 and 1/N terms in Ψ g N as In other words, we split Ψ g N in such a way that Ψ g sv,N contains ln j N, j = 0, 1, · · · terms while Ψ g nsv,N contains terms of the form (1/N ) ln j N, j = 0, 1, · · · . It is well known that Ψ g sv,N takes the following form: The exponents g g 0,i and g g i are listed in the Appendices B.1 and B.3. . The g g i coefficients are universal whereas the constant g g 0 contains only N independent pieces resulting from the Mellin moment of Ψ g D (µ 2 F , z). The constant g g 0 is ambiguous and is defined by the condition that Ψ g sv,N − ln(g g 0 ) = 0 when N = 1. This ambiguity can be exploited to define various resummation schemes. The N -independent coefficients C g 0 and g g 0 are related to the coefficientsg g 0 given in the paper [100,101] using the following relation, which is expanded in terms of a s (µ 2 R ) as, whereḡ g i (ω) and h g ik (ω) are presented in the Appendices B.4 and B.5, respectively. We find that in each of the exponents g g i (ω), g g i (ω) and h g ik (ω) from SV as well as NSV, we are resumming in the Mellin space "order one" ω terms to all orders in perturbation theory. In the large N expansion, one finds that the ln N terms are always associated with Euler-Mascheroni constant, γ E . As these γ E terms can spoil the reliability of perturbative expansion, they are often combined with ln N . This results in a change in the definition of order one term ω as This alters the numerical predictions considerably, see [102]. We call the former one as standard N-exponentiation and the latter by N -exponentiation. In addition to these γ E terms, one encounters N independent pieces that are mostly proportional to irrational terms namely Euler zeta terms ζ i . Unlike γ E , there exists no systematic way to sum them up. In other words, there exists ambiguity in whether one should keep them in the exponent or not. While there are many ways we can treat these terms, we adopt two distinct ways namely soft exponentiation and all exponentiation. In soft exponentiation, the Mellin moment of entire soft-collinear function [17,103] which contains all the + distributions as well as the δ(1 − z) terms is exponentiated. In the latter scheme , we exponentiate the form factor along with the softcollinear function in the Mellin space. In [44,45], this approach was used to study the inclusive Higgs boson production in gluon fusion. They found that the numerical predictions are less sensitive to the choice of factorisation and renormalisation sales compared to the standard threshold approach. This approach is theoretically supported by the fact that the Sudakov K+G differential equation [17,[103][104][105][106][107] that the form factors satisfy give solutions of the form of exponential of N independent terms.
In [61,101], we have presented various schemes that can deal with N independent terms for the numerical study of Drell-Yan process. We briefly describe them below for completeness: • Standard N exponentiation: In this scheme, we exponentiate only the large-N pieces that contribute to the CFs while keeping all the N -independent pieces outside the exponential. This will only alter the Ψ g sv,N whereas the Ψ g nsv,N will remain same as it contains only N -dependent terms. Therefore we write: where NSV part is defined in (2.37). And the SV part is given as, which can be obtained from the Mellin moment of Ψ g D and keeping only those terms that vanish when N = 1.
• Standard N exponentiation: Here, we compute Mellin moment in the large N = N exp(γ E ) limit that means the entire γ E dependent terms are exponentiated through N . In this scheme, the CFs ∆ g,N can be expressed as where g g 0 results from the N -independent part of the Mellin transformation of the distributions and the γ E 's are all absorbed into the N -dependent functions. Here the N -dependent part of SV, i.e., G g sv,N , is obtained by setting Ψ g sv,N − ln g g 0 = 0 when N = 1. Hence, G g sv,N and Ψ g nsv,N are given by Numerically, this can make difference at every logarithmic accuracy. In [102], it was found that the N -exponentiation shows a faster convergence compared to the N -exponentiation for the charged and neutral DIS processes. This is the standard approach adopted all along this paper and we compare other prescriptions with the N -exponentiation scheme in later sections.
• Soft exponentiation: Here, we exponentiate complete finite part of soft-collinear function Φ g to obtain Here the exponents g g,Soft i contain the complete finite part coming from the softcollinear function. The remaining N -independent terms coming from finite part of form factor and AP kernels contribute tog g,Soft 0 , whose expansion in powers of a s is given as:g The N -independent constantsg g,Soft

0i
(m 2 H ) and the resummed exponents g g,Soft i (ω) are listed in Appendix D.
• All exponentiation: This scheme arises in light of (2.12) which is a consequence of the differential equations satisfied by each of the building blocks leading to an all-order exponential structure for ∆ g,N . Here we study the numerical impact of the entire result taken to the exponent as given below.
The exponent in (2.49) contains both N dependent and independent terms which are listed in Appendix E. In [44,45] this scheme was used to study the inclusive cross section for the production of Higgs boson in gluon fusion at the LHC. For similar studies on DY and DIS processes in MS schemes, see [108].
In [61,101], we had studied how various schemes discussed so far can affect the predictions of invariant mass distribution of lepton pairs in the DY process. In the present work, we perform a similar analysis with the resummed NSV exponent for the Higgs boson production and report their numerical impact. In order to distinguish between SV and SV+NSV resummed results, we denote the NSV included resummed results by N n LL. We use fixed order results up to NNLO level and the resummed results up to NNLL accuracy. The resummed results are matched to the fixed order result in order to avoid any double counting of threshold logarithms. The resummed result at a given accuracy, say N n LL, is obtained by taking the difference between the resummed result and the same truncated up to order a n s . Hence, it contains contributions from the SV and NSV terms to all orders in perturbation theory starting from a n+1 s on wards: (2.51) We use PDFs known to NNLO accuracy and their Mellin space counterparts (f i,N ) can be obtained using QCD-PEGASUS [109]. However, we follow [14,43] to directly deal with PDFs in the z space. In the above eq. (2.51) the second term represents the resummed result truncated to N n LO order. In Table 1, we list the resummation coefficients which are required to obtain the resummed ∆ g,N at a given logarithmic accuracy. In the next section, we present the numerical results for the production Higgs boson in gluon fusion at the LHC.
Logarithmic Accuracy Resummed Exponents Table 1: The set of resummation coefficients g g 0,i , g g i (ω), g g i (ω), h g i (ω) which is required to compute ∆ g,N at a given logarithmic accuracy.

Numerical Analysis
In this section, we discuss the numerical impact of resummed soft virtual plus next-to-soft virtual (SV+NSV) corrections for the Higgs boson production in the gluon fusion channel at the LHC. The leading order process is initiated via gluon-gluon fusion through topquark loop, which is integrated out in the large m t -limit. Throughout the computation we use, the top quark mass m t = 173.2 GeV, the Higgs mass m H = 125 GeV and the Fermi constant G F = 4541.63 pb. The fixed order numerical predictions up to NNLO have been obtained using an in-house FORTRAN code and we have validated our results with the publicly available fixed order code ggHiggs [45]. We have performed the inverse Mellin transformation of the resummed N -space result using an in-house FORTRAN code. Minimal prescription [110] has been used to deal with the Landau pole in the Mellin inversion routines. We use the MMHT2014(68cl) [111] parton distribution set and LHAPDF [112] interface to provide the strong coupling constant a s with n f = 5 active massless quark flavours throughout. The detailed analysis is done for 13 TeV LHC, however in the later sections we have extended our analysis to other collider energies as well.
We begin by investigating the importance of NSV logarithms in comparison to the SV distributions, present in the fixed order predictions for the gluon induced Higgs boson productions. Hence, we restrict ourselves to the partonic coefficient function ∆ g given in Eq. 2.12 which consists of SV and NSV terms. The SV part consists of δ(1 − z) and D k (z) while the NSV part consists of collinear logarithms ln k (1 − z). We compute the percentage contribution of the NSV logarithms to the hadronic cross section to analyse their relevance. At NLO, we find the SV comprises 73.16% of the Born contribution while the NSV contribute to only 45.81%. However, from NNLO onwards the contributions coming from NSV logarithms become more significant than that of SV. For instance, at NNLO the NSV logarithms give 58.91%, whereas SV 15.81% of the Born cross-section. Similar trend follows at N 3 LO level where NSV contribution is 25.83% and SV has a negative contribution of −2.29% of the total Born cross section. In Table 2, we provide the individual contributions coming from different powers of SV distributions and NSV logarithms at NNLO accuracy. The corresponding percentages for SV and NSV contributions at N 3 LO accuracy are given in [48]. From this analysis, one can conclude that with increasing order of perturbative corrections, NSV logarithms become phenomenologically more relevant than the SV distributions.
Moreover, in case of Higgs production through gluon fusion, the diagonal gg-channel constitutes the dominant contribution to the hadronic cross section. For example, the contribution from the gg channel is 48.35% of the full NLO cross section whereas the correction from the qg-channel is only −0.79%. Similarly the a 2 s correction from the gg channel is 19.71% and from the qg is around 2% of the NNLO cross section. The other sub-dominant channels are more suppressed in comparison to the qg-channel. Since there is a significant increase in the percentage contribution of NSV logarithms over SV from NNLO accuracy on wards from the diagonal channel, which is also the main contributor to the hadronic cross section, it would be interesting to see how phenomenologically important the resummed NSV logarithms are, in addition to the well established SV resummation.
In the context of threshold resummation in [45], it was found that the resummation of the SV distributions to N 3 LL accuracy change the NNLO results by 20% at µ R = µ F = m H and 5.5% at µ R = µ F = m H /2. Compared to fixed order predictions, the SV resummed results were found to depend moderately on the renormalization scale while there was increase in the dependence on the factorization scale. Recently in [61], the SV resummation was extended to include the resummed effects from the NSV logarithms for the lepton pair production in Drell-Yan process. It was found that owing to the large coefficients of the NSV logarithms, the resummation of such logarithms lead to a considerable amount of enhancement of the fixed order cross section. The renormalization scale dependency were also seen to reduce compared to the fixed order results. However, the uncertainty from factorization scale increases as we did not include the resummed contributions of NSV terms, namely the collinear logarithms from the off-diagonal channels.
In summary, we find that in the dominant gg channel, as the order of perturbation increases, the numerical coefficients of collinear logarithms increase leading to contributions that are larger than those from SV distributions. In addition, there exists theoretical results for the resummation of NSV terms to all orders up to NNLL accuracy [60] and hence using the latter, it is desirable to study their numerical importance. In this context, we ask the following questions that can shed light on the importance of NSV terms and also unravel the role of beyond NSV terms in the context of Higgs boson production in gluon fusion at the LHC.
• How much is the enhancement of the resummed cross section in comparison to the fixed order ones, owing to the large coefficients of the NSV logarithms ?
• Since gg-channel is the dominant contributor to the hadronic cross section, what is the impact of resummation of NSV terms from the collinear logarithms on the unphysical scales, namely µ F and µ F ?
• How different is the impact of SV+NSV resummation from the well established threshold resummation ?
We address the above questions in the following sections.

Fixed-order vs Resummed results
In this section, we study the numerical relevance of resuming the NSV logarithms along with the SV distributions at LL, NLL and NNLL accuracy matched to the corresponding fixed order results using Eq. 2.51. We investigate the dependence of SV + NSV resummed cross section on the choices of renormalisation µ R and factorisation µ F scales to study their perturbative uncertainties. We start by plotting the inclusive cross section as a function of both µ R and µ F simultaneously in a 3D graph for centre-of-mass energy 13 TeV. at LO + LL, NLO + NLL and NNLO + NNLL accuracy at (a), (b) and (c) subpanels respectively. The scales µ R and µ F are varied in the range 0.
The fixed order predictions contain contribution from all the channels while the resummed result consists of logarithms coming only from the diagonal gg channel. There is a comprehensible reduction in the scale dependency of the cross section when we go from LO to NNLO. Thus as expected, the scale dependence is reduced when the higher-order corrections are included in fixed order result. Whereas, from the resummed predictions on the right panel of Fig.1, we find that the scale dependency increases in comparison to the corresponding fixed order predictions. However, if we consider the region in the vicinity of {µ R , µ F } = { m H 2 , m H 2 } and vary these scales in such a way that the ratio µ R µ F is not larger than 2 or smaller than 1/2, then we find that the scale dependency of the resummed result becomes comparable to the corresponding fixed order result at NNLO accuracy. The overall scale uncertainty using the canonical 7-point variation method, ranges between (+8.92%, −10.12%) for NNLO whereas it is (+11.90%, −8.32%) for NNLO + NNLL which validates our former statement. In addition to this, there is a significant enhancement of 4.76% to the NLO cross section due to the addition of resummed result at NLL accuracy. Similarly, the inclusion of the SV+NSV resummed result decreases the NNLO cross section by 3 with an overall scale uncertainty comparable to the fixed order ones. There is also a better perturbative convergence of the resummed result at this particular choice of central scale in comparison to the corresponding fixed order results. Henceforth, in rest of the paper, because of the above arguments, we keep the central scales at We would also like to mention that the resummation scheme chosen for the rest of the paper is N exponentiation. The reasons for this particular choice of scheme will be discussed in the last section of this paper.
We now proceed to study the quantitative impact of the SV + NSV resummed corrections against the fixed order results using the K-factor, defined by  Few comments are in order from the above table. We notice that the percentage enhancement in the cross section at NNLO + NNLL is 9.59% over NLO + NLL whereas it is 18.55% when we go from NLO to NNLO accuracy. Therefore, inclusion of resummed SV+NSV logarithms improves the reliability of perturbative predictions. Also, the K-factor values are closer for NNLO and NNLO + NNLL as compared to NLO and NLO + NLL. This suggests that the resummed contributions to the fixed order cross section decrease as we go to higher orders in perturbation theory.
From the above discussion on the K-factors, we have seen that the SV + NSV resummed contribution enhances the cross section as compared to the corresponding fixed order predictions till NLO accuracy. The 3D graph in Fig.1 shows the overall increase in the µ R and µ F scale variations when we add the resumed predictions till NLL accuracy. However, at NNLO + NNLL accuracy, the cross section decreases and the overall scale uncertainty gets closer to the corresponding fixed order prediction in the vicinity of {µ R , µ F } = {m H /2, m H /2}. This suggests that the truncation at NNLO improved by resummed contributions at NNLL accuracy provides a reliable perturbative predictions. In the following, we study the impact of µ R and µ F scales separately by keeping one of them fixed. We start by plotting the cross sections at various orders of fixed order as well as resumed result as a function of µ R and µ F in the range (0.2, 3)m H in Fig.2. The quantitative values of uncertainty is calculated in the range (1/4, 1)m H such that the ratio of µ R µ F is not greater than 2 or smaller than 1/2. In the first panel, the cross section is varied w.r.t µ R , keeping µ F held fixed at m H 2 . The uncertainty ranges between (+17.21%, −14.55%) for NLO + NLL whereas it is between (+23.21%, −17.18%) for NLO. Similarly, it ranges between (+6.54%, −8.32%) for NNLO + NNLL whereas for NNLO, it is between (+8.85%, −10.12%). Thus, we observe that µ R uncertainty decreases when the SV + NSV resummed contributions are added to the fixed order predictions. Also, it decreases significantly when we go from NLO + NLL level to NNLO + NNLL level. The maximum increments and minimum decrements in the cross section from the value at the central scale µ R = µ F = m H 2 in the range µ R = {1/4, 1}m H are given in Table 4 for both fixed order as well as resummed results. From Fig.1, we observed that the scale uncertainty arising by varying both the scales simultaneously was larger for the resummed cross section compared to the corresponding fixed order results. However, here we notice that µ R dependency is less for resummed cross sections indicating the role of µ F in the scale uncertainty of Fig.1  The second panel of Fig.2 shows the dependence of cross section on factorisation scale µ F keeping µ R held fixed at m H /2. We notice that the fixed order predictions till NNLO order have negligible dependence on µ F . This is also reflected from the maximum increments and minimum decrements from the central value (for µ R = µ F = m H /2) for fixed order predictions given in Table 5. However, in the case of resummed cross section, Table 5 shows significant increments and decrements around the central value. For instance, for the LO + LL, the uncertainty lies between (+4.51%, −6.92%) while for NLO + NLL it changes to (+32.17%, −18.62%). At NNLO + NNLL accuracy, the µ F uncertainty comes down to (+11.90%, −7.56%). We have seen previously in Fig.1 that the overall scale uncertainty due to renormalisation and factorisation scales was maximum for NLO + NLL case. This suggests that the scale uncertainties for the resummed result is driven by µ F variation. But we also know that µ F variation entangles contribution arising from different channels. Hence, the naive expectation is that the resummed predictions for the off-diagonal collinear logarithms may be accounted as the compensating factor for the µ F dependency coming from diagonal channels. However, this can not account for µ F uncertainty in the case of Higgs production through gluon fusion because the off-diagonal channels have minuscule contributions at each order in perturbation theory. We will interpret this in the next section by studying the behaviour of SV resummed results which contain only distributions and delta functions.  In the third panel of Fig.2, we vary µ R and µ F simultaneously by setting µ R = µ F = µ. This plot basically is the diagonal line from the axis in the 3D graph of Fig.1 where µ R = µ F . We observe that for fixed order predictions, the overall uncertainty is dictated by the behavior of µ R uncertainty. On the other hand, in the case of resummed predictions, it is the µ F dependency which plays the major role. In Table 6 given below, we have presented the values of cross section for fixed order as well as resummed predictions for central scale value µ R = µ F = m H /2.
The maximum increment and decrements around central scale value are calculated by setting µ R = µ F = µ and varying µ.
In summary, we studied the behaviour of resummed SV + NSV logarithms w.r.t µ R and µ F scales and compared it against fixed order results. We found that in comparison to the fixed order the µ R uncertainty decreases for the resummed results. But for the µ F variation we found that the uncertainty is more at NLO + NLL than at NNLO + NNLL with larger factorization scale dependence compared to their fixed order counterparts. In addition, we compared the K-factors of resummed results against the fixed order predictions to understand the perturbative convergence. In the next section we focus on the role of resummed NSV terms against SV resummed contribution.

SV resummation vs SV + NSV resummation
In the previous section, we compared the SV + NSV resummed results with the fixed order predictions. The K-factors in Table 3 show that there is a significant enhancement till NLO+NLL order. We find that the contribution from resummed SV+NSV terms decreases the fixed order NNLO predictions when matched at NNLO + NNLL level leading to better perturbative convergence. We also observed that µ R scale uncertainty gets improved at each order upon including the resummed corrections. On the other hand, uncertainty from µ F scale increases significantly by the inclusion of resummed cross section especially at NLO + NLL accuracy. Now, in this section we will perform a detailed analysis on the inclusion of resummed NSV logarithms by comparing it with the SV resummed results. We will also try to reason out the behaviour of SV + NSV resummed results w.r.t µ F and µ R uncertainty variations. NLO We begin with the K-factors for SV resummed results as well as for SV + NSV resumed results at NLO and NNLO accuracy given in Table 7. We find that the inclusion of resummed NSV logarithms enhances the cross section by 8.01% when we go from NLL to NLL accuracy. Whereas, at NNLO + NNLL the cross section decreases slightly by 0.42%. This shows the reduction of resummed corrections while taking into account the higher logarithmic accuracy. We also observe that the K-factors are closer for NLO + NLL and NNLO + NNLL as compared to NLO + NLL and NNLO + NNLL. This suggests better perturbative convergence when NSV logarithms are taken into consideration.
We now move on to study the scale uncertainties arising from the NSV logarithms. For this, we do a comparative study of µ R and µ F variations separately between fixed order, resummed SV and resummed SV + NSV results. In the first panel of Fig. 3, the cross section is plotted as a function of µ R keeping µ F = m H /2. The µ R uncertainty for NLO lies in the range (+23.21%, −17.18%), for NLO + NLL between (+18.57%, −14.99%) and for NLO + NLL, it lies between (+17.21%, −14.55%). Similarly, at second order, the uncertainty lies between (+8.85%, −10.12%) for NNLO, (+9.58%, −10.02%) for NNLO + NNLL and for NNLO + NNLL, it ranges between (+6.54%, −8.32%). We notice that at NLO accuracy, inclusion of resummed SV as well as resummed SV + NSV contributions improve the renormalisation scale uncertainty. Thus as expected, the inclusion of higher order logarithmic corrections within a particular channel reduces the µ R uncertainty. We also observe that the uncertainties at NLO + NLL and NLO + NLL are comparable. For instance, the fixed order result at NLO consists of 73.16% of SV contribution and 45.81% of NSV contribution. Thus, NSV being the sub-dominant contributor does not allow much improvement in the uncertainty when resummed NSV logarithms are added to the resummed SV distributions. Now, at NNLO accuracy, the inclusion of resummed SV contribution does not improve the µ R scale dependency significantly over the fixed order result. However, at the NNLO + NNLL, there is a comprehensible reduction in µ R scale variation. We find that the SV distributions at two-loop contribute to 15.81% of the Born cross section, while the NSV logarithms at the same order contribute to overall 58.91% to 0.5  the fixed order result at NNLO. This suggests a significant improvement in the µ R scale variation while including the resummation effects of NSV logarithms. From the above analysis we find that the resummed result leads to a reduction of the renormalization scale uncertainty. Let us now take a deeper look into this reduction. In general any resummed result carries informations of the all-order corrections arising from the quantity we are resumming. Further, owing to the "inexact" Mellin inversion, it also carries certain spurious terms which are beyond the precision of the resummed quantity. For instance, the resummed NSV logarithms at NNLO + NNLL accuracy incorporates to the fixed order results the all-order corrections arising from the summation of the nextto-next-to leading towers of NSV logarithms and certain spurious terms which are beyond NSV. Now looking at the reduction of the scale uncertainty from NNLO to NNLO + NNLL to NNLO + NNLL, it is clear that the all-order corrections present in the resummed results leads to the improvement in renormlization scale uncertainty. Now this reduction is more at NNLO + NNLL than at NLO + NLL owing to the significant contribution of the NSV logarithms at NNLO. In Table 8 Having established the reasons for the behaviour of resummed NSV logarithms w.r.t µ R variation, we proceed to study the scale uncertainty w.r.t the factorisation scale µ F . The second panel in Fig.3, shows the variation in cross section as a function of µ F keeping µ R = m H /2. Let us first consider the plot in Fig.4 where we depicted the fixed order result truncated to SV+NSV accuracy against µ F variation. We find that there is a significant increase in scale dependency when we go from NLO to NNLO accuracy. Recall that the µ F dependence in the fixed order predictions given in Fig. 3 is extremely mild. This suggests that the large µ F dependence that we see in in Fig.4 is expected to get compensation from beyond NSV terms in the threshold expansion. This is because the fixed order analysis confirms the significant contributions arising from beyond SV terms with increase in the order of perturbative expansion [48].
Our next observation is on the behaviour of resummed SV results against µ F variation. We find that the uncertainties at NLO + NLL and NNLO + NNLL lie in the range (+21.68%, −14.84%) and (+5.90%, −6.18%) respectively. Interestingly, we observe that the resummed SV cross section at NLO + NLL is more sensitive to µ F compared to fixed order NLO results. It is mainly due to the spurious terms that contributes beyond SV, arising from the "inexact" Mellin inversion of the N -space resummed result. Unlike NLO + NLL, at NNLO + NNLL, the µ F dependence drops down significantly but is still more than the fixed order counterpart. Hence the resummation of next-to-next-to leading SV distributions, compensates the spurious terms of NLO + NLL with higher order logarithmic corrections. But again the reduction in the uncertainty will not be as much as compared to the fixed order due to the presence of the residual spurious terms arising from the Mellin inversion at NNLO + NNLL accuracy. Now let us try to understand the behaviour of NSV logarithms based on our findings in fixed order and resummed SV predictions. Recall that the resummed SV results suggest that spurious beyond SV terms spoil the µ F scale dependence. Also, the behaviour of fixed order results w.r.t µ F and the Fig.4 suggest that µ F dependence of both NSV and beyond NSV terms increase with order of perturbative expansion and hence any truncated result will have more uncertainties than the fixed order corrections. This explains the µ F dependence of SV + NSV resummed prediction given in the second panel of Fig.3.
In addition, we also observe in Fig.3 that the µ F dependency goes down from (+32.17%, -18.61%) at NLO + NLL to (+11.89%, -7.56%) at NNLO + NNLL. If we compare the uncertainty of the SV resummed result, which was in the range of (+21.68%, −14.84%) at NLO + NLL we can conclude that the bulk of the uncertainty in SV+NSV resummation at this particular accuracy arises from the SV resummation. Now at the next logarithmic accuracy, the uncertainty coming from the SV resummation, which is (+5.90%, −6.18%) at NNLO + NNLL, is comparable to the uncertainties arising from the resummation of the NSV logarithms. It has be noted that resummation of NSV logarithms was never supposed to compensate the uncertainties arising from the resummation of SV distributions. But owing to the "inexact" Mellin Inversion, the spurious terms which are developed in the SV resummation at the NSV accuracy gets corrected through NSV resummation while beyond NSV spurious terms from both the distributions and logarithms will remain uncompensated and aid to the µ F uncertainty. Also similar to the SV resummation, higher order logarithmic corrections from the NSV terms will compensate for the lower spurious ones thereby reducing the uncertainties even more. In Table 9, we have provided the fixed order as well as resummed cross section values at the central scale for various perturbative orders. We have also given the maximum increment and decrements at each order around the central value by varying µ F in the range { 1 4 , 1}m H keeping µ R = m H /2 fixed. These are used to calculate the corresponding percentage uncertainties given above.
In summary, we find that the resummation of NSV terms reduces the uncertainty resulting from µ R when the logarithmic accuracy is improved from next-to-leading to nextto-next-to-leading order. As we know, any change resulting from the µ R variation at a given perturbative order will be compensated by the remaining higher order terms. Hence the resummation naturally reduces the dependence on this unphysical scale, as we sum up certain SV and NSV logarithms to all orders, which include higher order µ R dependent  terms as well. However, this is not case for the µ F sensitivity. Unlike µ R dependent terms, the µ F dependent terms do not cancel between different perturbative orders. On the contrary, we find that at each perturbative order SV, NSV and beyond NSV terms have different µ F dependence. Hence truncation of higher order threshold terms, or equivalently approximation based on threshold expansion can lead to residual µ F scale dependency at a particular order. This gets amplified further due to resummation where only certain SV and NSV logarithms are resummed to all orders due to approximation in the threshold expansion. This also hints towards the importance of beyond NSV terms for better accuracy of predictions at each level. In the next section, we study the importance of SV + NSV resummed contributions to inclusive cross section at various collider energies.

Numerical results for different collider energies
From Fig.1, we deduced that the appropriate central scale for the analysis would be µ R = µ F = m H 2 . In literature, one finds the choice µ R = µ F = m H as well to be the central scale. We will do a short analysis in this section to justify our choice of central scale. For this, we will vary the cross section as a function of hadronic centre-of-mass energy. It will also help us to understand the behaviour of SV + NSV resummed cross section at different collider energies.
In Fig.5, we plot the SV + NSV resummed cross section against the hadronic centre-ofmass energy from 7 TeV to 14 TeV. The bands in the plot correspond to the scale variations obtained by using the canonical 7-point variation where {µ R , µ F } is varied in the range { mx 2 , 2m x }, keeping the ratio µ R µ F not larger than 2 and smaller than 1/2. For the left panel, m x = m H and for the right panel, m x = m H /2. At 13 TeV LHC, the enhancement in the cross section from the central scale µ R = µ F = m H values while going from LO + LL to NLO + NLL is 37.99% and it is 5.60% from NLO + NLL to NNLO + NNLL. The large increment in the cross section for NLO+NLL over LO+LL hints towards inclusion of more and more higher order terms and therefore results in loss of credibility of the perturbative prediction. If we look at these enhancements for central scale µ R = µ F = m H /2 values, it is −4.04% while going from LO + LL to NLO + NLL and 9.59% for NLO + NLL to NNLO + NNLL. Also, we observe that for the plot corresponding to the central scale µ R = µ F = m H /2, the NLO + NLL scale uncertainty band is included within the LO + LL band and the NNLO+NNLL scale uncertainty band is included within that of NLO+NLL. However, this is not the case with the plot corresponding to central scale µ R = µ F = m H . Thus, the above two observations regarding the enhancement and the uncertainty band plot indicate that the perturbative expansion of the hadronic cross section is more convergent and therefore reliable for the right panel.
From the above plots, we also observe that the uncertainty due to renormalisation and factorisation scales increases with the energy of the collider for both the choices of central scales. The reason for this large uncertainty at high E CM could be the lack of knowledge of the PDF sets at these energies.
Let us now look at this theoretical uncertainty related to the choice of central scale for SV + NSV resummed cross section by comparing it with the fixed order and SV resummed results. Tables 10 and 11 present the cross sections and related scale uncertainty for fixed order, SV and SV + NSV resummed results at NNLO accuracy for different collider energies for the central scales µ R = µ F = m H and µ R = µ F = m H /2 respectively. The scale uncertainty in the tables has been calculated by varying the renormalisation and factorisation scales and using the canonical 7-point variation approach.
We start by calculating the quantitative effect of including the resummed SV and SV + NSV logarithms to the fixed order prediction at the above mentioned central scales. We find that at 13 TeV LHC, the inclusion of resummed SV + NSV result increases the NNLO prediction by 14 Table 11: Fixed-order+ Resummed gluon fusion cross section at different centre of mass energies at the LHC. The scale uncertainty has been estimated using 7-point scale variation around the central scale µ R = µ F = m H /2.
Next, we observe that the scale dependence due to µ R and µ F is less in the case of fixed order as well as SV and SV + NSV resumed results if the central scale is µ R = µ F = m H /2. Though the uncertainties don't differ comprehensibly for the fixed order and SV resumed results but the difference is significant in the case of SV + NSV resumed result between the two different choices of central scales. At 13 TeV LHC, for the central scale µ R = µ F = m H , the uncertainty for NNLO + NNLL accuracy lies between +19.65%, −15.10% whereas this goes down significantly to +11.90%, −8.32% for the central scale µ R = µ F = m H /2.
From the tables(10,11) given above, we can also see that the relations between fixedorder, SV and SV + NSV resummed results and their uncertainty bands are similar for all the collider energies considered in the table. For better understanding, in Fig.6, we present the plots for fixed-order, SV and SV + NSV resumed cross sections at NNLO accuracy, at various collider energies for both the central scale choices. Having established the reasons for the choice of central scale µ R = µ F = m H /2, we now proceed to the last part of our numerical analysis where we will see the effect of SV + NSV resummation in different schemes namely, N exponentiation, N exponentiation, all exponentiation and soft exponentiation.

SV + NSV resummation in different schemes
In this section, we will explore why we chose N exponentiation to do the numerical analysis in the previous sections. We will also do the comparative study between N exponentiation, N exponentiation, all exponentiation and soft exponentiation to see how the enhancements in the cross section at different orders and scale uncertainty gets modified due to SV + NSV resummation across these schemes.
We start by studying the numerical impact of N and N exponentiation. As we know, for the case of N exponentiation, we resum γ E terms along with ln N whereas for N exponentiation we only resum the ln N terms. We have already done detailed analysis of SV + NSV resummation effects in N exponentiation in the previous sections. Let us now see how the behaviour changes by the inclusion of resummed SV + NSV results in the N exponentiation.  For this, we first compare the K-factor values for these two schemes at various orders given in Table 12. We observe that for N exponentiation, there is an increment of 46.09% for NLO + NLL cross section over LO + LL and 10.51% while going from NLO + NLL accuracy to NNLO + NNLL. On the other hand, for N exponentiation as we have seen earlier, there is a reduction of 4.04% in the cross section at NLO+NLL compared to LO+LL and an enhancement of 9.59% while going from NLO + NLL accuracy to NNLO + NNLL. These percentages show that the perturbative convergence for SV + NSV resummed result in N exponentiation is better than N exponentiation.
Next, we look at the uncertainty due to renormalisation(µ R ) and factorisation(µ F ) scales. From Table 13, we find that the uncertainty calculated using 7-point canonical approach by varying µ R and µ F is less for N exponentiation compared to N exponentiation till NLO + NLL accuracy. However, at NNLO + NNLL accuracy, the uncertainty of SV + NSV resummed result for N exponentiation(lies between {+9.19%, −8.68%}) becomes comparable to N exponentiation(lies between {+11.90%, −8.32%}). Hence, the good perturbative convergence is achieved by the resummed result in N exponentiation which improves the reliability of the perturbative predictions and it becomes the deciding factor in choosing N exponentiation over N exponentiation for numerical analysis in the previous sections.  Table 13: A comparison of SV + NSV resummed cross sections (in pb) between different resummation schemes at central scale µ R = µ F = m H /2 for 13 TeV LHC.
In the end we discuss on numerical impact of SV + NSV resummed effects by giving a short analysis on the other two schemes namely, All exponentiation and Soft exponentiation.
In all exponentiation scheme, we exponentiate the Mellin moment of complete soft-collinear function along with the form factor contributions whereas in soft exponentiation scheme, we only exponentiate the soft-collinear part. We find that the SV+NSV resummed result in soft exponentiation is closest to our best prediction given in N exponentiation. From Table  13, we see that the maximum increments and decrements around the central scale value is comparable for soft exponentiation and N exponentiation schemes till NNLO + NNLL accuracy. Also, the enhancement while going from NLO + NLL to NNLO + NNLL accuracy is 3.64% in soft exponentiation scheme which is close to 9.59% increase in the case of N exponentiation scheme as compared to other schemes. We observe a striking feature in the case of all exponentiation scheme. There is an increase in cross section at NNLO + NNLL compared to NLO + NLL accuracy in all the schemes except all exponentiation. In this case, the cross section decreases by 11.28% at NNLO + NNLL accuracy. Also, the scale uncertainties coming from µ R and µ F variations are found to be least in all exponentiation scheme at NLO + NLL and NNLO + NNLL level with latter lying between (+11.40%, −4.26%).

Conclusions
Extensive study on the Higgs bosons produced at the LHC has been going on both in experimental as well as in the theoretical fronts. They provide not only stringent tests of the SM but also constraints on models beyond SM. Variety of observables involving Higgs bosons measured with unprecedented accuracy have been compared against state-of-the art predictions. One such observable is the inclusive cross section of single Higgs boson production and the predictions are known in perturbative QCD to N 3 LO accuracy in fixed order and to N 3 LL in resummed framework. The electroweak corrections at NNLO further improves the predictions. These estimates computed using effective theory approach are mildly improved when exact mass dependence is included.
Further, partial results beyond third order for SV part of the cross section are known and they are found to improve the predictions. In addition, there are remarkable developments to include terms from next to SV contributions beyond third order and also to systematically sum them up to all orders, through resummation framework. Recent numerical studies with resummed NSV contributions at leading logarithmic approximation strongly justify their phenomenological importance. The present paper incorporates NSV terms beyond leading logarithmic approximation, in particular for the dominant diagonal channel beyond NLO level and reports their numerical impact at the 13 TeV LHC. Note that the theoretical set up for the resummation of NSV terms from the non-diagonal channels is still a open problem beyond LL approximation. We expect that their numerical impact will be subdominant compared to the diagonal channel. Fixed order predictions already hints the importance of not only NSV terms but also beyond NSV terms and our study confirms this in the resummed framework. Our numerical predictions on the K factors taking into account the resummed contributions from NSV terms beyond leading logarithmic approximation imply better perturbative convergence compared to fixed order or SV resummed contributions. For instance, the cross-section at NNLO + NNLL accuracy reduces by 3.15% as compared to the NNLO result for the central scale µ R = µ F = m H /2 at 13 TeV LHC. In addition, the dependence on the renormalisation scale gets reduced upon the inclusion of NSV resummed making the prediction more reliable. However, we find that the sensitivity to factorisation scale increases in the presence of resummed NSV terms implying the importance of beyond NSV terms within the resummed framework. It is well known that the SV resummed predictions depend on how we treat N independent terms leading to different schemes. We have compared numerical predictions from these schemes in the presence of NSV terms and found that N exponentiation scheme shows better perturbative convergence and less sensitive to unphysical scales. In summary, we have presented a detailed study on the numerical importance of resummed NSV terms for the production of Higgs bosons at the LHC. We find that the inclusion of resummed NSV terms improves perturbative convergence and reduces the uncertainty from the choice of renormalisation scale.

Acknowledgements
We thank Claude Duhr, Moch and Bonvini for useful discussions throughout this project. We also thank Claude Duhr and Bernhard Mistlberger for providing third order results for the inclusive reactions. We would like to thank L. Magnea and E. Laenen for their encouragement to work on this area. In addition we would also like to thank the computer administrative unit of IMSc for their help and support.

A Anomalous dimensions
Here we present all the anomalous dimensions used in performing the resummation.

Cusp anomalous dimensions A g i
In the following we list the cusp anomalous dimensions A g i till four-loop level: where n f is the number of active quark flavours in the theory. The quadratic Casimirs C F and C A are given by n 2 c −1 2 nc and n c respectively. The quartic Casimirs are given by The collinear anomalous dimensions B g i are given till three-loop as, Soft anomalous dimensions f g i The soft anomalous dimensions f g i till three-loop are given as The NSV anomalous dimensions C g i and D g i till three-loop are given as (A.14)

SV Threshold Exponents
The function G g SV a s (q 2 (1 − z) 2 ), given in 2.16 is related to the threshold exponent D g a s (q 2 (1 − z) 2 ), via Eq.(46) of [17] where the universal D g i coefficients, till three-loop, are given as, The N -independent coefficients g g 0,i in Eq. 2.34, till three-loop, are given by Here, γ E is the Euler-Mascheroni constant and ζ i 's are the Riemann zeta functions. In the aforementioned equations, L qr = ln q 2 µ 2 R and L f r = ln The N -independent coefficientsg g 0,i in Eq. 2.36, till three-loop, are given bỹ γ E L f r ζ 3 + 15008 9 γ E L f r ζ 2 The resummation exponents g g i in Eq. 2.33, till three-loop, are given by Here, L ω = ln(1 − ω) and ω = 2β 0 a s (µ 2 R ) ln N . And β i 's are the QCD β functions which are given by B.4 The NSV resummed exponents g g i The resummation exponents g q i in 2.37, till three-loop, are given by

C Resummation Coefficients for the N exponentiation
For the case of N exponentiation, all the N -dependent resummed exponents namely g g i (ω), g g i (ω), and h g i (ω) given in Eqs 2.43 and 2.45, respectively can be obtained from the corresponding exponents in standard N -approach through setting all the γ E terms to zero and replacing all the ln N terms by ln N as well as all the O(1) ω terms by ω as mentioned in Sec. 2.3. The N -independent constantsḡ g 0 given in Eq. (2.42) can be obtained from their counterparts in standard N -approach by simply putting the γ E terms equal to zero.

D Resummation Coefficients for the Soft exponentiation
For the case of soft exponentiation, all the terms coming from the soft-collinear function Φ g are exponentiated and hence this means all the contribution to the finite (N -independent) piece from the soft-collinear function is also being exponentiated. The resummation coefficients for the soft exponentiation denoted byg g,Soft 0i and g g,Soft i in Eqs (2.48) and (2.47), respectively can be obtained from thier counterparts in standard N scheme as described below.
The N -independent constantsg g,Soft   where the coefficients ∆ g,Soft g i are given as,

E Resummation Coefficients for the All exponentiation
For the case of All exponentiation, the completeg g 0 is being exponentiated along with the large-N pieces. This brings into modification only for the resummed exponent compared to the 'Standard N exponentiation'. We write the modified resummed exponents denoted by g g,All i in Eq. (2.50) in terms of exponents in standard N as, g g,All 1 =ḡ g 1 , g g,All 2 =ḡ g 2 + a s ∆ g,All g 2 , g g,All 3 =ḡ g 3 + a s ∆ g,All g 3 , (E.1) where ∆ g,All g i terms are found from exponentiating the completeḡ g 0 and they are given as, ∆ g,All g 2 =ḡ g 01 , ∆ g,All where the coefficientsḡ g 0i can be obtained as described in C.