Resummed Higgs cross section at N$^3$LL

We present accurate predictions for the inclusive production of a Higgs boson in proton-proton collisions, via gluon-gluon fusion. Our calculation includes next-to-next-to-leading order (NNLO) corrections in perturbative QCD, as well as the resummation of threshold-enhanced contributions to next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy, with the inclusion of the recently-determined three-loop constant coefficient (sometimes referred to as N$^3$LL' accuracy). Our result correctly accounts for finite top, bottom and charm masses at leading order (LO) and next-to-leading order (NLO), and includes top mass dependence at NNLO. At the resummed level the dependence on top, bottom and charm mass is accounted for at NLL, while only the top mass at NNLL. The all-order calculation is improved by a suitable choice of the soft terms, dictated by analyticity conditions and by the inclusion of subleading corrections of collinear origin, which improve the accuracy of the resummation away from the threshold region. We present results for different collider energies and we study perturbative uncertainties by varying renormalization and factorization scales. We find that, at current LHC energies, the resummation corrects the NNLO result by as much as 20 % at $\mu_R=\mu_F=m_H$, while the correction is much smaller, 5.5 %, at $\mu_R=\mu_F=m_H /\,2$. While the central value of NNLO+N$^3$LL result depends very mildly on the scale choice, we argue that a more reliable estimate of the theoretical uncertainty is found if the perturbative scales are canonically varied about $m_H$.


Introduction
The resummation of soft-gluon (or threshold) logarithms in QCD plays an important role in precision phenomenology at hadron colliders, and in particularly at the LHC. Examples include Higgs boson production in gluon fusion, e.g. [1], top-pair production, e.g. [2] and supersymmetric particles, e.g. [3]. Soft gluon resummation improves the accuracy of the predicted cross section, leading, for instance, to a reduced scale dependence. This is particularly important in the case of Higgs production in gluon-gluon fusion. QCD corrections are fully known up to next-to-next-to-leading order (NNLO) accuracy [4,5] and N 3 LO calculations are underway [6]. Very recently, the first term in the soft expansion of the full N 3 LO cross section has been obtained [7]. The perturbative behavior of this series is very poor and thus logarithmically enhanced soft terms, predicted to all orders by soft-gluon resummation, provide a powerful tool to include (and check) higher order terms in the series.
Strictly speaking, soft-gluon resummation is needed when the partonic subprocess is close to threshold: being M the mass of the tagged final state (the Higgs boson mass, for instance) and √ŝ the center-of-mass energy of the partonic subsystem, in the limit z = M 2 /ŝ → 1 the QCD perturbative expansion of the partonic cross section is unstable and the resummation of the entire series is necessary. Whether this is the case in the computation of the physical hadron-level cross section depends both on hadron-level kinematics and on the shape of parton distribution functions (PDFs), since the physical cross section is a convolution of the partonic cross section and PDFs.
In most cases, and particularly for inclusive observables at the LHC, the partonic region z → 1 for which perturbativity is lost gives only a moderate, often negligible, contribution to the physical cross section [8]. In these cases, soft-gluon resummation is no longer needed ; however, it might still be advisable. Indeed, an intermediate range of values of z for which the soft terms approximate well the full partonic cross section usually exists [8,9]; although in this region the series is behaving in a perturbative way, inclusion of higher order terms from soft-gluon resummation leads to a more accurate and stable prediction for the partonic cross section. When this intermediate range dominates the physical cross section, soft-gluon resummation provides then a powerful way of including (the dominant part of) higher order terms in the perturbative expansion.
The ability of all-order calculations to capture the region of intermediate z strongly depends on the actual form of the soft terms that are being resummed [8,9]. Indeed, while the soft limit determines the large-z (or large-N , being N the conjugate variable of z upon Mellin transformation) behavior of the soft terms, it does not fix their functional form. Traditionally, the N -space resummation of the soft terms is organised in terms of powers of log N and constants. However, using analyticity arguments [9], we arrived at the conclusion that this choice is not optimal in several respects, chiefly because powers of log N exhibit a branch cut at finite N , in contrast to the pole structure of fixed-order coefficient functions. A form of the resummation that respects these analyticity properties is advisable, in particular if we aim to capture the dominant behavior in the region of intermediate N , i.e. intermediate z.
Following our previous studies [9,10], we consider a functional form for the soft terms that respects the analyticity properties of fixed-order results. We improve on that work by implementing this formalism in an all-order resummation formula and by computing the Higgs production cross section at NNLO+N 3 LL, for different collider energies. Our result includes all the information from the N 3 LO soft-virtual calculation of Ref. [7], and reproduces to order α 3 s the soft part of the N 3 LO approximate prediction of Ref. [10]. We also study a different form of soft terms, which has the advantage of respecting the aforementioned analyticity conditions, while having, at the same time, a fast numerical implementation.
Finally, we note that a form of threshold resummation with the correct singularity structure at finite N is a necessary step towards the construction of a double-resummed cross section in which threshold and high-energy (BFKL) logarithms are simultaneously accounted for to all orders.
2 Soft-gluon resummation 2.1 Generalities A physical (hadron-level) inclusive cross section at hadron colliders can be written in the factorized form 2) and i, j run over all parton flavours. Without loss of generality, we can suppress the flavour indices and concentrate on the dominant channels for soft resummation (gg for Higgs). For ease of notation, we also suppress factorization scale µ F and renormalization scale µ R dependence. The partonic cross sectionσ is related to the so-called dimensionless coefficient function C bŷ where σ 0 is the leading order (LO) partonic cross section, so that the coefficient function is normalized to δ(1 − z) at leading order: (2.4) and z = M 2 /ŝ is the variable already mentioned in the introduction. In terms of this coefficient function the cross section Eq. (2.1) reads which has the form of a Mellin convolution, and factorizes in Mellin space Note that we have used the same symbols, with different arguments, for a function and its Mellin transform; note also that, for convenience, we have indicated with σ(N, M 2 ) the Mellin transform of σ(τ, M 2 )/τ . Soft-gluon resummation is generically performed in N -space, where the multiple gluon emission phase-space factorizes. The N -space resummed coefficient function (for Higgs and Drell-Yan production) has the form [11] C res (N, α s ) =ḡ 0 (α s ) expS(α s , N ), where α s = α s (µ 2 R ), andḡ 0 (α s ) does not depend on N , but depends implicitly on µ F /M and µ R /M . The function A(α s ) (also called cusp anomalous dimension Γ cusp ) is the numerator of the divergent part of the relevant 1 diagonal Altarelli-Parisi splitting function, 11) and D(α s ) is a process-dependent function. A given logarithmic accuracy is obtained including the functions A(α s ), D(α s ) andḡ 0 (α s ) up to a given order in Eq. (2.7), according to Table 1. Table 1 shows two notations for the counting of logarithms, usually adopted in different contexts. In particular, in the two notation what is called N k LL for k > 0 represents two different accuracies, and therefore can lead to some confusion. In what we call Notation , the N k LL accuracy without decoration corresponds to a logarithmic counting on log C res , where the inclusion of an additional order inḡ 0 (α s ), corresponding to N k LL , does not increase the formal accuracy. However, as shown Notation* Notation Refs. [1,12] Table 1. Orders of logarithmic approximations and accuracy of the predicted logarithms L = log N . See Ref. [14]. Note that the four-loop contribution to A(αs) is yet unknown and in our study we take a Padé approximation [15].
explicitly in the table, N k LL predicts an additional term to all orders in the tower of logarithms contained in C res . In practice, when the process is not very close to the physical threshold (i.e., τ is not close to 1), as in all relevant cases, the additional logarithm included at N k LL improves the actual accuracy of the result. Therefore, in a wide literature, the N k LL is simply called N k LL, as shown in what we call Notation*, where the lower accuracy is denoted with a *, although this notation is not widespread. In what follows, we will refer to a logarithmic accuracy according to Notation*. The three-loop coefficients of A(α s ) and D(α s ) have been known for while (see for instance Refs. [15,16]), while the O α 3 s contribution toḡ 0 (α s ) has been recently computed in the infinite top-mass limit [7]. The function A(α s ), however, is needed at four loops in order to achieve full N 3 LL accuracy. This contribution is yet unknown; however, a Padé estimate [15] can suggest the size of its value, and a numerical analysis shows that its impact in a resummed result is essentially negligible. Furthermore, we note that, even without this contribution, the expansion of the N 3 LL resummation to third order in the strong coupling completely reproduces all the soft and constant contributions up to N 3 LO, i.e. the four loop coefficient of A(α s ) only enters at order α 4 s .

