On the Higgs cross section at N$^3$LO+N$^3$LL and its uncertainty

We consider the inclusive production of a Higgs boson in gluon-fusion and we study the impact of threshold resummation at next-to-next-to-next-to-leading logarithmic accuracy (N$^3$LL) on the recently computed fixed-order prediction at next-to-next-to-next-to-leading order (N$^3$LO). We propose a conservative, yet robust way of estimating the perturbative uncertainty from missing higher (fixed- or logarithmic-) orders. We compare our results with two other different methods of estimating the uncertainty from missing higher orders: the Cacciari-Houdeau Bayesian approach to theory errors, and the use of algorithms to accelerate the convergence of the perturbative series. We confirm that the best convergence happens at $\mu_R=\mu_F=m_H\,/\,2$, and we conclude that a reliable estimate of the uncertainty from missing higher orders on the Higgs cross section at 13 TeV is approximately $\pm4$%.


Introduction
The major success of the first run of the CERN Large Hadron Collider (LHC) was the discovery of the Higgs boson [1,2]. Run-I studies of this new resonance established its mass and quantum numbers, as well as its coupling strengths to many Standard Model particles [3,4]. Moreover, measurements of fiducial cross section and distributions have also been performed [5][6][7][8]. The start of the LHC Run II marked the beginning of the Higgs precision era. Now is the time to study the properties of this new particle in detail so that, nearly fifty years after its proposal, the Brout-Englert-Higgs mechanism for electroweak symmetry breaking can undergo its final tests. One key element in the LHC Run-II rich physics program is the measurement of the Higgs production cross section and differential distributions, and their comparison to state-of-the-art theoretical predictions. Thus, identifying the different sources of uncertainties and performing more refined calculations to improve on them is of paramount importance. In this paper we consider inclusive production of a Higgs boson in hadronhadron collisions, which is chiefly driven by the gluon-fusion mechanism. QCD corrections to Higgs production in gluon fusion are known to be very important, leading to poor convergence properties of the perturbative series. Indeed, it is known that the next-to-leading (NLO) contribution corrects the Born cross section by roughly 100%. Moreover, NNLO corrections are also significant. Recently, a milestone calculation of third-order perturbative coefficient was performed [9][10][11][12]. N 3 LO corrections appear to be small, indicating that the perturbative series is finally manifesting convergence.
Threshold resummation is a reorganization of the perturbative expansion that accounts for logarithmically enhanced contributions to all-orders in the strong coupling α s . In the same way as fixed-order calculations can be performed at different (N k LO) accuracy, resummed calculations can be systematically improved by including next k -to-leading logarithmic (N k LL) corrections. Threshold resummation for Higgs production is currently known to N 3 LL [13][14][15][16]. In this work we update the results of Refs. [13], in view of the recently computed N 3 LO result [9][10][11][12]. We discuss the impact of the resummation on the central value and we suggest a robust way of estimating the theoretical uncertainty due to missing higher orders, which goes beyond traditional scale variation.
The paper is organized as follows. We start in Sect. 2 by reviewing the results of Refs. [9][10][11][12] and using them to assess the validity of the fixed-order approximation based on analyticity properties of perturbative coefficient functions of Refs. [17,18]. In Sect. 3, we then review basic results in threshold resummation as well as describe the improvements suggested in Ref. [13]. Sect. 4 contains the main results of this study, namely the matched N 3 LO+N 3 LL cross section and a detailed study of the perturbative uncertainty. Finally, in Sect. 5 we compare our findings to different methods to assess the size of missing higher orders, namely the Cacciari-Houdeau method [19][20][21] and the application of convergence acceleration algorithms to the perturbative series, as suggested by David and Passarino [22].