Large N limit -N -soft
The Mellin transform in Eq. (2.8) is ill defined, because z ranges from 0 to 1, forcing the argument of α s to be arbitrarily small, therefore crossing the Landau pole. However, the integral can be made convergent by using the explicit solution for the running coupling to any finite order. In this way,S can be formally written as 2.15) are the usual plus-distributions and (2.16) their Mellin transform. The Mellin transform in Eq. (2.13) is usually computed, to any finite logarithmic accuracy, in the large-N limit, leading to an expression of the form where we have introduced a new notation (N -soft) for the resummed coefficient function to stress the fact that the large-N limit has been taken. Note that The functions g i , i = 1, 2, 3, 4 can be found explicitly for many processes in Ref. [15]. A given N k LL accuracy is obtained from Eq. (2.17) including g i up to i = k + 1, and g 0 (α s ) up to the same order asḡ 0 (α s ), see Table 1. The function S(α s , log N ) can be written as where the coefficients b n,k are the same as in Eq. (2.13), and the functions D log k (N ) are the large-N limit of D k (N ) expressed in terms of log N , and neglecting constant terms and terms suppressed by powers of 1/N : Here Γ (j) (x) is the j-th derivative of the Euler Gamma function. For completeness, we also report the functional form of the momentum space conjugates of the D log k (N ) functions: The relation between the function g 0 and the constant that multiplies the resummed exponent in Eq. (2.7), namelyḡ 0 , is (for further details, see Ref. [9]).

Prescriptions for the Landau pole
The resummed coefficient function C N -soft , Eq. (2.17), although well defined in N space, cannot be directly used for computing the corresponding hadron-level cross section because its inverse Mellin transform does not exist. Indeed, the functions g i (α s log N ) have a branch-cut for real N > N L = exp 1 2β0αs originating from the Landau pole of the running coupling, while Mellin transformation has a convergence abscissa [17]. On the other hand, if C N -soft is expanded in powers of α s , the inverse Mellin transform exists to any finite order, but the resulting series is divergent [18]. Therefore, a prescription is needed to compute physical observables from Eq. (2.17).
The most used prescription is the so called Minimal Prescription (MP), proposed long ago [17]. It consists on a simple modification of the Mellin inversion integral, and has the advantage of having a fast numerical implementation. More details are given in App. A.1. More recently, another prescription based on a Borel summation of the divergent series of the order-by-order inverse Mellin transform of C N -soft has been first proposed in Ref. [18], and refined and extended in Refs. [19][20][21][22]. This prescription, called Borel Prescription (BP), is typically slower but more flexible. More details are given in App. A.2.
It turns out that the numerical difference between the two prescriptions is small, being totally negligible for Higgs phenomenology [22]. The reason is that, for the kinematic configurations typical of high-energy colliders, i.e. τ 1, the series is behaving in a perturbative way. Hence, the all-order nature of the series does not play any role, and, a fortiori, the way the divergence of the series is dealt with is immaterial. We believe that expanding C N -soft to a sufficiently large, but finite, order in α s and inverting exactly would lead to a result virtually identical to the all-order MP or BP results.
One reason why the effective equivalence of MP and BP does not appear clearly in the literature is the fact that, within the BP, the form of the soft terms that are resummed can be easily modified, and this has been always done for phenomenological application, thereby giving a result which differs from the MP one. A discussion on the form of the soft terms will be performed in Sect. 3. Here, we just want to mention that, for practical applications, the choice of the prescription will be mainly dictated by its flexibility and numerical efficiency.
Details on the practical implementation of the prescriptions, as well as details on a new version of the BP acting directly on the Sudakov exponent S(α s , log N ), are given in App. A.

Soft terms and analyticity conditions
Soft gluon resummation fixes the coefficients b n,k of the soft terms in Eq. (2.22), however it does not fix the functional form of these contributions. In particular, choices that at large N only differ by terms suppressed by factors of 1/N are equally acceptable.
It has been pointed out [9,23] that the actual form of the soft terms is very important, and different choices would lead to very different accuracies, in particular when the considered process is far from the physical threshold.
In our previous analysis [9], we studied the analyticity properties of coefficient functions in N space and we arrived at an optimal choice of the soft terms. We used this improved large-N behavior, together with the knowledge of the rightmost singularity at finite N from high-energy resummation to compute an approximate expression for the N 3 LO Higgs production cross section [9,10].
In this section, we briefly review the two main theoretical ingredients that go into the construction of our resummation formula, namely a choice of soft terms that respects the singularity structure of coefficient functions and the improvement related to the inclusion of collinear contributions.

Functional form of the soft terms
The N -soft resummed exponent Eq. (2.22) is an infinite sum of contributions each of which has a logarithmic branch cut starting at N = 0, which is not compatible with the known singularity structure of coefficient functions. However, this problem is an artefact of the large N approximation which we have employed in going fromS, Eq. (2.13), to S, Eq. (2.22). Indeed, the resummed expo-nentS is written as an infinite sum of D k (N ) functions, whose singularity structure is compatible with the one of fixed-order calculations.
Furthermore, as we discussed at length in Ref. [9], we can improve on the use of D k (N ) by noticing that the correct kinematic limit of the µ 2 integration in Eq.
This consideration leads to the following choice for the soft terms in momentum spacê from which we can easily compute the functionsD k (N ), which enter our resummation formula: Note that in Eq. (3.1) we have chosen to apply the plus prescription only to the first term, singular in z = 1, which is the natural choice in fixed order calculations. In this way,D k (N ) differs from D k (N ) only by terms vanishing at large N : Adopting this form for the soft terms in all the terms generated in Eq. (2.8) (hence also those generated by the D(α s ) term), we arrive at the expression which represents the first improved version of Eq. (2.13). This is not yet our final formula, as we are going to discuss.