The N 3 LO cross section
In order to fix the notation, we write the hadron-level cross section for the production of a Higgs boson with mass m H in hadron-hadron collisions as where L ij (z, µ 2 ) is a parton luminosity and i, j run over all parton flavours. In the dominant gluon-fusion production mechanism, the Higgs is generated by the fusion of two gluons through a quark loop. For simplicity, we have only considered a top quark with mass m t running in the loop. Moreover, for ease of notation, the dependence on factorization scale µ F and renormalization scale µ R is often left understood. The partonic cross sectionσ ij is then related to the dimensionless coefficient function C ij bŷ σ ij (z, m 2 H , m 2 t , α s ) = z σ 0 (m 2 H , m 2 t ) C ij (z, m 2 H , m 2 t , α s ). (2.3) where σ 0 is such that the coefficient function has the perturbative expansion in the strong coupling α s While the NLO coefficient C (1) ij is known exactly [23] and C (2) ij is known as an expansion in m 2 H /m 2 t matched to the (full-theory) small-z limit [24][25][26][27][28][29], the third order coefficient has been computed only recently in the large-m t effective theory [9][10][11][12].
In the large-m t effective field theory (EFT), the heavy top is integrated out and consequently the dependence on the top mass only appears in a Wilson coefficient squared W . The EFT is usually improved with a rescaling of the cross section by the ratio of the exact LO over the LO in the EFT, leading to the so-called rescaled effective theory (rEFT), which is known to be a very good approximation for m H m t and for not too high collider energies, essentially because the large-s, i.e. small-z, behaviour of the EFT coefficient functions is double-logarithmic [30], while the full theory only exhibits single high-energy logarithms [24]. In the EFT, the coefficient function further factorizes whereC ij has an expansion analogous to C ij , Eq. (2.4), and W = 1 + O(α s ). Thus, in the EFT, the coefficient C (3) ij is given by where W (k) are the coefficients of the expansion of the Wilson coefficient squared W in powers of α s . The lower order coefficientsC (1) ij andC (2) ij are fully known, whileC ij is the recently computed N 3 LO contribution published in Ref. [12].
For its computation, the third order coefficientC (3) ij has been decomposed into contributions proportional to ln k (1 − z) with k = 0, . . . , 5. For k = 3, 4 and 5 the exact result is known [10], while for lower powers 0, 1 and 2, the result has been expressed in terms of a soft expansion in 1 − z up to order (1 − z) 37 . 1 While not included in the result of Ref. [10], we stress that the leading small-z logarithm, 1 z ln 5 z, is known from small-z resummation [30] (see App. A.1 for details). The construction of C (3) ij , Eq. (2.6), has then some degree of arbitrariness. The formal accuracy is set by the soft expansion ofC (3) ij , so in principle one could take a soft expansion of all its ingredients, in particular the lower ordersC (1) ij andC (2) ij . Alternatively, one can retain as much available information as possible, thereby using the exact expressions forC (1) ij andC (2) ij and adding the leading smallz logarithm. These (and any other intermediate) options are all perfectly valid and the difference among them could be considered a measure of the uncertainty related to the soft expansion. We have found that the residual uncertainty on the soft expansion of C (3) ij mostly comes from the small-z logarithms. Thus, for all phenomenological applications where the small-z logarithms are negligible, the soft expanded result for C (3) ij (in any practical incarnation) is perfectly reliable and accurate, in full agreement with the analysis of Ref. [10]. On the other hand, the uncertainty related to the soft expansion becomes larger as we increase the centre-of-mass energy, leading to sizeable effects at energy scales of Future Circular Colliders (FCC), √ s 100 TeV. Furthermore, even if we controlled the full z dependence of C (3) ij , the EFT result itself becomes unreliable at very large energies because of the appearance of double logarithmic contributions, as previously discussed. Therefore, rather than improving the EFT result with EFT small-z logarithms, one can try to improve it with finite m t effects. One way to do it is including those m t -dependent contributions which are predicted by all-order calculations in the threshold and high-energy limits. As we will see explicitly in Sect. 3, the m t dependence of the soft logarithms appears in a factorized form, and it is therefore possible to account for the exact m t dependence of all the logarithmic terms, i.e. all the plus-distributions, that appear at N 3 LO. The situation in the opposite, small-z, limit is less satisfying because the resummation is known strictly speaking only at the first non-trivial order [17,24], although a large class of subleading running-coupling corrections can also be resummed [31,32]. Both effects have been implemented 2 in the code ggHiggs, version 3.1, publicly available from the webpage [33]. Further details on the ggHiggs implementation are given in App. A.1. In order to illustrate our findings, we now focus on the dominant gg channel, and consider the "n-th order K-factor" K (n) gg , defined as the contribution to the K-factor coming from C (n) gg only, namely such that In Fig. 1 we show the results for the NLO, NNLO and N 3 LO K-factors as a function of the collider energy √ s, for fixed Higgs mass m H = 125 GeV, and top pole mass m t = 172.5 GeV. We use the PDF4LHC15_nnlo_100 PDF set [34][35][36][37][38]. Results obtained in the large-m t EFT are shown in dashed blue, while solid black curves also contain finite-m t corrections (we remind the Reader that the m t dependence is fully accounted for at NLO, while it is treated as a power expansion at NNLO). In all cases, differences between large-m t and finite-m t is small, but it increases with √ s, as expected. At N 3 LO, we computed the EFT result using the exact expressions forC (1) ij andC (2) ij in Eq. (2.6), whileC (3) ij is soft-expanded. 3 The "m t -improved EFT" curve (black solid) is constructed as described above: in particular, the coefficients of all the large-z logarithms have the correct m t -dependence, while in the opposite, small-z, limit we have included the appropriate leading logarithmic term 1 z ln 2 z. In all plots of Fig. 1 we additionally show the approximate result of Ref. [17,18] in red, with its uncertainty, in order to assess the validity of that approach. The approximation is based on a combination of the soft and high-energy behaviours, including finite m t effects. The soft part of the approximation, which gives the dominant contribution, presents several improvements with respect to standard soft approximations (see e.g. [39]); in particular, two versions of the soft approximation, denoted soft 1 and soft 2 , have been considered, and the area between the two has been considered as the uncertainty on the soft result, using the average as central prediction. This improved soft definition has been later extended to all orders in Ref. [13]; in that work, two of us noticed that soft 2 (denoted A-soft 2 there) performed much better than (A-)soft 1 . Since then, we decided to centre the approximation on soft 2 , symmetrizing the difference between soft 1 and soft 2 about soft 2 to get the uncertainty (which is then twice as large as in the original version). This is how we now compute the (soft part of the) error in Fig. 1. This prescription was also used for the prediction published in Ref. [7].
As known from Ref. [17], the approximation is very accurate at NLO and NNLO, the exact result lying well within the uncertainty band, as shown in the first two panels of Fig. 1. In the third panel, we show for the first time the approximate prediction of the N 3 LO compared with the full EFT result of Ref. [12] (in dashed blue) and to the m t -improved EFT (in solid black).
As expected, at lower collider energies, where the process is closer to threshold, the approximation is very good and perfectly compatible with the exact result within uncertainties. Note that at these energy scales the missing terms beyond the soft expansion at N 3 LO are completely negligible, so the full N 3 LO result can be regarded as exact. At very high (FCC) energies, the comparison between EFT (blue) and approximation (red) is not significant. Perhaps surprisingly, we also find that our approximation and the m t -improved EFT prediction differ at very high energy (although they are still compatible within errors) despite having the same leading logarithmic behaviour. The difference is due to subleading high-energy terms 1 z ln z and 1 z , which are included in the approximate result because they are driven by the resummation of the splitting functions [17], but not in the m t -improved EFT because they are not fully determined. The effect of these subleading terms is rather large (larger than at NNLO) and a full determination of at least the 1 z ln z contribution is highly desirable. Finally, we see that in the intermediate region, relevant for LHC, our approximation works less well and the full result lies at the lower edge of the uncertainty band of the approximate result. This seems to suggest that, at N 3 LO, the contributions which are neither soft nor high-energy are more important than at previous orders, a fact which was not taken into account when estimating the uncertainty from these terms in Ref. [17]. Nevertheless, in the gg channel considered here, the overall agreement of the approximate result with the full result remains rather good.