Altarelli-Parisi contributions
We have already discussed the origin of the coefficient A(α s ) in the resummation formula Eq. (2.8): it is the coefficient of the soft-enhanced part of the relevant Altarelli-Parisi splitting function Eq. (2.11), where P = P gg in the case of Higgs production. In resummed calculation, the coefficient B(α s ) in Eq. (2.11) is also retained at the appropriate accuracy because it corresponds to a constant term in N -space, while contributions that vanish as z → 1 are usually neglected in Eq. (2.8). However, an important class of these subleading corrections can be taken into account to all orders, essentially because the full leading order anomalous dimension exponentiates [23,24]. However, the 1/z pole present in the LO gluon-gluon splitting function would introduce spurious singularities in the resummed coefficient function at N = 1 [9]. Nevertheless, the expansion of gg (z) being the LO gluon-gluon Altarelli-Parisi splitting function) in powers of 1 − z to any finite order is not singular in z = 0, and therefore does not affect the singularity structure around N = 1. The expansion up to second order reads 5) where A 1 = C A /π is the first coefficient in the expansion of A(α s ), Eq. (2.10). Note that the third order term in the expansion is accidentally zero. Upon Mellin transformation, multiplication by the factor 2 − 3z + 2z 2 results into the shift where we have introduce the Altarelli-Parisi (AP 2 ) operator to second order for future convenience. The A 1 term in Eq. (2.8) is responsible for the tower of LL terms, namely terms α n s D 2n−1 (z) to all orders n in C(z, α s ). When A 1 is replaced with the expansion Eq. (3.5), logarithmic terms suppressed by powers of (1−z) are generated to all orders. The towers of LL suppressed logarithms, namely terms of the form with k running from one up to the order of the expansion of (1 − z)P (0) gg (z), are correctly predicted to all orders [25]. However, the simple inclusion of additional information from Altarelli-Parisi splitting functions is not enough to predict terms beyond these LL towers.
Nevertheless, we have shown in Ref. [9] that including the expansion of (1 − z)P (0) gg (z) up to second order for all the terms in the Sudakov exponent Eq. (2.14) leads to an approximation of the exact fixed-order terms which is very good in a wide range of N values, down to values where high-energy terms (not considered in this work) start being relevant. This is achieved when the soft terms introduced in Sect. 3.1 are used. Therefore, the inclusion of such second order expansion for all the terms in Eq. (2.8) Note that the inclusion of terms of order (1 − z) 4 and higher in the expansion of (1 − z)P (0) gg (z) does not affect the results significantly [9]. At resummed level, the operator AP 2 in Eq. (3.7) can be applied directly to the Sudakov exponent.

A-soft 2
We are now ready to present our resummed expression. We name it 2 A-soft 2 , because it takes into account the analyticity properties discussed in Sect. 3.1 and it includes the first two terms of the (1 − z) expansion of the LO splitting function, as detailed in Sect. 3.2: Note that because the difference between the functionsD k (N ) and D(N ) vanishes at large N , Eq. (3.3), the constant termḡ 0 in Eq. (3.8) coincides with the one in Eq. (2.7). As in the case of Nsoft resummation previously discussed, the series that defines C A-soft2 in Eq. (3.8) is divergent. This series is summed using the Borel Prescription 3 detailed in App. A.2, which leads to the following result , (3.9) where W is a cut-off, minimally set to W = 2. The inverse Mellin transform of the above result is finite, and can then be used to compute hadron-level cross sections. We have noticed in previous papers [9,10] that one of interesting consequences of using A-soft 2 instead of N -soft is a better perturbative behavior, which is partially due to the resulting constant multiplying the resummation,ḡ 0 or g 0 respectively. In fact, while the perturbative expansion of the function g 0 is known to be poor, driving large fixed-order corrections to the Higgs production, the functionḡ 0 has a much more stable perturbative expansion [9,10]. The relation between the two constants is given in Eq. (2.25), from which it is clear that by usingḡ 0 , we are effectively exponentiating and, hence, resumming, part of the constant contribution. This is similar in spirit to the so-called π 2 -resummation discussed in Refs. [26,27] for Higgs and Drell-Yan processes.
We can also go a step further and try to exponentiate the whole constant term [27]. To N 3 LL accuracy this amounts to replacing in the resummed expressionsḡ 0 (α s ) with (3.10) By construction the difference betweenḡ 0 andḠ 0 is O α 4 s and hence beyond N 3 LL accuracy considered here. Because of the good convergence of the perturbative expansion ofḡ 0 we do not expect the result obtained withḠ 0 to be much different compared to the one obtained withḡ 0 . We will further comment on this in Sect. 4, where we present our hadron-level results.

ψ-soft 2
Before turning our attention to phenomenology, we discuss an alternative option for the soft terms, which is very similar to the one used in our improved result A-soft 2 . We first writeD k (N ) in Eq. (3.2) as [22] where ψ 0 (N ) = Γ (N )/Γ(N ) is the DiGamma function. Thus, except for the constant term, the functionsD k (N ) are equivalent to the corresponding D log k (N ) after the replacement log N → ψ 0 (N ), up to corrections of order 1/N 2 . This means that Eq. (2.17) can be upgraded by simply replacing log N with ψ 0 (N ), thereby restoring the analyticity properties of the coefficient function.
Therefore, we propose a new prescription, called ψ-soft 2 , where we also include Altarelli-Parisi improvement, which has the advantage of having almost all the good properties captured by A-soft 2 , while being numerically very fast, because it can be computed using the MP. In fact, any existing soft-gluon resummation code can be easily upgraded to ψ-soft 2 . However, Eq. (3.12) has a clear disadvantage with respect to A-soft 2 , Eq. (3.9), namely the presence of the function g 0 rather thanḡ 0 . We have already commented in Sect. 3.3 about the poor perturbative behavior of g 0 compared to that ofḡ 0 . Here we limit ourselves to mention that, in this case, writing g 0 in exponential form G 0 (α s ) = exp α s g 0,1 + α 2 s g 0,2 −  as done forḡ 0 in Eq. (3.10), can improve significantly the perturbative stability of the resummed result. We anticipate that a phenomenological study shows that, after exponentiation of both g 0 andḡ 0 , results obtained with A-soft 2 and ψ-soft 2 are very similar.

Hadron-level results
In this section we present numerical results for the Higgs production cross section in the gluon fusion channel at the LHC. Our result correctly accounts for finite top mass (m t = 172.5 GeV) at NNLO [5] and NNLL. We also include bottom and charm masses in the LO prefactor (thus changing the overall normalization) and in the NLO contributions, using the results of Ref. [28]. The dependence on bottom and charm masses in the resummation is NLL accurate. For consistency with the NNLO set of parton distribution functions NNPDF2.3 [29] with α s (m Z ) = 0.118, that we adopt for our phenomenological analysis, we use m b = 4.75 GeV and m c = 1.41 GeV. We refer the reader to Ref. [30] for a discussion about the role of N 3 LO parton densities. We use m H = 125 GeV.

A study of the scale dependence
In order to make contact with our previous work [9], we start by showing in Fig. 1 the cross section as a function of µ R at LO, NLO, and NNLO with only top included in the loop (left panel) and with also bottom and charm (right panel). The results have been computed using the code ggHiggs. The collider energy is √ s = 8 TeV. We use different colors for different values of µ F in all curves, for µ F = {2, 1, 1/2, 1/4}m H . We show four choices because the typical scale variation is by a factor of two about its central value, but the central value is sometimes suggested to be m H (e.g. [1]) and sometimes m H /2 (e.g. [31]). It is interesting to observe that fixed-order results with only the top quark in the loop, left-hand plot of Fig. 1, barely depend on µ F : this is due to a compensation between different channels (mainly gg and qg channels) [9]. However, this cancellation is not as perfect at NNLO, when we introduce bottom and charm contributions, right-hand plot of Fig. 1. This happens because in our framework, these corrections are correctly implemented only at NLO and there are no m b , m c dependent contributions at O α 2 s to compensate the NLO µ F dependence. We observe that the main effect of including bottom and charm in the loop is to significantly reduce the cross section.   We now move to resummation. In order to study the effect of different logarithmic orders, we show in Fig. 2 the resummation at LL, NLL, NNLL and N 3 LL accuracy 4 , always matched to the same NNLO contribution, as a function of µ R , for fixed µ F = m H . We also show, for comparison, LO, NLO and NNLO curves. The fixed order results have been computed using the code ggHiggs, while for the resummation we have written a new code called ResHiggs. The plots show our best prediction, A-soft 2 , withḡ 0 (left panel) and its exponentiated versionḠ 0 (right panel). It is interesting to observe that exponentiatingḡ 0 leads to a flatter resummed result, thereby suggesting that its exponentiation is probably improving the convergence of the series. We also observe that, in any case, the N 3 LL result is very similar in both cases over a wide range of scales, so the exponentiation ofḡ 0 does not change significantly the final result, as we have anticipated at the end of Sect. 3.3. In both cases, we note that the inclusion of soft-gluon resummation at N 3 LL significantly reduces the µ R scale uncertainty of fixed-order results and of previous resummed orders.
In Fig. 3 we concentrate on NNLO+N 3 LL and also show the effect of varying µ F . Since the resummation involves only the gg channel, the resummed result depends more significantly on the scale µ F , although formally such dependence is of order α 3 s with respect to the Born cross section. Over a range of roughly a factor of 2 about µ R = m H /2 the results with (right panel) or without (left panel) exponentiation ofḡ 0 are very similar, while they differ (and are more sensitive to µ F ) for more extreme choices of µ R (especially at small µ R ). In these regions, the result obtained exponentiatinḡ g 0 looks more sensible and stable, suggesting, once again, that exponentiatingḡ 0 provides a more stable result. Moreover, we notice that NNLO+N 3 LL result with µ F = m H / 2 barely depends on µ R . We also observe that resummed curves for different values of µ F approximately coincide for a value of µ R slightly smaller than m H /2.
In Fig. 4 we show the same plots as in Fig. 3, but this time obtained with the ψ-soft 2 prescription. Since now the constant function in front of the exponential is g 0 rather thanḡ 0 , we can expect a result different from that of A-soft 2 , when g 0 is not exponentiated (left panel). However, the result with G 0 (right panel) is very similar to the analogous result with A-soft 2 . It follows that ψ-soft 2 provides an acceptable alternative to our best choice A-soft 2 , provided that G 0 is used, i.e with g 0     exponentiated. Since the numerical implementation of ψ-soft 2 is much faster than that of A-soft 2 , its usage can be convenient. We show the result of the more traditional N -soft resummation in Fig. 5. This is interesting because the value of the Higgs production cross section which is currently recommended by the Higgs Cross Section Working Group [32] is based on Refs. [12], which includes N -soft resummation to NNLL accuracy. In this case the difference between N 3 LL and NNLL is bigger than it was for A-soft 2 , as it is shown on the left-hand plot of Fig. 5. This difference is partly due to the fact that g 0,3 is much larger thanḡ 0,3 [9,10]. The dependence of renormalization and factorization scale at N 3 LL is comparable ψ-soft 2 when g 0 is not exponentiated; however, the typical increase of the cross section is smaller than what we find with our result.   Figure 6. Ratios of different resummed results to our best prediction A-soft2 with the exponentiated constantḠ0, plotted as a function of µR, for different choices of µF.
A quantitative comparison between the different resummed results is shown in Fig. 6, where ratios to our best prediction, namely A-soft 2 with the exponentiated constantḠ 0 , are plotted as a function of µ R , for different choices of µ F . As previously observed, we confirm here quantitatively that the result obtained with ψ-soft 2 with g 0 exponentiated (solid red line) is almost identical to our best prediction, the difference being always below 1%, and confirming that this prescription can be indeed used as a numerically convenient alternative to A-soft 2 withḠ 0 . We also observe that for a wide choice of scales not exponentiatingḡ 0 in A-soft 2 (dashed black line) leads to a result which only differs from the result withḠ 0 by a few percent. In contrast, the difference between resummed results with g 0 , e.g. dotted red curve, or its exponentiated version G 0 , e.g. solid red curve, is more pronounced. Finally, we note that the difference between our best prediction and N -soft is about 10% for µ R = µ F = m H . In Figs. 7 and 8, we show the results for A-soft 2 and ψ-soft 2 with collider energy √ s = 13 TeV. The shape of the NNLO+N 3 LL cross section in terms of renormalization and factorization scales, is very similar to the ones we have found at 8 TeV, both for A-soft 2 and ψ-soft 2 , and the analogous ratios showed in Fig. 6 for 8 TeV look almost identical at 13 TeV. Thus, most of the points we have discussed for the 8 TeV case, also apply at 13 TeV. We just note that the NLO cross section exhibits a somewhat larger µ F dependence at 13 TeV than at 8 TeV.