Threshold resummation(s)
Soft-gluon resummation is usually performed in Mellin (N ) space, where the multiple gluon emission phase-space factorizes. The cross section Eq. (2.1) has a simpler structure in Mellin space: where we have used the same symbols, with different arguments, for a function and its Mellin transform. Note that threshold resummation only affects the gg channel: we therefore suppress the flavour indices and implicitly focus on the gg channel. We will later comment on the role of the quark channels. The N -space resummed coefficient function has the form (see [13] and reference therein): whereḡ 0 (α s , µ 2 F ) does not depend on N . We note that in the full theory, all the top-mass dependence is inḡ 0 . Furthermore, under the rEFT assumption, its expression factorizes as where nowg 0 (α s , µ 2 F ) does not depend on the top mass. Note that we have restored explicit scale dependence and we have chosen the factorization scale µ F as the scale of the running coupling α s = α s (µ 2 F ). The three-loop coefficients of A(α s ) and D(α s ) have been known for a while (see for instance Refs. [39][40][41]), while the O α 3 s contribution tog 0 has been recently computed [9]. The four-loop contribution to A(α s ), which is needed to achieve full N 3 LL accuracy, is unknown. However, a Padé estimate [40] can suggest the size of its value, and a numerical analysis shows that its impact in a resummed result is essentially negligible.
The integrals in Eq. (3.3) can be computed at any finite logarithmic accuracy by using the explicit solution of the running coupling, in terms of α s at a given reference scale, which we can also choose to be µ F in first place. At this point we have a result which depends on a single scale µ F , with α s always computed at µ F (note that, while the µ F dependence ofS is explicit, the one ofḡ 0 can be recovered by imposing µ F -independence of the full cross section). In order to write the result in a canonical way, we further evolve α s from µ F to µ R using the explicit solution of the running coupling equation at sufficiently high order, and propagating the resulting logarithms in the various terms at each fixed-order (inḡ 0 ) and logarithmic-order (inS) accuracy. Then, the final result explicitly depends on both µ R and µ F .
The computation of the integrals in Eq. (3.3) is rather cumbersome when performed exactly. The resulting expression was called A-soft in Ref. [13]. The computation is much simpler when performed in the large-N limit, where the result of the integrals is written as a function of ln N only. We call the result in this limit N -soft. Explicit expressions forS in the N -soft limit up to N 3 LL are given in Ref. [40] 4 with full µ F and µ R dependence.
In Ref. [13] two of us proposed a variant of the N -soft resummation based on the simple replacement