Numerical results for different collider energies
In Table 2 we summarise the results for our best prediction, namely A-soft 2 withḠ 0 , for different collider energies and different values of the Higgs mass. We estimate the theoretical uncertainty by independently varying the scales up and down, by a factor of two m x /2 < µ R , µ F < 2m x , with the  Table 2. Values of the NNLO and NNLO+N 3 LL (A-soft2 withḠ0) gluon fusion cross section for selected values of the collider energy. We use NNPDF2.3 [29] with mH = 125 GeV, mt = 172.5 GeV, mb = 4.75 GeV and mc = 1.41 GeV. The central value is for µR = µF = mH (on the left) and µR = µF = mH/ 2 (on the right). As detailed in the text, we recommend to evaluate theoretical uncertainties by scale variation around µR = µF = mH. Electro-weak corrections are not included.
where m x = m H (on the left) and m x = m H /2 (on the right). Several comments can be made on the results in Table 2. We first discuss on the effect of the resummation with respect to NNLO. If µ R = µ F = m H is chosen as the central scale, we find that the inclusion of the resummation increases the NNLO prediction by 22% at √ s = 7, 8 TeV and 21% at √ s = 13, 14 TeV. This is consistent with the 16% increase we previously found just for N 3 LO approx [9,10]. We remind the reader that the absolute value of the cross sections reported in Table 2 is lower than the one of our previous studies because we now include bottom and charm quarks in the loops. On the other hand, if µ R = µ F = m H / 2 is chosen as the central scale, we find that the resummation only corrects the NNLO result by roughly 5.5% at √ s = 7, 8 TeV and √ s = 13, 14 TeV. It is important to note that, while the NNLO cross sections at µ R = µ F = m H and µ R = µ F = m H / 2 differ by more than 10%, the NNLO+N 3 LL is much more stable and only varies by 3% (in the opposite direction with respect to the change of the fixed-order). Thus, as expected, the resummation of soft-enhanced contributions reduces the theoretical uncertainty related to the choice of the hard scale.
However, while the NNLO+N 3 LL central values in Table 2 are rather similar, their uncertainties differ significantly. The total uncertainty band 5 at µ R = µ F = m H is around 12-13%, while at µ R = µ F = m H / 2 is 4-4.5%. This dramatic reduction of the scale dependence can be understood thanks to the study we have performed in Sect. 4.1, where we have noticed, for instance, that A-soft 2 (and ψ-soft 2 ) curves with µ F = m H / 2 barely depend on µ R . However, in Sect. 4.1 we have also noticed that resummed curves for different µ F all meet to a point which is not too far away from µ R = m H / 2, resulting into an artificially small scale dependence, which is probably not representative of the true theoretical uncertainty.
Therefore, in order not to underestimate the uncertainty of the NNLO+N 3 LL result, we recommend to vary the scales around µ R = µ F = m H , which results into a fairly conservative 12-13% uncertainty band, yet significantly smaller than the 20-22% band of the NNLO result. One can also imagine implementing less canonical scale choices. For instance, the plots in Sect. 4.1 seem to suggest that varying the scales around µ F = m H /2 and µ R = m H would lead to a less conservative, but still reliable, estimate of the theoretical uncertainty.
We also compute the cross section for a future, very energetic run of the LHC, with √ s = 33 TeV 6 . The relations between fixed-order and resummed results, and their uncertainty bands, are qualitatively similar to the ones found at lower energies. However, we note that at √ s = 33 TeV we may be becoming sensitive to high-energy (BFKL) corrections, that we are not resumming here. These corrections have negligible impact on the inclusive Higgs cross section at current LHC energies, but can play a significant role for future high-energy runs.
Our analysis only includes QCD corrections. Electro-Weak contributions [33] can be also taken into account at NLO and if one assumes that they factorize from QCD corrections, they typically lead to a few percent increase of the cross section. Moreover, our calculations assume that the Higgs boson is produced on its mass shell. However, off-shell effects can be relevant for precision Higgs physics, even if the Higgs is light [34]. We leave the inclusion of these effects, as well as a study of PDF uncertainties, to future work.
Finally, it is interesting to compare our result to the lower-order prediction NNLO+NNLL obtained with N -soft, which is the accuracy of the QCD calculation of the cross section currently recommended by the Higgs Cross Section Working Group [32], for the central scale µ F = µ R = m H . At 8 TeV, for instance, the cross section we obtain 7 at this accuracy is 18.57 +2.06 −1.81 pb, with a total uncertainty band of 15%. Hence, our calculation leads to a 14% increase with respect to what is currently used by the LHC experiments, and reduces the theory uncertainty.

Conclusions
Accurate theoretical predictions for Higgs production in gluon-gluon fusion are of primary importance for the LHC physics programme. In this paper we have studied the all-order resummation of enhanced contributions due to the emissions of soft gluons.
Our resummation (A-soft 2 ) differs from the traditional one (called here N -soft), in that it is build in such a way to respect the analyticity properties of coefficient functions [9,10]. We improve on our previous work by computing all-order N 3 LL matched to NNLO cross sections for proton-proton collisions at different center-of-mass energies. Our result therefore accounts for all the contributions also present in the N 3 LO soft-virtual approximation of Ref. [7]. Matching to full N 3 LO, when this becomes available, will be straightforward.
Because of its analyticity properties, our results can be also easily matched to the resummation of high-energy (BFKL) contributions, which we leave for future work. However, high-energy resummation by itself has a relatively small direct numerical impact, at current LHC energies, so we do not expect this inclusion to significantly alter the results presented in this paper. This situation can change, if future, very high-energy colliders are considered.
We have also shown that the traditional resummation formula can be minimally modified (ψsoft 2 ) in such a way that the final result is almost identical to A-soft 2 , when constants are exponentiated in both cases. This provides a fast alternative to implement our improved soft-gluon resummation results.
The resummation presented here has been implemented in a code called ResHiggs which can be interfaced with the fixed-order code ggHiggs. Both codes are publicly available at the website http://www.ge.infn.it/∼bonvini/higgs/.
We have shown that, at energies relevant for LHC Run-I and Run-II, the resummation (Asoft 2 ) corrects the NNLO result by as much as 20% at µ R = µ F = m H , while the correction is much smaller, 5.5%, at µ R = µ F = m H / 2. However, the central value of NNLO+N 3 LL results depends very mildly on the scale choice. Moreover, the result obtained at NNLO+N 3 LL (A-soft 2 ) at the permille level. 7 The number presented here slightly differs to the one in Ref. [32]. Most of the difference comes from the presence of Electro-Weak contributions in the result of Ref. [12], which leads to an increase of about 5% [33]. The remaining difference (about 1%) is due to the different PDF set used and the lack of finite top mass correction at NNLO in Ref. [12]. leads to a 14% increase with respect to the N -soft result at NNLO+NNLL [12] at µ R = µ F = m H , corresponding to the accuracy recommended by the Higgs Cross Section Working Group [32] and currently used by the LHC experiments. We have also argued that theoretical uncertainties are better estimated by scale variations about µ R = µ F = m H , with a resulting 12-13% total band.
Our results have a moderate dependence on the renormalization scale, which instead drives the theoretical uncertainty of fixed-order calculations. The dependence on the factorization scale instead is somewhat larger than in fixed-order calculations. This is mainly due to the fact that we only resum the gg channel. This situation is likely to improve, if we were able to resum at least the leading logarithms in the qg channel, which are of the form α k s log 2k−1 (1 − z), using techniques similar to the ones developed in Ref. [35] for deep-inelastic scattering or by resumming the whole class of the so-called next-to-eikonal contributions [36].
converges numerically only for Re N > 0, while the MPc contour involves values of N with negative real part. Therefore, for an efficient practical realization of the MP, the Mellin transform Eq. (A.3) must be computed analytically, meaning that we need a functional form for L (x).
A very efficient way of approximating L (x) is obtained expanding on a basis of Chebyshev polynomials the function on the range u min < u < 0, with u min ≤ log τ , and where β and γ are parameters aimed to make f (u) as flat as possible (by default we use β = 1 and γ = 0). Fast routines are available for computing the coefficient of the expansion of f (u) on a basis of Chebyshev polynomials. After straightforward manipulations [22] we are able to write, for integer γ, where n is the order of the Chebyshev approximation and c k are coefficients which depend on µ F . Its Mellin transform is suggesting that a small value for γ is advisable to reduce the number of terms in the sum. This luminosity can then be used in Eq. (A.2) for numerical evaluation. The whole procedure leads to a very fast numerical implementation: therefore, where the MP is applicable, its usage is preferred.