7)
ψ 0 (N ) being the Euler digamma function. This prescription was called ψ-soft. It has the advantage that it reproduces the functionS up to corrections of O(1/N 2 ), while N -soft reproducesS only up to O(1/N ) corrections. Note that this does not mean that ψ-soft resummation is accurate at next-to-soft (NS) O(1/N ) level, because the original expression Eq. (3.2) was not. However, as pointed out in Ref. [13], a class of NS terms can be predicted to all orders by adding a "collinear improvement" to Eq. (3.3). This is achieved by recalling that the function A(α s (µ 2 ))/(1−z) is nothing but the divergent part (in the z → 1 limit) of the Altarelli-Parisi splitting function P gg (z, α s ). One can therefore keep more terms in the soft-expansion of P gg about z = 1. As shown in Ref. [13], by using the LO gluon splitting function up to order (1 − z) k−1 one accounts for the leading-logarithmic (LL) N k S terms correctly to all orders. In Ref. [13] two variants of the collinear improvement, AP1 and AP2, were considered, obtained by expanding the LO P gg to first and second order in 1 − z, respectively. Since these corrections correspond to extra powers of z, the effect is to shift the value of N . Extending the shifts to the wholē S function, we therefore have 5 We have verified that AP2 leads to more reliable results, in the sense that expanding the ψ-soft resummed expression with AP2 to fixed-order reproduces to a good accuracy the exact results (see Ref. [17,42], and discussion in Sect. 2 which extends the validation to the third order). In fact, comparing AP2 with AP1 also allows to estimate the uncertainty due to missing 1/N terms, and constructing an uncertainty from the difference between AP2 and AP1 was indeed successful (see again Ref. [17] and Sect. 2). We finally turn to discussing another source of uncertainty at resummed level coming from subleading logarithmic terms. The functionS, Eq. (3.3), contains on top of logarithmically enhanced contributions (terms which grow logarithmically at large N ) also constant (N -independent) terms, analogous to those included inḡ 0 . In standard N -soft resummation (see e.g. Refs. [14,43,44]) all the constants are usually removed from the exponent and collected in the function in front, which then changes name to g 0 (without bar; consequently alsoS changes into a new function S which contains only logarithmic terms). Alternatively, all constants can be moved into the exponent [13], where lnḡ 0 is meant to be expanded to the appropriate order (which is O(α 3 s ) for N 3 LL accuracy). 6 Up to the working logarithmic accuracy, the position of the constants does not make any difference. However, beyond the working logarithmic accuracy, moving constants produces, by interference, different subleading term. Therefore, one can consider Eq. (3.10) and Eq. (3.11) as two opposite options which treat in a maximally different way these subleading terms, and use them to assign an uncertainty to the default (most natural) expression Eq. (3.2). Note that, since constants are known to play an important role for Higgs production [45][46][47], these variations provide a robust way to estimate the perturbative uncertainty.
Having described the various prescriptions available for the threshold resummation, we now move to a description of how we propose to use them to improve the central value and most importantly the uncertainty from missing higher (fixed-or logarithmic-) orders for the inclusive Higgs cross section.
First, we discuss how threshold resummation is matched to a fixed-order calculation. The coefficient function C res (N, α s ) contains all orders in α s but it is accurate only in the soft limit. Assuming we have available the exact result for the coefficient function up to O(α k s ), to maximally use information from both exact fixed order and soft all orders, one should use the fixed-order result up to (and including) O(α k s ), and the resummed result for the remaining terms from α k+1 s onwards. Matching is achieved by simply adding together the two calculations and subtracting the expansion of the resummation to O(α k s ), in order to avoid double counting. We then define res (N ) the coefficients of the expansion of C res (N, α s ) in powers of α s . Therefore, the matched cross section is written as is the inverse Mellin transform of the gluon-gluon luminosity times ∆ k C res , Eq. (4.1). By construction, is done through the public code TROLL [48], formerly ResHiggs, which changed name after the inclusion of Drell-Yan and DIS resummation [49]. Some details on the validation of the code are given in App. A.2. As the name implies (the meaning of TROLL is TROLL Resums Only Large-x Logarithms), the code does not compute the fixed order, but only the resummed contribution, Eq. (4.3), so the fixed order has to be supplied by an external code. In this work, we use the code ggHiggs [33], where we have implemented the new N 3 LO result [9][10][11][12]. We stress that any other fixed-order code can be used (e.g., future versions of ihixs [50]), provided the same setting is used in both codes. A new version of TROLL, v3.1, which interfaces directly with ggHiggs v3.1 is publicly available from the webpage [48].
For simplicity, and for disentangling effects coming from different sources, we work in the clean environment of the (rescaled) large-m t effective theory (rEFT), using the top mass m t = 172.5 GeV in the pole scheme. We take the Higgs mass to be m H = 125 GeV, and use the PDF4LHC15_nnlo_100 PDF set [34][35][36][37][38]. As far as α s is concerned, we take α s (m 2 Z ) = 0.118 from the PDF set and evolve it at four loops to µ R . We focus on LHC at √ s = 13 TeV first, and consider different collider energies later.
To show the stability of the resummed result, we consider four options for the central common factorization and renormalization scale µ 0 : We then vary the scales µ R and µ F about µ 0 by a factor of 2 up and down, keeping the ratio µ R /µ F never larger than 2 or smaller than 1/2: we call this a canonical 7-point scale variation. This results  The resummed results at N 3 LO+N 3 LL, obtained as the sum of the last column of Tab. 1 and the resummation contribution ∆ 3 σ N 3 LL computed with TROLL, are shown in Tab. 2.
It is well known [43] that threshold resummation introduces a dependence on the factorization scale which can be larger than what is obtained in fixed-order perturbation theory. This is essentially due to the fact that threshold resummation predicts only the (dominant) gg channel, while factorization scale dependence is compensated among different channels, as DGLAP evolution mixes quarks and gluons. Moreover, at fixed order the factorization scale dependence for Higgs production is very mild (and much milder than renormalization scale dependence, see Tab. 1), so the factorization scale dependence of the resummed result is visibly larger. The quark channels, of which qg gives the most important contribution, give rise to logarithmic terms that are suppressed by 1/N , in the large N limit. Including a prediction of this channel to all orders should compensate most of the factorization scale dependence. The resummation of the leading logarithms of this class of NS contributions has been performed in Ref. [51]. However, these contributions are not yet implemented in the current version of TROLL.
We now turn to our proposal for the perturbative uncertainty of our resummed results. We consider ψ-soft with AP2 and with the natural choice for the constants, Eq. (3.2), as our best option for threshold resummation. However, the other variants of ψ-soft have the same formal accuracy and allow us to estimate the uncertainty from 1/N terms and subleading logarithmic terms. We therefore suggest to consider, for each of the central scales µ 0 in Eq. (4.4), the envelope of the canonical 7-point scale variations and the 6 variants of ψ-soft resummation (we exclude N -soft from the computation). This corresponds to a total of 7 · 6 = 42 cross section points, 7 from which one takes the highest and the lowest cross sections as the maximum and minimum of the uncertainty band. As an example, we highlighted in Tab. 2 those 42 cross sections entering in the error band computation for the central scale µ 0 = m H /2. We conventionally take our default best option (shown in red in Tab. 2) as the central prediction.  Table 3. Let us first comment the fixed-order results. Ignoring the LO which contains too few information for being predictive, we can investigate the convergence pattern of the fixed-order perturbative expansion when going from NLO to NNLO and to N 3 LO, relative to the scale uncertainty. For "large" central scales, µ 0 = m H and µ 0 = 2m H , NNLO is a large correction and its central value is not covered by the NLO uncertainty band. The N 3 LO is a smaller correction, a sign that the series is converging (at least asymptotically), but for µ 0 = 2m H its central value is not covered by the NNLO uncertainty band. For µ 0 = m H /2, the convergence pattern is improved, now with the central NNLO contained in the NLO band, and the central N 3 LO contained in the NNLO band. However, for instance, the central N 3 LO and its band are not contained in the NLO band (they do not even overlap). At µ 0 = m H /4 the convergence pattern seems further improved, however the N 3 LO error is very asymmetric and large (same size of the NNLO error). Additionally, the N 3 LO results at the four central scales shown in Table 3 are barely compatible (if one had chosen µ 0 = 4m H the result would not be compatible with the one at µ 0 = m H /2). This analysis shows that the estimate of the uncertainty from missing higher orders using canonical 7-point scale variation is not reliable at fixed order. This is perhaps not surprising, as scale variation provides a very crude estimate of the uncertainty from missing higher orders, since it is based on arbitrary variation of a (not necessarily significant) subset of known coefficients.
On the other hand, resummation allows for a different way of estimating the effect of missing higher orders, which is not purely based on scale variation. We observe that, for each choice of the central scale µ 0 , the uncertainty of the resummed results from NLO+NLL onwards covers the central value and at least a portion of the band of the next (logarithmic) order. In fact, with the exception of the choice µ 0 = m H /4 (the pathological behaviour of which seems to be driven by the N 3 LO contribution), the NNLO+NNLL band is fully contained in the NLO+NLL band, and the N 3 LO+N 3 LL band is fully contained in the NNLO+NNLL band. We also note a systematic reduction of the scale uncertainty  when going from one logarithmic order to the next. We also observe that the resummed results at each order are all compatible among the different choices of the central scale µ 0 , thereby showing little sensitivity on µ 0 . It is true that at extreme choices of µ 0 the error bands become very asymmetric and lead to higher values of the cross section at large µ 0 and to lower values of the cross section at small µ 0 ; nevertheless, a region of overlap always exists.
We note that our observations on the behaviour of the resummed results would still hold if one considers a less conservative option, namely our default ψ-soft resummation with AP2 and the natural choice for the constants, Eq. (3.2), corresponding to the red dots and the thick red bands in the plots. In fact, with this option the errors would look more natural, especially at large µ 0 where the AP1 variation (thinner band) increases the size of the error band dramatically. It is interesting to observe that the different options for the position of the constants, while giving a large spread at NLO+NLL, is of little importance at higher orders, especially at N 3 LO+N 3 LL. This is particularly true for smaller µ 0 , while at larger µ 0 the version with all constants in g 0 , Eq. (3.10), gives a slightly smaller cross section.
We finally observe that, in many respects, the choice µ 0 = m H /2 seems optimal, in full agreement with previous analyses, e.g. [12]. The convergence of the fixed-order is already quite good, and the convergence of the resummed result is very good. The error band at N 3 LO+N 3 LL is smaller than for other central scales, but compatible with the results computed at different values of µ 0 . On top of these a posteriori observations, one could determine a priori an optimal choice for the scale by requiring that the partonic coefficient functions do not contain large logarithms, such that possible logarithmic enhancements are minimized. Factorization and renormalization scales typically appear together with threshold logarithms, in the form (in Mellin space) ln(µ 2 N 2 /m 2 H ). From a saddle point analysis [52], we know that the Mellin space cross section is dominated by a single value of N = N saddle , leading to an optimal choice for the scales µ 0 m H /N saddle . In the present case N saddle 2, so we find that the scale that minimizes the size of the logarithms is close to µ 0 = m H /2. Similar conclusions have been obtained from z-space arguments [53] and within an effective theory framework [54,55]. As the process gets closer to threshold, N saddle grows and the optimal scale gets correspondingly smaller: we have indeed verified that closer to threshold (larger m H or smaller collider energy) the perturbative convergence is much improved at smaller scales.
Given that the way of estimating the uncertainty is very conservative, and successful at previous orders, the uncertainty on the N 3 LO+N 3 LL at µ 0 = m H /2 seems reasonably trustful. To be even more conservative, one can symmetrize the error, so rather than 48.5 +1.5 −1.9 pb we can quote 48.5 ± 1.9 pb as most reliable prediction for the inclusive Higgs cross section (in the rEFT setup). Note also that at µ 0 = m H /2 the effect of adding the resummation to the N 3 LO on the central value is rather small, +0.4 pb, corresponding to +0.8%, which is not covered by the fixed-order uncertainty. However, we find the asymmetric error on the N 3 LO hard to trust. Indeed, this is due to the vicinity of a stationary point in the scale dependence, which in fact shows that such a scale error is not reliable. Thus, at the very least, one should symmetrize it, leading to 48.1 ± 1.8 pb, which is then compatible with what we obtain from the resummation procedure.
To emphasize the robustness of our method for the estimate of the perturbative uncertainty, we repeat our analysis for different collider energies: √ s = 2, 8, 14 and 100 TeV. We show in Fig. 3 the analogous of the µ 0 = m H /2 plot of Fig. 2 for the aforementioned collider energies. We observe the same pattern found at 13 TeV. Note that for a FCC-like energy of 100 TeV the smaller values of z accessible at that energy will make the correct inclusion of small-z effects (at fixed-order or resummed levels) very important. We stress in particular that the rEFT is not able to predict the small-z behaviour correctly, and one should revert to the full theory for an appropriate description.