A.2 Borel Prescription
A prescription that deals directly with the divergent nature of the series of the order-by-order Mellin inversion of C N -soft has been proposed in Refs. [18][19][20][21][22]. This prescription adopts a Borel method for summing the divergent series, and it is therefore called Borel Prescription (BP). We briefly review here the derivation, while referring to the original works for a detailed discussion. The k-th derivative of a function can be written as where in the first line we have used the Cauchy formula (the integration encloses the singular point ξ = 0) and in the second line we have rewritten the k! as an integral. As a result, the k-th derivative of f has been translated into an integral operator acting on the k-th power of the variable w/ξ. This operator can be used to translate a power of log N in Eq. (2.17) into so that we can rewrite Notice that the function C N -soft (e −w/ξ , α s ) has a cut in −2β 0 α s w < ξ < 0, due to the Landau pole. Therefore, the ξ integration contour, which must encircle the cut, extends to minus infinity for w → ∞, where the oscillatory behavior of the integrand makes the integral divergent. This shows that the series is not Borel summable. The Borel Prescription is now formulated as [18][19][20][21][22] C N -soft (N, α s ) BP = 1 2πi W 2β 0 αs 0 dw e −w dξ ξ N −ξ C N -soft (e −w/ξ , α s ), (A.10) where the w integral has been cut off, and W is minimally W = 2 (corresponding to the inclusion of twist 4 terms). This cut-off makes the integral convergent, because the cut extends always on a finite range. Since the BP is applied directly to the coefficient function C N -soft , the physical cross section maintain its convolution structure, without any violation of factorization. As stressed already in Sect. 2.3, the numerical result of Eq. (A.10) is virtually identical to that obtained with the MP Eq. (A.2), since for phenomenologically relevant kinematic configurations the series is behaving perturbatively and the way the divergence of the series is tamed is immaterial [22].
One of the advantages of the BP is that the inverse Mellin transform can be computed analytically, since the N dependence is confined in the generating function N −ξ . The convolution with the parton luminosity can hence be computed directly in z space, without the need of any approximation (although in practice approximations might be convenient numerically, see Ref. [22]). The price to pay is a slower implementation, since the BP consists of two integrals, while the MP needs just one integration.
Another advantage, which is the main reason to consider the BP, is the fact that the form of the soft terms is completely under control. Indeed, in the derivation of Eq. (2.17) from Eq. (2.7), a large-N limit of the Mellin transform of the soft terms Eq. (2.15) is taken, in order to transform the complicated derivatives of Γ functions into simple powers of log N , so that a closed form for the functions g i (α s log N ) can be found. The BP makes possible to skip this large-N limit, since the derivatives can be translated into powers through the operator Eq. (A.7). This can be obtained by simply replacing the generating function N −ξ in Eq. (A.10) with the generating function of the derivatives in Eq. (3.2), which corresponds to introducing soft terms with the right analyticity properties (see discussion there).
Moreover, in deriving Eq. (2.17) all constant terms from the Mellin transform of plus distributions have been (by choice) removed from the exponent and put into the function g 0 (α s ), Eq. (2.25).
Here we have the opportunity to also restore these constants in the exponent, since they were there in the original expression Eq. (2.7).
However, the proper way to restore the original logarithms and constants is to let the BP operator act directly on the exponent S(α s , log N ). Therefore, we propose a new version of the BP, S α s , − w ξ , (A.11) which can be further improved by letting the Altarelli-Parisi operator, Eq. (3.6), act on the generating function, as in Eq. (3.9). Since the conversion of soft terms is done at the exponent, the two integrals of the BP have to be computed for each value of N , and then the inverse Mellin transform has to be computed numerically. In practice, we can use the same numerical technique used for the MP (an integral of the form Eq. (A.2), with the luminosity approximated as Eq. (A.6)), with the exception that now the function Eq. (A.11) does no longer have a cut, and hence the inverse Mellin transform exists. The resulting numerical implementation is slow, compared to the MP. However, this procedure is, to our knowledge, the only viable solution to reproduce the structure of Eq. (2.7) to all orders in a numerical way.