Other ways to estimate the uncertainty from missing higher orders
In the previous section we have proposed a robust way of estimating the uncertainty from missing higher orders using scale variation and variation of subleading terms in the all-order resummation. In this section we want to explore alternative strategies to estimate the uncertainty from missing higher orders which do not rely on arbitrary variations of the perturbative ingredients. In particular, we will consider in Sect. 5.1 the Cacciari-Houdeau method [19][20][21] for estimating the theory uncertainty using a Bayesian approach to infer the uncertainty from the progression of the perturbative expansion.
Then, in Sect. 5.2, we will follow the idea of Ref. [22] to apply convergence acceleration algorithms to the perturbative series (either the fixed-order or the resummed one) to estimate the truncation error. We shall then compare the findings of these methods with our results of Sect. 4.

The Cacciari-Houdeau method
In Ref. [19] Cacciari and Houdeau proposed a statistical model for the interpretation of theory errors, from which one can compute the uncertainty on the truncated perturbative series for a given degree of belief (DoB) given the first terms in the expansion. Among the assumptions of the model, there is the fact that all coefficients c i of the expansion are bounded by a common valuec. To account for potential power growth of the coefficients, the expansion parameter α s is rescaled, giving the expansion In Eq. (5.1) we have factored out the LO cross section σ 0 such that the sum starts from k = 0, with c 0 = 1. The coefficients c k depend on the rescaling factor λ, which should be determined such that the boundc exists. We come back on the determination of λ later in this section. This is the original method, denoted CH. It has been noted [19] that the assumption that all c k are bounded is surely broken by the presence of known factorial growths in the coefficients, such as those due to renormalons. However, it was pointed out that the growth will set in at a high order, so for practical applications at low orders ignoring it is harmless. However, in Ref. [21] this factorial growth is instead included in the definition of the series expansion, where the new coefficients b k do not contain the factorial growth anymore. This is the modified Cacciari-Houdeau method, denoted CH. We included an explicit offset k 0 in the factorial: this is useful since different observables will have the factorial growth starting at different orders. Note that in the original works [19,21] there is a distinction between observables starting at different orders α l s . In our case, having factored out the LO, the complications arising from this distinction go away. Factoring out a power of α s does not change the error estimate in the original CH method, so our Eq. (5.1) gives identical results as those obtained with the original formulation of Ref. [19]. On the other hand, due to the presence of the factorial, factoring out a power of α s in the modified CH method would make a difference if the factorial is left unchanged. The inclusion of the offset k 0 in the factorial in Eq. (5.2) accounts for this difference and allows to exactly reproduce the results of Ref. [21] for the specific choice k 0 = l − 1.
Given the set of N k LO coefficients c 0 , . . . , c k (or b 0 , . . . , b k in the modified version) one can construct the credibility interval [−d k ] in which the remainder of the perturbative series is expected to lie with DoB p. We have [19,21] CH: d Note that these simple analytic expressions assume that the remainder of the perturbative series is dominated by the next higher order; this approximation is good if the expansion parameter α s /λ is small enough, but breaks down for small values of λ. In such cases, one should use the full result [19], which is not expressible in closed analytic form, or interpret the resulting uncertainty as arising just from the next order. It now remains to determine the values of k 0 and λ. Regarding k 0 , in Ref. [21] it is argued that for a series starting at order α l s the proper value is k 0 = l − 1. This follows from the observation that for weak processes starting at order α 0 s the renormalon factorial growth behaves as α k s (k − 1)!. Therefore, for Higgs production which starts at order α 2 s the suggested value is k 0 = 1. However, we point out that what matters here is where the perturbative corrections to gluon propagators start rather than the power of α s at LO. Indeed, for processes like Drell-Yan, the gluon appears first at NLO, and the correction to the gluon propagator only at NNLO, so indeed the factorial growth is delayed by one order and k 0 = −1. However, for Higgs production, there are gluons already at LO, so the first correction appears at NLO, from which one can conclude that k 0 = 0. We adopt here the latter choice (k 0 = 0) rather than the default option of Ref. [21] (k 0 = 1).
We now turn to the determination of λ. This parameter is in principle a free parameter, but it must be such that there exists a bound for the coefficients of the expansion. In Ref. [21] λ is determined from a survey over several processes, giving λ = 0.6. However, different perturbative expansions can behave in very different ways, and in particular the perturbative expansion for Higgs production in gluon fusion is badly behaved, so a separate treatment is to be preferred. Ideally, the value of λ should be determined according to the asymptotic behaviour of the perturbative coefficients. Since this is unfortunately not known, we follow the proposal of Ref. [20], where the value of λ is fitted such that the perturbative coefficients (in absolute values) are all of the same size. In fact, as we also confirmed, it is convenient to exclude the first coefficient from the fit, on the ground that the LO result is not in line with the next orders (it is much smaller), and the fact that this fit aims at guessing the asymptotic behaviour of the coefficients. In the results that follow, we will then use for each method (CH and CH) and for each central scale the value of λ obtained by such fit. As we have observed, however, this determination of λ risks being ad hoc and therefore may result in a somewhat biased error.
In Fig. 4 we show the four results at LO, NLO, NNLO and N 3 LO for the four scales µ F = µ R = m H /4, m H /2, m H , 2m H , each with the two versions (CH and CH) of the Cacciari-Houdeau uncertainty. We observe that the CH uncertainty is larger than the CH one at LO and NLO, but is smaller at NNLO and N 3 LO: this effect originates from the factorial contribution, which changes the relative weight of the individual orders in the determination of the uncertainty. In this respect, the CH uncertainty at N 3 LO is more conservative than the CH one. We also note that the 68% DoB uncertainty (thicker band) is able to cover the next order only at NNLO, while for lower orders only the 95% DoB uncertainty (thinner band) works (except at LO for CH). We also see that for small scales µ F = µ R = m H /4, m H /2 the uncertainty shrinks considerably as the perturbative order increases, an indication that the series is converging. For larger scales, µ F = µ R = m H , 2m H , the observed pattern is much worse and, as a consequence, the uncertainty band of the N 3 LO is still large.
In Tab. 4 we report the value of N 3 LO cross section together with its uncertainty as obtained with CH and CH, for four different choices of the central scale. We indicate both 68% and 95% DoBs (the latter in brackets). If we focus on the default choice m H /2, we note that the estimate of the uncertainty due to missing higher orders as obtained combining scale and resummation uncertainties (see N 3 LL+N 3 LO in Tab. 3) is in agreement with the CH uncertainty estimate at 95% DoB (while CH provides with a smaller uncertainty). In contrast, as previously noted, scale variation on its own (e.g. N 3 LO in Tab. 3) provides us with highly asymmetric error, which appears to underestimate the   upper portion of the uncertainty band.

Convergence acceleration algorithms
In this last section we want to explore another method to gain information on the uncertainty from missing higher orders by estimating the sum of the perturbative series. Our strategy is based on convergence acceleration algorithms: given a sequence which converges to some limit there exist several algorithms which transform the sequence into new ones with possibly faster convergence. Non-linear sequence transformations usually provide faster convergence than linear transformations. We follow an idea by David and Passarino [22] and we apply some of these methods to the sequence of partial sums of the perturbative expansion of the cross section.
Following Ref. [56], we consider the following sequence transformation: where s n = n i=0 a i is the n-th partial sum, and a i the series coefficients. For k > 0, it provides a non-trivial transformation of the original sequence. The non-linear transformation Eq. (5.6) depends on the function q m , and for particular forms of it reduces to transformations widely studied in the literature [56], e.g. q m = β > 0 gives a Levin transformation L k (β, s n , ω n ) transformation, and q m = β + (m − 1)/α gives a transformation C (n) k (α, β, s n , ω n ) which interpolates between L and S depending on the value of α. 8 These transformations further depend on the function ω n , which is related to the remainder estimate of the truncated sequence. In our case we make no assumptions about the asymptotic behaviour of the perturbative series of the Higgs cross section, mainly because the knowledge of only the first 4 terms is not sufficient to guarantee that the asymptotic behaviour has set in. Therefore, we consider various forms of the remainder estimate, leading to a number of variants of the aforementioned methods [56]: u-type variant given by ω n = (β + n)a n (with an overall minus sign for the M transformation), t-type variant given by ω n = a n , d-type variant given by ω n = a n+1 , and v-type variant given by ω n = a n a n+1 /(a n − a n+1 ).
For u-and t-type variants, the transformation G (n) k requires the knowledge of k + n + 1 terms in the sequence, while for d-and v-type variants k + n + 2 terms are needed. At such low orders, several transformations turn out to be redundant; after removing equivalent transformations, the non-trivial remaining ones are u L 2 ), known as Weniger δ transformation. We recall that each of these variants further depends on the variable β, and the C transformations also depend on α.
In this study, we use all the non-equivalent transformations listed above and scan for various values of β and α. Our idea is that, given the unknown higher-order behaviour of the perturbative series, it is impossible to choose a particular algorithm over the other ones or tune the parameters of the transformations. Indeed, non-linear sequence transformations are in general not guaranteed to work in all cases, and privileging a specific algorithm would require the knowledge of the asymptotic behaviour of the series, which we have not. However, having at hand several different acceleration algorithms, we can judge a posteriori how stable the estimate of the sum of the series is by comparing among the different predictions: when most results cluster around the same value, we can expect it to be a reliable estimate of the real sum.
As previously mentioned, there are contributions to the perturbative expansion (at parton level) which grow factorially (due e.g. to renormalons). Even though some of the aforementioned sequence transformations proved to be effective also in case of factorially divergent series, we also consider here a more standard method based on Borel summation, where the sum of the divergent perturbative 25   series is defined as When a finite number of coefficients is known, the Borel sum Eq. (5.7) becomes an identity. To all orders, the left-hand side series can diverge due to renormalons, while the sum on the right-hand side is typically convergent, due to the factorial that cures the growth. The integral can then be computed: if the integral is finite, 9 then the right-hand side of Eq. (5.7) can be used to define the (Borel) sum of the divergent series. In our case, we can benefit from the Borel summation by applying a convergence acceleration algorithm to compute the all-order sum on the right-hand side, thereby curing possible divergences of the original series. Hence, for each of the methods listed above, we can produce a Borel variant, doubling the available options.
Having implemented all these variants, we apply them to the fixed-order and resummed 10 ex- 9 In general, the integrand on the right-hand side of Eq. (5.7) can have poles on the integration range. This problem can be solved by modifying the integration contour and avoiding the singularity from above or below, leading to a well-known renormalon ambiguity. In our case, the ambiguity is always very small, so we ignore it. 10 For the resummation we use ψ-soft AP2, with default option for the constants. Resummed expansion 48.9 ± 0.5 48.9 ± 0.6 50.2 ± 1.0 52.6 ± 1.6 Table 5. Mean and standard deviation of the estimates of the all-order sum (in pb) of the fixed-order (first row) and resummed (second row) expansions, based on the set of convergence acceleration algorithms described in the text, for mH = 125 GeV at LHC with √ s = 13 TeV.
pansions for each of the previously considered central scales µ 0 = m H /4, m H /2, m H , 2m H . For each method, we compute the results for several values of β; applying these methods to a number of known series we identified β = 0.01, 1, 2, 5 as a sensible set. For the C transformations we also scan over α = 0.5, 2, 10, 100: this is motivated by the fact that a value of α between 1 and +∞ interpolates between the S and L transformations, and we also considered a value outside this range. Note that u L do not depend on β, so these are counted only once. We collect the results for the estimated sum of the fixed-order and resummed series in the form of histograms in Fig. 5.
We immediately observe that the distribution of results is narrower for resummed results than for fixed-order ones. This can be expected since the resummed expansion is much better behaved than the fixed-order one, as commented extensively in Sect. 4. Additionally, we observe that the distributions for µ 0 = m H /4 and µ 0 = m H /2 are narrower than those at higher scales, a fact which is especially true for the fixed-order expansion. This confirms the faster convergence of the perturbative series at µ 0 = m H /2 (and also µ 0 = m H /4) observed in previous sections and discussed in the literature. The fact that for µ 0 = m H and µ 0 = 2m H the spread of the results is rather wide shows that the algorithms considered here are sufficiently different among each other. This is important because it validates the set of algorithms we have chosen, which in turn validates the rather precise results obtained at lower scales.
To make the discussion more quantitative, we report in Tab. 5 the mean and standard deviation for each histogram of Fig. 5. We warn the Reader that the statistical interpretation of these results is not solid: in particular, the standard deviations are likely to depend on the set of convergence acceleration algorithms considered.
All the numbers in Tab. 5 come from estimates of the all-order sum of the series, which should be then the same for all scales and for both the fixed-order and the resummed expansions. They are indeed all compatible within the quoted errors, except the resummed result at µ 0 = 2m H which is higher than most of the other results: this is just a consequence of the limited statistical meaning of the error estimates, which does not take into account the shape of the distribution of the results, which is rather asymmetric in this case. The smaller standard deviation on the resummed results shows once again that the resummed series converges faster, as well as the smaller standard deviation on the results at lower scales indicates that using µ 0 = m H /2 or µ 0 = m H /4 leads to a faster convergence, in agreement with the findings of our previous sections. It is interesting to observe that at both scales µ 0 = m H /2 or µ 0 = m H /4 and for both fixed-order and resummed expansions the estimate of the sum is basically the same (48.7 from fixed-order and 48.9 from resummation), both with a small standard deviation (of the order of 1 pb). Keeping in mind the limitation of this analysis, we are tempted to consider this result as a good candidate for the all-order sum of the series. Interestingly, this result is perfectly compatible with (and very close to) our best N 3 LO+N 3 LL result at µ 0 = m H /2 well within its ±1.9 pb uncertainty. This provides another valuable validation of our proposal for estimating missing higher-order uncertainty from resummation. On the other hand, it is not compatible with the N 3 LO result within its asymmetric scale-variation band, while it is considering a CH error already at 68% DoB (see Tab. 4).

Conclusions
We have presented threshold-resummed results for the inclusive Higgs cross section in gluon fusion at N 3 LL, matched to an implementation of the recent N 3 LO result [9][10][11][12]. We have considered several variants of the resummation as a portal to carefully estimate subleading effects at higher orders. We have proposed a conservative estimate of the uncertainty from missing higher orders based on the envelope of the resummed predictions obtained using the various resummation variants, as well as canonical scale-variation. We have demonstrated that resummed results with this conservative error manifest a good perturbative convergence, as opposed to the fixed-order expansion, the convergence of which is very poor relative to the uncertainty coming from a canonical 7-point scale variation.
Despite the conservativeness of our method, we find that the Higgs cross section at 13 TeV, for the central scale µ R = µ F = m H /2, has a small (yet reliable) uncertainty of ±1.9 pb, which corresponds to ±4%. We stress that all choices of central scales considered in this work (m H /4, m H /2, m H and 2m H ) yield results which are compatible within the quoted uncertainty. The shift in the central value and the uncertainty, though computed within the framework of the (rescaled) large-m t effective theory, are likely to remain unchanged after inclusion of quark mass effects and Electro-Weak corrections. For the most reliable predictions the inclusion of quark mass effects is important, and can be performed straightforwardly at resummed level [13,16] with TROLL. Moreover, a fully consistent resummed result would require the use threshold-improved parton distribution functions, which have recently become available [49].
We have compared our proposal with different methods for estimating the uncertainty from missing higher orders. Our findings are summarized in the following from accelerated resummed series First, we have considered the Cacciari-Houdeau Bayesian approach, which employs the known perturbative orders to construct a probability distribution for the subsequent unknown order. In its modified incarnation (CH), the method gives an uncertainty of ±2 pb at 95% degree of belief, fully compatible with the estimate obtained from resummation, and similar to the fixed-order scale variation uncertainty if the latter is symmetrized. Second, we have considered several algorithms to accelerate the convergence of the perturbative series, based on non-linear sequence transformations. By performing a survey of different algorithms, we found that both the fixed-order and resummed series exhibit good convergence properties at m H /2 (and also at m H /4). Noticeably, the mean of each distribution is very close to the N 3 LO+N 3 LL prediction. In conclusion, these tests provide a solid support to our method, and let us conclude to a high degree of belief that the all-order Higgs cross section in the rEFT lies within the quoted uncertainty of our N 3 LO+N 3 LL result.
which has been computed analytically, is indeed correct. To do so, we compared the analytic fixedorder expansion of the resummation to a numerical expansion of the resummation itself. We report here the results of this comparison. Recalling Eq. (4.1) we have that the resummed expression can be expanded in powers of α s as For matching to N 3 LO the expansion coefficients C (1) res (N ), C res (N ) and C (3) res (N ) are needed, and they have been computed analytically and coded in TROLL. Each coefficient can be also extracted numerically according to the formula ): therefore, the → 0 limit suppresses higher-order corrections and isolates the O(α k s ) term. In Fig. 6 we show the analytic results for the cross section contributions corresponding to α k s C (k) res (N ) for k = 1, 2, 3 together with the numerical expression Eq. (A.3) as a function of . We see that at small the numerical expansion reproduces the analytic result for several non-trivial combinations of the scales µ R and µ F . This represents a strong check of the implementation of the matching to fixed-order for any µ F and µ R . In particular, it ensures that the contribution from the resummation is always one order higher than the fixed order we are matching to.
In order to cross check the scale dependence of the resummed result, we now verify the explicit scale dependence of the expansion of the resummation to fixed order, which, as we have just verified, is consistent with the all-order expression. To do so, we compare the results obtained from the internal implementation of the scale dependence with an external implementation. To be precise, we compare the contributions to the cross section (where for simplicity we removed overall factors from the definition of σ (k) ) to the analogous computed with µ F = µ R = m H plus scale dependent terms In this way, in the second expression Eq. (A.5) all the explicit logs of µ R and µ F in the partonic coefficients are set to zero, and the scale dependent terms are provided with the additional contribution ∆σ ) anomalous dimensions in the soft limit (explicit expressions are given in Ref. [57]). We have verified that the two expressions Eq. (A.4) and Eq. (A.5) give identical results for any combination of µ R and µ F at orders k = 1, 2, 3. This represents an independent cross-check of the scale dependence of the resummed expression.