High precision determination of the gluon fusion Higgs boson cross-section at the LHC

We present the most precise value for the Higgs boson cross-section in the gluon-fusion production mode at the LHC. Our result is based on a perturbative expansion through N$^3$LO in QCD, in an effective theory where the top-quark is assumed to be infinitely heavy, while all other Standard Model quarks are massless. We combine this result with QCD corrections to the cross-section where all finite quark-mass effects are included exactly through NLO. In addition, electroweak corrections and the first corrections in the inverse mass of the top-quark are incorporated at three loops. We also investigate the effects of threshold resummation, both in the traditional QCD framework and following a SCET approach, which resums a class of $\pi^2$ contributions to all orders. We assess the uncertainty of the cross-section from missing higher-order corrections due to both perturbative QCD effects beyond N$^3$LO and unknown mixed QCD-electroweak effects. In addition, we determine the sensitivity of the cross-section to the choice of parton distribution function (PDF) sets and to the parametric uncertainty in the strong coupling constant and quark masses. For a Higgs mass of $m_H = 125~{\rm GeV}$ and an LHC center-of-mass energy of $13~{\rm TeV}$, our best prediction for the gluon fusion cross-section is \[ \sigma = 48.58\,{\rm pb} {}^{+2.22\, {\rm pb}\, (+4.56\%)}_{-3.27\, {\rm pb}\, (-6.72\%)} \mbox{ (theory)} \pm 1.56 \,{\rm pb}\, (3.20\%) \mbox{ (PDF+$\alpha_s$)} \]

Abstract: We present the most precise value for the Higgs boson cross-section in the gluon-fusion production mode at the LHC. Our result is based on a perturbative expansion through N 3 LO in QCD, in an effective theory where the top-quark is assumed to be infinitely heavy, while all other Standard Model quarks are massless. We combine this result with QCD corrections to the cross-section where all finite quark-mass effects are included exactly through NLO. In addition, electroweak corrections and the first corrections in the inverse mass of the top-quark are incorporated at three loops. We also investigate the effects of threshold resummation, both in the traditional QCD framework and following a SCET approach, which resums a class of π 2 contributions to all orders. We assess the uncertainty of the cross-section from missing higher-order corrections due to both perturbative QCD effects beyond N 3 LO and unknown mixed QCD-electroweak effects. In addition, we determine the sensitivity of the cross-section to the choice of parton distribution function (PDF) sets and to the parametric uncertainty in the strong coupling constant and quark masses. For a Higgs mass of m H = 125 GeV and an LHC center-of-mass energy of 13 TeV, our best prediction for the gluon fusion cross-section is σ = 48.58 pb

Introduction
With the discovery of the Higgs boson [1,2], the Large Hadron Collider (LHC) has achieved a major landmark in science. Indeed, we now have conclusive evidence that space is filled with the Higgs field and that the mass of elementary particles is not an ad-hoc concept, but an elaborate outcome of the mechanism of spontaneous symmetry breaking. Moreover, with the Higgs boson the Standard Model is a mathematically self-consistent theory, and it can be used to formulate physically credible predictions at extremely high energies, many orders of magnitude higher than what we can probe with man-made experiments. This is a great triumph of theoretical physics.
Besides the success of the Standard Model as a theory of electroweak interactions, it is a phenomenologically incomplete theory, and it needs to be extended in order to obtain a satisfying description of all known physics, including cosmology. It is unclear at what energy the Standard Model will stop being a good theory and it will require the introduction of new laws of physics. If open questions, such as for example the origin of dark matter, are related to the question of the origin of mass of elementary particles, then it is likely that Higgs phenomena will differ quantitatively from Standard Model expectations. We should therefore view the Higgs boson discovery as the foundation of a long-term precision physics program measuring the properties of the Higgs boson. This program may yield direct or indirect evidence of physics beyond the Standard Model, and it requires the measurement of the mass, spin/parity, width, branching ratios and production rates of the Higgs boson. All of the above are predicted or constrained in the Standard Model and its viable extensions. The success of the program will rely crucially on the combination of highly precise experimental data with equally accurate theoretical predictions.
The purpose of this article is to supply a key ingredient to upcoming high-precision studies of the Higgs boson by providing the most accurate determination of the Higgs production cross-section in gluon fusion. Higgs production in gluon fusion is mediated mainly through a top-quark loop [3]. Despite the absence of a tree-level contribution (which makes gluon fusion a pure quantum process), it is the dominant production mode of the Higgs boson due to the large gluon luminosity and the size of the top-quark Yukawa coupling. Reliable predictions of this process require the inclusion of higher-order corrections, both from the QCD and the electroweak sectors of the Standard Model. The phenomenological importance of these corrections can be seen from the large size of the next-to-leading order (NLO) QCD corrections [4,5,6,7,8,9,10,11,12], which almost double the original leading order (LO) prediction [3]. The large size of the NLO corrections indicate potentially significant contributions from even higher perturbative orders, thus resulting in a substantial theoretical uncertainty on the gluon-fusion cross-section. This uncertainty is difficult to quantify conventionally by varying nuisance parameters in the theoretical prediction such as renormalization and factorization scales.
The past twenty years have seen substantial theoretical advances in the perturbative description of Higgs production in gluon fusion, using a multitude of techniques and aiming for various directions of improvement. The theoretical description of gluon fusion is rendered particularly challenging by the fact that the Born process is already a one-loop process involving two mass scales (the masses of the top quark and the Higgs boson), such that higher-order corrections will involve multi-scale multi-loop amplitudes. This challenge can be overcome by integrating out the top quark at the level of the Standard Model Lagrangian. This procedure results in an effective field theory (EFT) [13,14,15,16] containing a tree-level coupling of the Higgs boson to gluons. This EFT can be matched onto the full Standard Model in a systematic manner, resulting in corrections from higher orders in perturbation theory to the Wilson coefficients [17,18,19] and from subleading terms in the large mass expansion [20,21]. In the EFT framework the coefficient function for inclusive Higgs production in gluon fusion depends only on the ratio of Higgs boson mass m H to the partonic center-of-mass energy √ s, usually expressed through the variable z = m 2 H /s. By comparing the full NLO QCD expression for the gluon fusion cross-section with the EFT result, one observes that a very good approximation of the full NLO prediction can be obtained by re-weighting the EFT prediction with the ratio of LO predictions in the full and effective theories. This re-weighting is commonly applied to all predictions obtained within the EFT framework.
NNLO corrections in the EFT [22,23,24] turn out to be substantial, albeit smaller than at NLO. This slow convergence pattern entails the risk that the estimation of the uncertainty at NNLO, obtained by varying the renormalization and factorization scales, may be misleading. To settle the size of the QCD corrections, additional information about the behavior of the perturbative expansion beyond NNLO is necessary.
A part of these corrections arises from the region of phase-space where the Higgs boson is produced at or near to its kinematical threshold, z → 1. In this region, contributions from soft and collinear gluon emission can be resummed to all orders in the coupling constant, either using Mellin space methods [25,26,27,28,29] or within soft-collinear effective field theory (SCET) [30,31,32,33,34]. Always working in the EFT framework, either method enables the resummation of logarithmically-enhanced threshold corrections to the gluon fusion process to high logarithmic order [35,36,37]. Both methods agree to a given formal logarithmic accuracy, but their results can in principle differ [38,39] by non-logarithmic terms.
The combination of predictions at fixed order with threshold resummation (together with electroweak corrections [40,41,42,43], bottom and charm quark contributions through NLO and subleading mass corrections at NNLO [20,21,44,45]) provided the default predictions for the interpretation of Higgs production data in the LHC Run 1 [46,47,48]. Besides uncertainties from the parametrization of parton distributions and the values of the strong coupling and quark masses, these predictions were limited in accuracy by missing terms at N 3 LO in the fixed-order expansion. By expanding resummed predictions in powers of the coupling constant, logarithmic terms in perturbative orders beyond NNLO can be extracted. These were used (often combined with the knowledge of the high-energy z → 0 behaviour of the coefficient function [49]) to obtain estimates of the fixed-order gluon fusion cross-section at N 3 LO and beyond [50,51,52]. Although individual results for these estimates typically quoted a very small residual uncertainty, the scatter among different estimates was quite substantial, thereby putting serious doubts on the reliability of any such estimation procedure. These ambiguities were resolved by the recent calculation of the full N 3 LO QCD corrections to the gluon fusion process [53], which are the key ingredient to the work presented here.
The gluon-fusion cross-section in N 3 LO QCD in the EFT approach receives (besides a network of lower order renormalization and mass-factorization terms [54,55,56,57,58]) contributions from four types of processes, ranging at fixed sum of loops and external legs from three-loop virtual corrections to the ggH-vertex to triple real radiation corrections from processes like gg → Hggg, and denoted respectively as VVV, (RV) 2 , RVV, RRV, RRR. While the three-loop virtual corrections were already known for quite some time [59,60], new technical advances were needed in order to evaluate all the contribution from the different real-radiation subprocesses, either in closed from or as a high-order expansion around the threshold limit that is sufficient to precisely account for the full z-dependence. Based on the two-loop matrix elements for Higgs-plus-jet production [61], closed expressions for the RVV contributions could be obtained by direct integration [62,63]. In the same way, it was possible to derive closed expressions for the (RV) 2 contribution [64,65]. The major challenge in the RRV and RRR processes are the very intricate phase space integrals for double real radiation at one loop and triple real radiation at tree level. These phase-space integrals can be related to specific cuts of loop integrals using reverse-unitarity [22,66,67,68], allowing the application of modern integral reduction techniques [69,70,71] that express all relevant phase-space integrals by a limited set of master integrals, which are functions of z. With the same integral reduction techniques, differential equations [72,73,74] can be derived for the master integrals. By solving these differential equations (either in closed form [75], or as an expansion in z) for appropriate boundary conditions, the direct integration of the master integrals can be circumvented. The use of these techniques enabled the computation of the RRV and RRR [76] contributions at threshold. Combining them with the two-loop correction to the soft-gluon current [77,78] enabled first breakthroughs with the N 3 LO threshold cross-section [79,80,81] and the first beyond-threshold term [82]. More recently, the systematic expansion of the RRV [83] and RRR contributions to very high orders in z enabled the calculation of the full N 3 LO gluon fusion cross-section [53].
In this publication, we combine the N 3 LO cross-section in the EFT with the previously available state-of-the-art predictions for other types of corrections (electroweak, mass effects, resummation) to obtain a highly precise theoretical description of the inclusive Higgs production cross-section in gluon fusion. We also assess remaining theoretical uncertainties on the cross-section. Our results will form a cornerstone to precision studies of the Higgs boson in the upcoming high-energy and high-luminosity data taking periods of the CERN LHC.
Our paper is structured as follows: in Section 2 we present our setup and summarize the different contributions that we include into our prediction. In Section 3 we study the phenomenological impact of QCD corrections through N 3 LO in the large m t -limit. We investigate the missing higher-order effects and threshold resummation in the EFT in Section 4. Effects due to quark masses and electroweak corrections are studied in Sections 5 and 6, and we assess the uncertainty on the cross-section due to PDFs and the strong coupling constant in Section 7. In Section 8 we combine all effects and present our recommendation for the most precise theoretical prediction of the inclusive Higgs cross-section. In Section 9 we draw our conclusions. We include appendices where we present the coefficients appearing in the threshold expansion of the N 3 LO coefficient, as well as tables summarizing our results for a variety of different Higgs masses and collider energies.

Setup
The inclusive hadronic cross-section σ for Higgs production in gluon fusion can be calculated as the convolution integral where m H , s and S denote the Higgs mass and the squared partonic and hadronic centerof-mass energies. The convolution of two functions is defined as In the Standard Model (SM) the Higgs boson is predominantly produced through the annihilation of virtual top and bottom quarks, as well as W and Z bosons, produced in gluon fusion. All of these channels are greatly enhanced by QCD corrections, and also electroweak corrections are important. Hence, having good control over higher-order corrections in perturbation theory, both in the QCD and electroweak sectors, is of paramount importance to make precision predictions for Higgs production in the framework of the SM. We note that non-perturbative contributions to the inclusive Higgs boson production cross-section are suppressed by powers of (Λ/m H ), with Λ being the QCD scale. For the Drell-Yan process, the linear non-perturbative correction could be shown to vanish [84,85], such that the leading power correction term is quadratic and potentially relevant at low invariant masses. For the inclusive Higgs production cross-section, no thorough investigation of the power corrections has been performed up to now, but even the linear term would be a per-mille-level correction. The goal of this paper is to provide the most precise predictions for the inclusive hadronic Higgs production cross-section in the SM. We use state-of-the-art precision computations for electroweak and QCD corrections to inclusive Higgs production and combine them into the most precise theoretical prediction for the Higgs cross-sections available to date. The master formula that summarizes all the ingredients entering our prediction for the (partonic) cross-sections iŝ σ ij R LO σ ij,EF T + δ tσ N N LO ij,EF T + δσ ij,EW + δσ LO ij,ex;t,b,c + δσ N LO ij,ex;t,b,c .
(2.4) Equation (2.4) includes QCD corrections to the production cross-section in an effective theory where the top quark is infinitely heavy and has been integrated out. In this limit, the Higgs boson couples directly to the gluons via an effective operator of dimension five, where H is the Higgs boson field, G a µν is the gluon field strength tensor and L SM,5 denotes the SM Lagrangian with N f = 5 massless quark flavours. The Wilson coefficient C is obtained by matching the effective theory to the full SM in the limit where the top quark is infinitely heavy. In Appendix A we give its analytic expression through N 3 LO in the MS and OS schemes [17,18], in the five-flavour effective theory with the strong coupling constant decoupled.
QCD corrections to the production cross-sectionσ ij,EF T in the heavy-top limit have been computed at NLO [4,5,6] and NNLO [22,23,24]. Recently also the N 3 LO corrections have become available [53]. One of the main goals of this work is to combine the N 3 LO corrections in the large-m t limit with other effects that can provide corrections at a similar level of accuracy, in particular quark-mass effects and electroweak corrections. We also investigate the impact of the resummation of threshold logarithms, both within the frameworks of exponentiation of large logarithms in Mellin space and using soft-collinear effective theory (SCET).
While the production cross-section is known to high accuracy in the framework of the effective theory, reaching a similar level of accuracy when including quark-mass effects (also from bottom and charm quarks) is currently beyond our technical capabilities. Nonetheless, various quark-mass effects have been computed, which we consistently include into our prediction (2.4). First, it was already observed at LO and NLO that the validity of the effective theory can be greatly enhanced by rescaling the effective theory by the exact LO result. We therefore rescale the cross-sectionσ ij,EF T in the effective theory by the ratio where σ LO ex;t denotes the exact (hadronic) LO cross-section in the SM with a massive top quark and N f = 5 massless quarks. Moreover, at LO and NLO we know the exact result for the production cross-section in the SM, including all mass effects from top, bottom and charm quarks. We include these corrections into our prediction via the terms δσ (N )LO ij,ex;t,b,c in eq. (2.4), consistently matched to the contributions from the effective theory to avoid double counting. As a consequence, eq. (2.4) agrees with the exact SM cross-section (with massless u, d and s quarks) through NLO in QCD. Beyond NLO, we only know the value of the crosssection in the heavy-top effective theory. We can, however, include subleading corrections at NNLO in the effective theory as an expansion in the inverse top mass [20,21,44,45]. These effects are taken into account through the term δ tσ N N LO ij,EF T in eq. (2.4), rescaled by R LO .
Finally, we also include electroweak corrections to the gluon-fusion cross-section (normalised to the exact LO cross-section) through the term δσ ij,EW in eq. (2.4). Unlike QCD corrections, electroweak corrections have only been computed through NLO in the electromagnetic coupling constant α [40,41,42]. Moreover, mixed QCD-electroweak corrections, i.e., corrections proportional to α α 3 s , are known in an effective theory [43] valid in the limit where not only the top quark but also the electroweak bosons are much heavier than the Higgs boson. In this limit the interaction of the Higgs boson with the W and Z bosons is described via a point-like vertex coupling the gluons to the Higgs boson. Higher-order corrections in this limit can thus be included into the Wilson coefficient in front of the dimension-five operator in eq. (2.5).
In the remainder of this paper we give a detailed account of all the ingredients that enter our best prediction for the inclusive gluon-fusion cross-section. Furthermore, we carefully analyze the residual uncertainty associated with all of these contributions. In this way we obtain the most precise theoretical prediction for the Higgs production cross-section available to date.
We conclude this section by summarizing, for later convenience, the default values of the input parameters and the concrete choices for PDFs and quark-mass schemes used in our numerical studies. In particular, we investigate three different setups, which are summarized in Tab. 1-3. Note that we use NNLO PDFs even when we refer to lower order terms of the cross-section, unless stated otherwise. The values for the quark masses used are in accordance with the recommendations of the Higgs Cross Section Working Group [86], wherein the top-quark mass was selected to facilitate comparisons with existing experimental analyses at LHC, Run 1 1 .
The cross-section through N 3 LO in the infinite top-quark limit 3.1 The partonic cross-section at N 3 LO in the heavy-top limit In this section we discuss the contributionσ ij,EF T in eq. (2.4) from the effective theory where the top quark is infinitely heavy. This contribution can be expanded into a perturbative series in the strong coupling constant, in eq. (2.5), which admits itself a perturbative expansion in the strong coupling [17,18,19], C n a n s .
Here both the coefficients C n and the strong coupling are functions of a common scale µ. At LO in a s only the gluon-gluon initial state contributes, and we have QCD corrections beyond LO are also known. In particular, the perturbative coefficients η (n) ij are known at NLO [4,5,6] and NNLO [22,23,24] in QCD. Recently, also the N 3 LO corrections η (3) ij have been computed [53]. As they are the main new addition in our computation, we briefly review the N 3 LO corrections to the inclusive gluon fusion crosssection in the heavy-top limit in this section.
We follow the notation of ref. [82] and we split the partonic cross-sections into a singular and a regular part, η The singular contribution is precisely the cross-section at threshold, also known as the soft-virtual cross-section. It contains the contributions from purely virtual three-loop corrections as well as from the emission of soft gluons [80,79,81,87,88]. The regular term takes the form of a polynomial in log(1 − z), where the η for m = 5, 4, 3 have been given in closed analytic form in ref. [82]. For m = 2, 1, 0 no closed analytic expression is available in the literature so far (except for the qq channel [75]). In ref. [53] these coefficients were computed as an expansion around threshold to order 30 in z ≡ 1 − z. In Appendix C we present the numerical values for the first 37 coefficients of the expansion, setting the renormalization and factorization scales equal to the Higgs mass and substituting N f = 5 for the number of light quark flavors and N c = 3 for the number of quark colours. Moreover, it was shown in ref. [53] that a truncation of the series at order O(z 5 ) yields a good approximation to the hadronic cross-section. The first few terms of the expansion may be insightful for theoretical studies of perturbative QCD and discovering universality patterns in subleading terms of the soft expansion. We therefore provide the analytic results for the coefficients in the threshold expansion up to O(z 5 ) in Appendix D.
In the rest of this section we study the numerical impact of the N 3 LO corrections to the inclusive gluon fusion cross-section in the heavy-top limit. We start by studying the validity of approximating the cross-section by its threshold expansion and we quantify the uncertainty introduced by truncating the expansion after only a finite number of terms. We then move on and investigate the perturbative stability ofσ ij,EF T by studying the scale variation of the gluon-fusion cross-section at N 3 LO in the heavy-top limit.

Convergence of the threshold expansion at N 3 LO
As parts of the N 3 LO coefficient functions η (3,m),reg ij (z) have not yet been derived in closed analytic form and are only known as truncated series expansions inz, it is important to assess how well these truncated power series approximate the exact result. In other words, we need to establish how well the threshold expansion converges. Indeed, the partonic cross-sectionsσ ij,EF T need to be convoluted with the partonic luminosities, eq. (2.1), and the convolution integrals receive in principle contributions down to values of z τ 10 −4 . Hence, assessing the residual uncertainty due to the truncation of the series is of utmost importance.
In ref. [82,89] a method was introduced to study the convergence of the threshold expansion. We start by casting the hadronic cross-section in the large-m t limit in the form For n = 0, we recover precisely the usual QCD factorization formula. For n = 0, however, eq. (3.6) is a deformed, but equally valid and equivalent, formulation of the usual QCD factorization formula (2.1). Indeed, it is easy to check that the hadronic cross-section σ EF T is independent of the arbitrary parameter n. Expandingσ ij,EF T /z 1+n into a series around  z = 1, however, introduces a dependence on n order by order in the expansion, which only cancels once infinitely many terms in the series are summed up. Hence, if a truncated series is used to evaluateσ ij,EF T /z 1+n , the result will in general depend on n, and we can use the spread of the n-dependence as a quantifier for the convergence of the series. In Fig. 1 we show the N 3 LO contribution to the hadronic cross-section from the gg−channel. We observe that the hadronic cross-section is very stable with respect to the choice of the arbitrary parameter n after the first ∼ 5 terms in the threshold expansion. In ref. [53] we observed a mild growth of the cross-section at high orders of the threshold expansion (see inlay of Fig. 1 in ref. [53]). This is attributed to the presence of log z terms [49] (and for n > −1 also global factors of 1/z n ) which, after threshold expansion and convolution with the parton distributions, yield a small part of the cross-section. In Fig. 2 we show the convergence of the total cross-section, including all partonic channels, for a variety of different values of n, from the 20th term onwards in the expansion. While we observe good apparent convergence for n > −1, there remains a relatively large spread between the different curves for n ≤ −1. The qualitative difference between these two cases can be understood as follows: For n > −1, we absorb additional factors of 1/z into the partonic cross-sections and expand them around z = 1. This may result in a slower convergence of the partonic threshold expansion for small values of z. At the same time, however, the luminosities are multiplied with powers of z which suppress the contribution from the region z ∼ 0 in the convolution (3.6). The net effect is then a fast apparent convergence for n > −1. This has to be contrasted with the case n ≤ −1, where the luminosities are multiplied by factors of 1/z, which enhance the contribution from the region z ∼ 0 in the convolution (3.6). This leads to a slower apparent convergence, at least in the case where only a few terms are taken into account in the threshold expansion. While the spread between the different curves gives a measure for the quality of the convergence of the threshold expansion, we know of no compelling argument why any of this curves should be preferable over others at this order of the expansion. We observe, however, that the different curves agree among each other within a range of 0.1 pb, thereby corroborating our claim that the threshold expansion provides reliable results for the N 3 LO cross-section. In Fig. 3 we plot the N 3 LO corrections for the gg and qg channels 2 , as well as the total inclusive cross-section, as a function of the truncation order (for n = 0). The quark-initiated channels contribute only a small fraction to the inclusive cross-section. The convergence of the threshold expansion for these channels is less rapid than for the dominant gluon-gluon channel. This is better demonstrated in Fig. 4, where we plot the ratio Here, σ X,EF T (N ) denotes the contribution of the partonic channel X to the N 3 LO correction to the hadronic cross-section when computed through O(z N ) in the threshold expansion. N last (equal to 37) is the highest truncation order used in our current computation. Although the convergence of the quark-gluon and the quark channels is rather slow, the 2 We sum of course over all possible quark and anti-quark flavours. total cross-section and the convergence rate of the threshold expansion are dominated by the gluon-gluon channel. This enables us to obtain a reliable estimate of the cross-section for Higgs production via gluon fusion, even though we have only included a finite number of terms in the threshold expansion. We remark, however, that for quark-initiated processes such as Drell-Yan production a computation in closed form will most likely be necessary.
Besides studying the n-dependence of the truncated power series, we have another way to assess the convergence of the expansion. In ref. [82] it was shown that the knowledge of the single-emission contributions at N 3 LO [64,65,63,62] and the three-loop splitting functions [57,58] is sufficient to determine the coefficients η (3,m) ij in the N 3 LO cross-section (3.5) exactly for m = 5, 4, 3. Recently, also the double-emission contribution at one-loop has been computed in closed form [90]. Using a similar analysis as for m = 5, 4, 3 in ref. [82], it has now been possible to determine also the coefficients with m = 2, 1 exactly for all partonic subchannels. As a consequence, we know all the logarithmically-enhanced terms in eq. (3.5) in closed form, and we only need to resort to a truncated threshold expansion for the constant term, m = 0. We can thus study the convergence of the threshold expansion for the coefficients of η   To summarize, we have investigated the convergence of the threshold expansion at N 3 LO using two different methods. Both methods confirm our expectation that the threshold expansion provides a very good approximation to the exact result. The result of our analysis can be quantified by assigning a (conservative) uncertainty estimate to the truncation of the threshold expansion. We assign an uncertainty due to the truncation of the threshold expansion which is as large as 3 . (3.10) The factor 10 is a conservative estimator of the progression of the series beyond the first 37 terms. Note that the complete N 3 LO cross-section appears in the denominator of eq. (3.10), i.e., the uncertainty applies to the complete N 3 LO result, not just the coefficient of a 5 s .

Scale variation at N 3 LO and the omission of N 3 LO effects in parton densities
Having established that the threshold expansion provides a reliable estimate of the N 3 LO cross-section, we proceed to study the dependence of the cross-section on the renormalization and factorization scales µ R and µ F . In Fig. 5 we fix the factorization scale to µ F = m H /2 and vary the renormalization scale. We observe that the perturbative series in the strong coupling converges faster for 3 In the estimate of the various components of the theoretical uncertainty that we carry out in these sections, we always give numerical results for Setup I. When considering different parameters (Higgs mass or collider energy, for example), we re-assess these uncertainties. For example, δ(trunc) increases from 0.11% at 2 TeV to 0.38% at 14 TeV. Table 4: Renormalization scale variation of the cross-section as defined in eq. (3.11). The factorization scale is fixed to µ F = m H /2. small values of the renormalization scale. It is well known that the scale variation is very large at LO and NLO, and it is still significant at NNLO. To emphasize this point, we indicate in Fig. 5 by horizontal lines the range of predictions for the cross-section at each perturbative order when µ R varies in the interval [m H /4, m H ]. This interval seems to capture the characteristic physical scales of the process, as indicated by the convergence pattern of the series. We quantify the renormalization scale variation by looking at the spread around the average value of the cross-section in this interval, i.e., we define and similarly for σ min EF T,k . The results are shown in Tab. 4. Before we move on to study the dependence of the cross-section on the factorization scale, we note that we evolve the strong coupling α s (µ R ) at N 3 LO, and we use NNLO parton densities at all perturbative orders. The scale variation differs quantitatively from the above table and the convergence of the perturbative series is faster than what is displayed in Fig. 5 if one uses LO or NLO PDFs and α s evolution at the corresponding orders.
Let us now turn to the study of the factorization scale dependence of the N 3 LO crosssection. In Fig. 6 we fix the renormalization scale to µ R = m H /2 and we vary the factorization scale. We observe that at all perturbative orders the variation with the factorization scale is much smaller than with the corresponding variation of the renormalization scale. At N 3 LO, the factorization scale dependence is practically constant over a wide range of values of µ F . A comment is in order concerning the self-consistency of the factorization scale variation at N 3 LO. Traditionally, in a LO computation of a hadronic cross-section the partondensities are not taken to be constant, but they are evolved with the one-loop Altarelli-Parisi splitting functions P (0) . Similarly, at NLO and NNLO the P (1) and P (2) corrections to the splitting functions are included. Following this approach, one would be compelled to include the yet unknown P (3) corrections to the splitting functions in the evolution of the parton densities for our N 3 LO Higgs cross-section computation, which is of course not possible at this point. Nevertheless, our N 3 LO computation with corrections only through P (2) in the DGLAP evolution is consistent in fixed-order perturbation theory, since this is the highest-order splitting function term appearing in the mass factorization contributions. Including the P (3) corrections would be merely a phenomenological improvement (which is necessary for LO calculations in order to obtain qualitatively the physical energy dependence of hadronic cross-sections) but it is not formally required. An inconsistency may only arise due to the extraction of the parton densities from data for which there are no N 3 LO predictions. In fact, this problem has already arisen at NNLO where in global fits of parton distributions jet observables are fitted with NLO coefficient functions. When additional processes are computed at N 3 LO, it is expected that the gluon and other parton densities will be extracted with different values. To our understanding, the uncertainties assigned to the parton densities do not presently account for missing higher-order corrections, but merely incorporate the experimental uncertainties of the data from which they were extracted.
To assess this uncertainty we resort to the experience from the previous orders and present in Fig. 7 the NNLO gluon-fusion cross-section using either NNLO or NLO parton densities as a function of the factorization scale (for a fixed renormalization scale). We notice that the shape of the two predictions is very similar, indicating that differences in the evolution kernels of the DGLAP equation beyond NLO have a small impact. However, in the mass range [m H /4, m H ] the NNLO cross-section decreases by about 2.2 − 2.4% when NNLO PDFs are used instead of NLO PDFs. We can attribute this shift mostly to differences in the extraction of the parameterization of the parton densities at NLO and NNLO. Similarly, we can expect a shift to occur when the N 3 LO cross-section gets evaluated in the future with N 3 LO parton densities rather than the currently available NNLO sets. The magnitude of the potential shift will be determined from the magnitude of the unknown N 3 LO corrections in standard candle cross-sections used in the extraction of parton densities. Given that N 3 LO corrections are expected to be milder in general than their counterparts at NNLO, we anticipate that they will induce a smaller shift than what we observe in Fig. 7. Based on these considerations, we assign a conservative uncertainty estimate due to missing higher orders in the extraction of the parton densities obtained as 4 δ(PDF − TH) = 1 2 where σ (2),(N )N LO EF T denotes the NNLO cross-section evaluated with (N)NLO PDFs at the central scale µ F = µ R = m H /2. In the above, we assumed conservatively that the size of the N 3 LO corrections is about half of the corresponding NNLO corrections. This estimate is supported by the magnitude of the third-order corrections to the coefficient functions for deep inelastic scattering [92] and a related gluonic scattering process [93], which are the only two coefficient functions that were computed previously to this level of accuracy.
So far we have only studied the scale variation from varying µ F and µ R separately. The separation into a renormalization and factorization scale is to a certain extent conventional and somewhat artificial. Indeed, only one regulator and one common scale is required for  Table 5: Scale variation of the cross-section as defined in eq. (3.11) for a common renormalization and factorization scale µ = µ F = µ R . the treatment of both infrared and ultraviolet singularities. For a physical process such as inclusive Higgs production, where one cannot identify very disparate physical scales, large separations between the renormalization from the factorization scale entail the risk of introducing unnecessarily large logarithms. In Fig. 8 we present the dependence of the cross-section on a common renormalization and factorization scale µ = µ R = µ F . Through N 3 LO, the behaviour is very close to the scale-variation pattern observed when varying only the renormalization scale with the factorization scale held fixed. More precisely, using the same quantifier as introduced in eq. (3.11) for the variation of the renormalization scale only, the variation of the cross-section in the range [m H /4, m H ] for the common scale µ is shown in Tab. 5. We observe that the scale variation with µ R = µ F is slightly reduced compared to varying only the renormalization scale at NLO and NNLO, and this difference becomes indeed imperceptible at N 3 LO.
The scale variation is the main tool for estimating the theoretical uncertainty of a cross-section in perturbative QCD, and it has been successfully applied to a multitude of processes. However, in Higgs production via gluon fusion it underestimates the uncertainty both at LO and NLO. It is therefore a critical question to assess whether the scale variation uncertainty is a reliable estimate of the true uncertainty due to missing higher orders in perturbative QCD. We believe that this is most likely the case, because, at least for natural choices of the scales in the interval [m H /4, m H ], the N 3 LO cross-section takes values within the corresponding range of cross-section values at NNLO. Therefore, the progression of the perturbative series from NNLO to N 3 LO corroborates the uncertainty obtained by the scale variation. Indeed, for the central scale µ = m H /2 the N 3 LO cross-section is only ∼ 3.1% higher than at NNLO, i.e., the shift from NNLO to N 3 LO is of the same size as the scale variation uncertainty at N 3 LO. We will therefore take the scale variation uncertainty as our uncertainty estimate for missing higher-order QCD corrections at N 4 LO and beyond. In Section 4, we will also discuss the effect of missing higher orders through resummation methods. This will give additional support to our claim that the scale variation at N 3 LO provides a reliable estimate of missing higher orders beyond N 3 LO.
So far we have only discussed the scale variation for the total hadronic cross-section. It is also interesting and instructive to analyze the scale dependence of the cross-section for individual partonic channels. In Fig. 9 we present the scale dependence at N 3 LO of the gluon-gluon channel, the quark-gluon channel and the total cross-section. The quark-quark and quark-antiquark channels are very small and are not shown explicitly in the plot. We see that, while the gluon-gluon channel dominates over the quark-gluon channel, the latter is important in stabilizing the scale dependence of the total cross-section. Indeed, with the exception of extremely small values of µ, the quark-gluon channel has the opposite slope as the gluon-gluon channel, and therefore a somewhat larger scale variation of the gluon-gluon channel is getting cancelled in the total cross-section. This behaviour can be qualitatively understood from the fact that a change in the factorization scale modifies the resolution on quark-gluon splitting processes, therefore turning quarks into gluons and vice versa. We remark that this feature is not captured by approximate predictions of the cross-section based on the soft-approximations, which only include the gluon-gluon channel.
To summarize, we have identified in this section two sources of uncertainty for the N 3 LO cross-section in the limit of infinite top mass. We observe that the dependence on the factorization scale is flat over wide ranges of values of µ F , and the scale variation is dominated by the µ R variation. Moreover, we see that the inclusion of the quark-gluon channel plays an important role in stabilizing the scale dependence at N 3 LO. Our scale variation estimate of the uncertainty is 1.9% (according to our prescription in eq. (3.11)). We believe that at this order in perturbation theory this uncertainty gives a reliable estimate of missing higher-order corrections from N 4 LO and beyond. In the next section we give further support to this claim by analyzing the effect of various resummations beyond N 3 LO.

Corrections at N 4 LO and beyond in the infinite top-quark limit
In the previous section we have argued that the scale variation at N 3 LO gives a reliable estimate for the missing higher-order corrections to the hadronic gluon-fusion cross-section. In this section we corroborate this claim by investigating various other sources of terms beyond N 3 LO. We check that, if we restrict the analysis to the natural choice of scales from the interval [m H /4, m H ], the phenomenological effect of these terms is always captured by the scale variation at N 3 LO. We start by investigating higher-order terms generated by using an alternative prescription to include the Wilson coefficient C into a perturbative computation, and we turn to the study of higher-order effects due to resummation in subsequent sections. We note at this point that the effect of missing higher-order terms beyond N 3 LO was already investigated in ref. [52] by analysing the numerical impact of the leading N 4 LO threshold logarithms. The conclusions of ref. [52] are consistent with the findings in this section.

Factorization of the Wilson coefficient
The (partonic) cross-section in the effective theory is obtained by multiplying (the square of) the Wilson coefficient by the perturbative expansion of the coefficient functions η ij , see eq. (3.1). As the Wilson coefficient itself admits a perturbative expansion, eq. (3.2), eq. (3.1) takes the following form up to N 3 LO in perturbation theory, where σ 0 denotes the Born cross-section. Conventionally in fixed-order perturbation theory through N 3 LO, one only includes corrections up to O(a 5 s ) from the product in eq. (4.1) and drops all terms of higher order (the Born cross-section is proportional to a 2 s ). This is also the approach adopted in Section 3, where consistently only terms up to O(a 5 s ) had been included. In this section we analyze how the cross-section changes if all the terms shown in eq. (4.1) are included. In this way, we obviously include terms into our prediction that are beyond the reach of our fixed-order N 3 LO computation. We stress that the inclusion of these terms does not spoil the formal N 3 LO accuracy. They can lead, however, to sizeable effects that can be used as a quantifier for missing higher-order terms.
In Fig. 10 we show the value of the cross-section as a function of a common renormalization and factorization scale µ, obtained by either truncating the full cross-section (solid) or by multiplying the truncated Wilson coefficient and truncated coefficient function (dashed). We stress that the difference between the two curves stems entirely from terms at N 4 LO and beyond. However, we observe that if the scale is chosen to lie in our preferred range µ ∈ [m H /4, m H ], then the two curves agree within the scale uncertainty at fixed order N 3 LO, and the difference between the two cross-section values is always well below 2%. In particular, the two curves intersect for µ m H /2. Hence, if we choose the scale µ in the range [m H /4, m H ], both approaches give phenomenologically equivalent answers, and higher-order terms generated by the factorization of the Wilson coefficient only have a very mild phenomenological impact, which is captured by the fixed-order scale variation. We stress that this also supports the claim in Section 3 that the scale variation at N 3 LO gives a reliable estimate of missing higher-order terms in perturbation theory.

Threshold resummation in Mellin space
Fixed-order computations beyond N 3 LO are currently beyond our technical capabilities. Nevertheless, we can get some information on corrections at N 4 LO and beyond from resummation formulae, which allow one to resum certain logarithmically-enhanced terms to all orders in perturbation theory.
In this section we look in particular at higher-order corrections generated by the resummation of threshold logarithms in Mellin space. Before studying the phenomenology of the resummed inclusive Higgs cross-section at N 3 LO+N 3 LL, we give a short review of the formalism that allows one to resum large threshold logarithms in Mellin space.
The Mellin transform of the hadronic cross-section with respect to τ = m 2 H /S is In the following, we always work in the effective theory with an infinitely-heavy top quark. Since the Mellin transform maps a convolution of the type (2.3) to an ordinary product, the QCD factorization formula (2.1) takes a particularly simple form in Mellin space, with the Mellin moments where we suppressed the dependence of the PDFs and the partonic cross-sections on the scales. The Mellin transform is invertible, and its inverse is given by where the contour of integration is chosen such that it lies to the right of all possible singularities of the Mellin moments in the complex N plane. From the definition of the Mellin transform it is apparent that the limit z → 1 of the partonic cross-sections corresponds to the limit N → ∞ of the Mellin moments ofσ ij (N ). In the limit N → ∞ the partonic cross-section in Mellin space can be written as [35] where σ 0 denotes the LO cross-section in the large-m t limit andσ res (N ) is related to the Mellin transform of the soft-virtual cross-section. The constant and logarithmically-divergent contributions in the limit N → ∞ can be written in terms of an all-order resummation formula [25,35,94,95], where the function C gg contains all contributions that are constant for N → ∞. The function G H exponentiates the large logarithmic contributions to all orders and can be written as where β 0 denotes the LO coefficient in the QCD β function. The functions g (n) H are known exactly up to NNLL accuracy [28], i.e., up to g H , which requires knowledge of the cusp anomalous dimension up to three loops [35,87]. In order to perform resummation at N 3 LL accuracy [29,36], the function g (4) H is needed. This function depends on the four-loop cusp anomalous dimension, which is not yet known in QCD. We employ the Padé approximation of ref. [28] for the four-loop cusp anomalous dimensions to obtain a numerical estimate for g (4) H . The numerical impact of this approximation has been studied, e.g., in ref. [28] and we checked that by varying the Padé approximation up and down by a factor of 10, our results do not change 5 .
Let us now turn to the phenomenological implications of resummation at N 3 LL. We obtain N 3 LO + N 3 LL predictions for the cross-section by matching the resummation formula (4.7) to the fixed-order N 3 LO cross-section, i.e., by subtracting from eq. (4.7) its expansion through O(a 5 s ). In this way we make sure that the resummation only starts at O(a 6 s ), which is beyond the reach of our fixed-order calculation. We present our numerical method to perform the inverse Mellin transform (4.5) in Appendix B. In Fig. 11 we show the scale dependence of the resummed cross-section in comparison to the fixed-order crosssection. We see that the resummation stabilizes the scale dependence of the cross-section in comparison to the fixed-order result. At N 3 LO+N 3 LL, the value of the cross-section is essentially independent of the scale choice, and roughly equal to the value of the fixedorder cross-section at µ ≡ µ F = µ R = m H /2. In particular, at µ = m H /2 the effect of the resummation on the N 3 LO cross-section is completely negligible, and in the range µ ∈ [m H /4, m H ], the effect of the resummation is captured by the scale uncertainty at fixed order (albeit at the upper end of the uncertainty band). Hence, the fixed-order result at N 3 LO for µ ∈ [m H /4, m H ] contains the value of the cross-section at N 3 LO+N 3 LL. This corroborates our claim made at the end of Section 3 that at N 3 LO the scale uncertainty provides a reliable estimate of higher-order corrections at N 4 LO and beyond.
We conclude this section by studying the impact of changing the prescription of which terms are exponentiated in Mellin space. The resummation formula eq. (4.7) exponentiates the large-N limit of the fixed-order cross-section eq. (4.6), and as such it is only defined up to subleading terms in this limit. It is therefore possible to construct different resummation schemes that formally agree in the limit N → ∞, but that differ by terms that are suppressed by 1/N . In particular, we may change the exponent in eq. All these schemes are formally equivalent resummation schemes, because they agree in the large-N limit. However, the formally subleading corrections can have a significant numerical impact. In Fig. 12 we show the cross-section predictions for the four different resummation schemes discussed in this section. We observe that within our preferred range of scales, µ ∈ [m H /4, m H ], all four schemes considered in this paper give results that agree within the fixed-order scale variation at N 3 LO, giving further support to our claim that the scale variation at N 3 LO provides a reliable estimate of the remaining missing perturbative orders. We note, however, that outside this range of scales the different prescriptions may differ widely, and we know of no compelling argument why any one of these schemes should be more correct or reliable than the others. Based on these two observations, we are led to conclude that threshold resummation does not modify our result beyond its nominal theory error interval over the fixed-order N 3 LO prediction when the scales are chosen in the range [m H /4, m H ], and we will therefore not include the effects of threshold resummation in Mellin space into our final cross-section prediction.

Threshold and π 2 -resummation in Soft-Collinear Effective Theory
In this section we discuss an alternative way to represent the soft-virtual cross-section in Higgs production, based on ideas from Soft-Collinear Effective Theory (SCET) [30,31,96,97,98]. Just like in the case of threshold resummation in Mellin space, we start by introducing the necessary terminology and review the main ideas, in particular the resummation formula of ref. [37,99,100]. At the end, we combine the N 3 LO coefficient functions with the SCET resummation and study its phenomenological impact. In our analysis we closely follow ref. [99] (for a pedagogical review see ref. [34]). For a comparison to Mellin space resummation see for example refs. [38,39].
Using SCET factorization theorems, the partonic cross-sections at threshold can be factorized into a product of a hard function H, a soft functionS and the effective theory Wilson coefficient C, multiplied by the Born cross-section σ 0 . SCET provides a field theoretical description of these individual functions. They arise when effective field theory is systematically applied and degrees of freedom corresponding to various energy scales are integrated out. The individual coefficients are defined at the respective energy scales, but the total cross-section is independent of these scales. The idea of the SCET formalism is to exploit the factorization of degrees of freedom to derive and solve an evolution equation for each coefficient in the cross-section separately. Consequently, one solves the renormalization group equation for the hadronic cross-section, with the explicit aim to cure the dependence of the cross-section on the various scales. We refer to the scheme outlined above as SCET resummation.
A formula that achieves the aforementioned goals has been derived in ref. [37,99]. It reads, Here, C A = N c is the quadratic Casimir of the adjoint representation of SU (N c ). The hard function was computed through fourth order in refs. [37,60] (in particular, see ref. [60] eq. (7.6), (7.7) and (7.9)). The soft function was recently computed through N 3 LO in ref. [81,101]. We recomputed the soft function up to N 3 LO based on the soft-virtual Higgs cross-section at N 3 LO of ref. [79], and we confirm the result of ref. [81,101]. The definition of A γ and S γ are given for example in ref. [37]. γ cusp is the cusp anomalous dimension [28,102,103,104]. The anomalous dimension γ V can be extracted from the QCD form factor [60] and γ g corresponds to the coefficient of δ(1 − z) of the g → g splitting function [104] . In the previous expression, the soft scale µ s , hard scale µ h and topquark scale µ t are the energy scales of the soft function, the hard function and the Wilson coefficient. The function U mediates the evolution of the individual coefficient functions to the common perturbative scale µ. We note that if we choose the scales according to , and eq. (4.9) corresponds to the fixed-order soft-virtual cross-section. More precisely, if we expand the product of the soft and hard functions and the Wilson coefficient through order a n s , then we reproduce the fixed-order soft-virtual cross-section at N n LO up to terms that vanish in the threshold limit.
Next, let us discuss the choice of the values of the scales µ t , µ h and µ s . First, it is natural to choose the top-quark scale µ t to be the top-quark mass, because this it is the only other mass scale in the Wilson coefficient. Similarly, it seems natural to choose the Higgs mass to be the hard scale entering the hard function H. The cross-section depends on the hard function via its modulus squared, and the hard function depends on the Higgs boson mass via logarithms of the type log we have to analytically continue the logarithms, which then give rise terms proportional to π 2 , In ref. [99] it was observed that at NLO the π 2 term is responsible for a large part of the perturbative corrections at this order. It was suggested to analytically continue the hard function to the space-like region by choosing µ 2 h = −m 2 H . In this approach no π 2 is produced by the analytic continuation of the fixed-order hard coefficient, and the analytic continuation from the space-like to the time-like region is performed in the exponential of eq. (4.10). In the following we also adopt this procedure, which is sometimes referred to as π 2 -resummation, and we choose the hard scale as µ 2 h = −m 2 H . Finally, we have to make a suitable choice for the soft scale µ s . In ref. [97,99] two specific choices were outlined.
1. µ I : The value of µ s where the contribution of the second-order coefficient of the soft function to the hadronic cross-section drops below 15% of the leading order coefficient.
2. µ II : The value of µ s that minimizes the contribution of the second-order coefficient of the soft function to the hadronic cross-section.
Both of the above choices depend on m H and µ, and following ref. [37,97] we choose the average of both scales, We have now all the ingredients to study the phenomenological implications of the SCET resummation. We have implemented the resummation formula, eq. (4.9), into a C++ code and combined it with fixed-order cross-section through N 3 LO. We write the full SCET-resummed cross-section aŝ The above formula matches the resummation to our fixed-order cross-section at N 3 LO such as not to spoil our fixed-order accuracy, and resummation effects only start contributing from N 4 LO. In our implementation we followed closely the public code RGHiggs [37,99,100] and validated our results by comparison. Unfortunately, not all anomalous dimensions required for the evolution of the N 3 LO cross-section are known at this point. We therefore truncate all anomalous dimensions at the maximally known order. Note that already at NNLO the unknown four-loop cusp anomalous dimension would be required. We checked that the numerical dependence of the result on the four-loop cusp anomalous dimension is small and insignificant for phenomenological purposes.
In Fig. 13 we show the hadronic cross-section as a function of a common scale µ = µ R = µ F . We observe that at lower orders there are significant differences between fixed-order and SCET-resummed cross-sections. At N 3 LO, the scale dependence of the resummed cross-section is flat over a wide range of scales. The dependence of the SCET-resummed cross-section on unphysical scales is reduced overall. This can be regarded as a means to find an optimal central value for our prediction. Comparing fixed-order and SCET-resummed cross-section predictions at N 3 LO we find perfect agreement for µ = m H /2, which supports our preferred choice for the central scale. The upward bound of the uncertainty interval obtained by means of scale variation is comparable to the one obtained for the fixed-order cross-section. The lower bound of SCET-resummed cross-section scale variation interval is well contained within the fixed-order interval.
To conclude the analysis, we also need to assess the stability of our result under a variation of the soft, hard and top scales. We do this by varying these scales independently. scales is of the order of ±0.1% as noted already in ref. [37]. As the derived uncertainty intervals and the central values of the SCET-resummed and fixed-order cross-sections are in very good agreement, we will not consider the SCET-resummed cross-section in subsequent chapters.
Let us conclude this section by commenting on the validity of using the π 2 -resummation to predict constant terms at higher orders. Indeed, the exponentiation of the π 2 terms makes a prediction for terms proportional to powers of π 2 , and it is of course interesting to see if these terms capture the bulk of the hard corrections not only at NLO, but also beyond. In particular, we can compare the numerical size of the constant term at N 3 LO predicted by the exponentiation of π 2 to the exact soft-virtual cross-section at N 3 LO of ref. [79]. Since we are interested in fixed-order predictions, we start from eq. (4.9) and we choose the scales according to µ = µ t = µ s = m H and µ 2 h = −m 2 H . Note that choosing µ 2 h < 0 amounts to exponentiating π 2 terms to all orders. Next, let us assume that we know the hard and soft functions to some order in perturbation theory, say through O(a n s ), and all anomalous dimensions governing the evolution equations to one order higher than required to obtain a result that is correct through order n. If we expand the SCET-resummed crosssection in perturbation theory, then we will reproduce the exact soft-virtual cross-section through O(a n s ). By expanding the SCET-resummed cross-section to one order higher we obtain a prediction of terms proportional to powers of π 2 . We want to assess the quality of this prediction by comparing it to the known values of soft-virtual cross-section at low orders. For example, before the coefficient of δ(1 − z) at N 3 LO was computed, all plusdistribution terms of the soft-virtual cross-section at N 3 LO were already known [87,88]. We could thus have made a prediction for the coefficient of δ(1 − z) at N 3 LO based on π 2 resummation. If we denote by C (n) δ the coefficient of the distribution δ(1 − z) in the partonic soft-virtual cross-section accurate through O(a n s ), and where the term proportional to a n+1 s was obtained from the exponentiation of π 2 , we obtain the following sequence of predictions: with a s ≡ a s (m 2 H ). In the previous expressions the Wilson coefficient was set to unity and the number of colours and light flavours are N c = 3 and N f = 5 respectively, and we truncated all numerical results after two digits. We observe that, in the scenario where only the LO cross-section is known, we are able to predict the order of magnitude of the NLO correction, and this prediction would indeed suggest large corrections at NLO. At higher orders, however, the quality of the predictions deteriorates, and in C (2) δ even the prediction of the sign of the N 3 LO correction is wrong. Even if we include the coefficients of the other distributions contributing to the soft-virtual cross-sections, we observe a similar unsatisfactory pattern. We conclude that π 2 terms originating from the systematic exponentiation of the analytic continuation of the hard function constitute only one source of large perturbative corrections to the Higgs boson cross-section, and hence on its own this procedure of predicting higher orders does not provide reliable estimates of the missing dominant corrections.

Quark-mass effects
So far we have only considered QCD corrections to the effective theory where the top quark is infinitely heavy. In this section we discuss effects that are not captured by the effective theory, but that can still give rise to sizeable contributions. In particular, we discuss the inclusion of quark-mass effects from top, bottom and charm quarks, to the extent that these corrections are available in the literature. We start by discussing the effect of quark masses at LO and NLO, where it is possible to obtain exact results including all quark-mass effects. In order to stress the importance of including these effects, we remind that the cross-section changes by +6.3% already at LO if the exact top-mass dependence is taken into account. The exact mass dependence of the cross-section is also known at NLO [5,6,7,8,9,10,11,12,40], and we can thus include all effects from top, bottom and charm quarks up to that order. The value of the cross-section through NLO as we add quark-mass effects for the parameters of Setup 1 (cf. Tab. 1) is summarized in Tab. 6. Beyond NLO finite quark-mass effects are in general unknown, and they can at best be included in an approximate fashion.

pb
Let us start by analyzing finite top-mass effects. The exact NLO cross-section is approximated well by rescaling the EFT cross-section at NLO by the leading-order ratio R LO defined in eq. (2.6). For example, within Setup 1 we have R LO = 1.063, and we see from Tab. 6 that the rescaled NLO cross-section in the effective theory, R LO σ N LO EF T , reproduces the NLO cross-section σ N LO ex;t with full top-mass dependence within 0.65%. Because of this, it has become standard to multiply the EFT cross-section at NNLO by R LO , and we follow this prescription also for the N 3 LO coefficient.
In addition to this rescaling, in ref. [20,21,44,45] top-mass corrections at NNLO were computed as an expansion in m H /m t , after factorizing the exact LO cross-section. We include these corrections into our prediction via the term δ tσ N N LO ij,EF T in eq. (2.4). In particular, we include the contribution from the subleading 1/m t terms for the numerically significant gg and qg channels [45]. The gg channel increases the rescaled EFT crosssection at NNLO by roughly +0.8%, while the qg channel leads to a negative contribution of −0.1%, so that the total net effect is of the order of +0.7%. Note that the small size of these effects corroborates the hypothesis that the cross-section in the effective theory rescaled by R LO gives a very good approximation of the exact result.
Despite the fact that the approximation is good, these contributions come with an uncertainty of their own: the 1/m t expansion is in fact an expansion in s/m 2 t , and consequently it needs to be matched to the high-energy limit of the cross-section, known to leading logarithmic accuracy from k t -factorization. The high-energy limit corresponds to the contribution from small values of z to the convolution integral in eq. (2.1). Since this region is suppressed by the luminosity, a lack of knowledge of the precise matching term is not disastrous and induces an uncertainty of roughly 1%, which is of the order of magnitude of the net contribution. In conclusion, following the analysis of ref. [45], whose conclusions were confirmed by ref. [21], we assign an overall uncertainty of 1% due to the unknown top-quark effects at NNLO.
So far we have only discussed the effect of including top mass effects at NNLO. Despite their suppressed Yukawa couplings, the bottom and charm quarks also contribute to the Higgs cross-section, mainly through interference with the top quark. Indeed, we can easily see from Tab. 6 that the inclusion of bottom-quark effects at LO and NLO leads to sizeable negative contributions to the cross-section, and hence it is not unreasonable to expect this trend to continue at NNLO. Unlike the case of the top quark, however, the contributions of the bottom and charm quarks at NNLO are entirely unknown. We estimate the uncertainty of the missing interference between the top and light quarks within the MS as: With respect to the NNLO cross-section with the exact top effects described in the previous paragraph, this uncertainty is at the level of 0.6%, but it becomes slightly larger at lower energies. For example, at a 2 TeV proton-proton collider it increases to 1.1%. So far, we have assumed that all quark masses are given in the MS-scheme. We now analyze how our predictions are affected if we use the on-shell (OS) scheme. In Tab. 7 we summarize the values of the NLO cross-sections with the quark masses of Setup 1 (MS) and Setup 2 (OS) for a common scale choice µ F = µ R = m H /2. Moreover, the ratio R LO as well as the Wilson coefficient multiplying the cross-section are functions of the top mass, and so they are affected by the choice of the renormalization scheme.
First, let us comment on the use of the OS-scheme for the top-quark mass on the Wilson coefficient. The analytic expression for the Wilson coefficient in the two schemes is the same through NNLO but differs at N 3 LO (see Appendix A). However, this difference is compensated by the different values of the top-quark mass in the two schemes and the numerical value of the Wilson coefficient in the two schemes at N 3 LO agrees to better than a per mille (see penultimate line of Tab. 7). Next, let us turn to the scheme-dependence of R LO . For the top mass of Setup 1 (MS), the value of this ratio is R LO = 1.063, while for the top mass of Setup 2 (OS), we find R LO = 1.066, i.e., the scheme dependence of the rescaled EFT prediction is at the level of 0.3%. Figure 14: The dependence of the cross-section on a common renormalization and factorization scale µ = µ F = µ R in the EFT vs the EFT rescaled with the exact LO contribution in the MSscheme.
Since the top mass runs in the MS-scheme, the LO cross-section acquires its own scale dependence through the dependence of the top mass on the renormalization scale. In Fig. 14 we compare the two approximations as a function of the renormalization schale µ in the MS−scheme. We observe that the scale variation of the rescaled-EFT cross-section is slightly smaller. The variation of the rescaled N 3 LO cross-section in the scale range µ ∈ [ m H 4 , m H ] is ±1.3% (compared to ±1.9% in the pure EFT, cf. Section 3.3). Note that in the OS-scheme the scale uncertainty is the same for the rescaled and pure EFT cross-sections, because the ratio R LO is a constant in this scheme.
The largest scheme dependence appears at LO and NLO due to the non-negligible interference between top and light quarks (see Tab. 7). At LO, the results for the crosssection in the two schemes are in excellent agreement. However, including bottom and charm quark loops gives rise to substantial differences, which at LO are as large as −6.9%. While the difference between the two schemes is reduced to −2.1% at NLO, it still remains larger than the uncertainty estimate of eq. (5.1).
From Tab. 8 it becomes evident that the difference between the results in the two schemes originates from the light-quark contributions. The first line of Tab. 8 shows that, if we only include mass effects from the top quark through NLO, then the results in both schemes are in perfect agreement. Fortunately, light-quark contributions are suppressed in the Standard Model in comparison to the pure top-quark contributions. Indeed, if we set the top-quark Yukawa coupling to zero and only include contributions from bottom and charm quarks (see third line from the bottom of Tab. 8), we observe that the NLO crosssection in the MS scheme is only about a third compared to its value in the OS-scheme. Similarly, the value of the cross-section changes by an order of magnitude between the two  renormalization schemes if only the charm-quark loop is included and both bottom and top-quark Yukawa couplings are set to zero. The very large size of the differences between the two schemes shed some doubt on how well we control the perturbative corrections due to light-quark masses even at NLO. Nonetheless, it is at least encouraging that the scheme dependence of the contributions from bottom and charm quarks alone is significantly smaller at NLO than at LO. Moreover, since the cross-section is dominated by the top quark, the overall scheme dependence due to the light quarks is significantly reduced when the top-quark Yukawa coupling is set to its physical coupling. For our purposes, the most interesting contribution is σ t+b+c − σ t in the last line of Tab. 8, which is the difference between the exact cross-section and the cross-section when all Yukawa couplings, except for the Yukawa coupling of the top quark, are set to zero. This part of the cross-section is only known through NLO and is not captured (at least not in any direct or trustworthy way) by existing NNLO computations 6 . We observe that the NLO K-factor of this contribution is smaller than the NLO K-factor of pure topquark contributions in the cross-section. Therefore, we anticipate that the estimate of the magnitude of the σ t+b+c − σ t correction at NNLO, based on the size of the top-only NNLO K-factor in eq. (5.1), is a conservative estimate within the MS-scheme. However, as we notice from the value of R N LO scheme , there is a scheme dependence of ∼ 30% at NLO. Our preferred scheme is the MS-scheme due to the bad convergence of the perturbative series for the conversion from an MS mass to a pole mass for the bottom and charm quarks [106,107]. To account for the difference with the OS scheme, we enlarge the uncertainty on σ t+b+c −σ t , as estimated via eq. (5.1) within the MS scheme, by multiplying it with a factor of 1.3, Let us conclude this section by commenting on the amount by which the cross-section changes when the values of the quark masses used as input vary from those of Setup 1. As argued in the previous section, the dependence on the rescaled EFT cross-section on the top-quark mass is extremely mild. We will therefore focus in this section on the exact QCD corrections (including the light quarks) through NLO, and we study the variation of the cross-section when the quark masses are varied following the internal note of the HXSWG [86], which either conforms to the PDG recommendation or is more conservative (see Tab. 9). We see that the parametric uncertainties are entirely negligible, at the level of 0.1% or below. Finally, the parametric uncertainty on the ration R LO does not exceed 0.1%. For this reason, we will not consider parametric uncertainties on quark masses any further.

Electroweak corrections
So far we have only considered higher-order QCD corrections to the gluon fusion crosssection. However, in order to obtain precise predictions for the Higgs cross-section also electroweak (EW) corrections need to be taken into account. The EW corrections to the LO gluon fusion cross-section have been computed in ref. [40,41,42]. For a Higgs mass of m H = 125 GeV, they increase the LO cross-section by 5.2%, and we take these corrections into account in our cross-section prediction.
Given the large size of the NLO QCD corrections to the Higgs cross-section, we may expect that also the EW corrections to the NLO QCD cross-section cannot be neglected. Unfortunately, these so-called mixed QCD-EW corrections are at present unknown. The contribution from light quarks, which at O(a EW a 2 s ) is the dominant one accounting for ∼ 95% of the total EW corrections at that order, was computed in ref. [43] within an effective field theory approach where the W and Z bosons are assumed heavier than the Higgs boson and have been integrated out. This has the effect of introducing EW corrections to the Wilson coefficient describing the effective coupling of the Higgs boson to the gluons in eq. (2.5), where C QCD encodes the pure QCD corrections to the Wilson coefficient, λ EW denotes the EW corrections to the LO cross-section of ref. [41] and C 1w are the mixed QCD-EW corrections in the EFT approach of ref. [43]. The value of the coefficient C 1w is [43] C 1w = 7 6 .
Adopting the modification of the Wilson coefficient also for higher orders in a s leads to a total correction of 5.0%. We stress that the numerical effect of this correction is very similar to that of the 'complete factorization' approach to include EW corrections of ref. [41], which lead to an increase of the NLO cross-section by 5.1%. The effective theory method for the mixed QCD-EW corrections is of course not entirely satisfactory, because the computation of the EW Wilson coefficient assumes the validity of the m H /m V expansion, V = W, Z while clearly m H > m V . We thus need to carefully assess the uncertainty on the mixed QCD-EW corrections due to the EFT approximation. In the region m H > m V , we expect threshold effects to be important and one should not expect that a naive application of the EFT can give a reliable value for the cross-section. However, in eq. (6.1) the EFT is only used to predict the relative size of QCD radiative corrections with respect to the leading order electroweak corrections. This can only vary mildly above and below threshold. For phenomenological purposes, we expect that the rescaling with the exact λ EW in eq. (6.1) captures the bulk of threshold effects at all perturbative orders. To quantify the remaining uncertainty in this approach, we allow the coefficient C 1w to vary by a factor of 3 around its central value in eq. (6.2). We do this by introducing a rescaling factor y λ by Varying y λ in the range [1/3, 3], we see that the cross-section varies by −0.2% to +0.4%. We summarize the dependence of the cross-section on y λ in Fig. 15. Note that the result obtained by assuming complete factorization of EW and QCD corrections (marked by 'CF' in Fig. 15) lies in the middle of the variation range, slightly higher than the y λ = 1 prediction. Finally, we stress that the choice of the range is largely arbitrary of course. It is worth noting, however, that in order to reach uncertainties of the order of 1%, one needs to enlarge the range to y λ ∈ [−3, 6]. An alternative way to assess the uncertainty on the mixed QCD-EW corrections is to note that the factorization of the EW corrections is exact in the soft and collinear limits of the NLO phase space. The hard contribution, however, might be badly captured. At NLO in QCD, the hard contribution amounts to ∼ 40% of the O(a 3 s ) contribution to the cross-section, where we define the hard contribution as the NLO cross-section minus its soft-virtual contribution, i.e., the NLO contribution that does not arise from the universal exponentiation of soft gluon radiation (see Section 4). In the notation of Section 3 the hard contribution is defined as the convolution of the parton-level quantitŷ with the PDFs, which receive contributions from the gg, qg and qq initial state channels. The mixed QCD-EW corrections are 3.2% of the total cross-section. Even if the uncertainty of the factorization ansatz is taken to be as large as the entire hard contribution, we will obtain an estimate of the uncertainty equal to 0.4 × 3.2% = 1.3% with respect to the total cross-section. An alternative way to define the hard contribution is to look at the real emission cross-section regulated by a subtraction term in the FKS scheme [108]. We could then exclude the contribution of the integrated subtraction term, which is proportional to the Born matrix element, and hence of soft-collinear nature. We would then estimate the hard contribution as ∼ 10% of the O(a 3 s ) contribution to the cross-section, which would lead to an uncertainty equal to 0.1 × 3.2% = 0.32%.
We note that the different estimates of the uncertainty range from 0.2% to 1.3%. We therefore assign, conservatively, an uncertainty of 1% due to mixed QCD-EW corrections for LHC energies. This uncertainty decreases for smaller collider energies as the soft contributions become more important and the factorization ansatz becomes more accurate.
For example, at a 2 TeV proton-proton collider the most conservative estimate of the uncertainty is 0.8%.

PDF comparison
So far we have only discussed perturbative higher-order corrections to the partonic crosssections. The full hadronic cross-section is then obtained by convoluting the partonic coefficient functions by the parton distribution functions. In the last few years significant progress has been made towards the improvement of the PDF fits, also through the inclusion of new data from collider and fixed-target experiments. We refer to the analysis in the latest PDF4LHC working group paper [109] for a review of the updated sets ABM12 [110], CT14 [111], JR14 [112], MMHT2014 [113], NNPDF3.0 [114] and HERAPDF2.0 [115], which are available through NNLO, as well as the NLO set CJ12 [116]. In this Section, we will compare the predictions from various pdf sets using Setup 1 and the partonic cross-sections derived in the rescaled EFT through N 3 LO for a factorisation and renormalisation scale µ = m H /2.
The three sets that enter the PDF4LHC fit (CT14, MMHT14 and NNPDF3.0) and HERAPDF2.0, are provided at the same value of the strong coupling constant as the global PDF4LHC15 combination [109], This value is consistent with the PDG average [117]. In Fig. 16 we compare the 68% C.L. predictions from CT14, MMHT2014 and NNPDF3.0 with those from the PDF4LHC15 combination. For comparison purposes, in this section we combine (potentially asymmetric) PDF and α s uncertainties in quadrature 7 , δ ± (P DF + α s ) = δ ± (P DF ) 2 + δ ± (α s ) 2 .
(7.2) From Fig. 16, we observe that the predictions obtained from the three sets that enter the PDF4LHC15 combination lie well within 1% of each other over the whole range of center-of-mass energies from 2 to 15 TeV. In particular, MMHT2014 and NNPDF3.0 agree at the per mille level. The combined PDF+α s uncertainty is at the level of 3 − 4% for LHC energies, and it captures very well the small differences in the predictions among the different sets.
Good agreement with the PDF4LHC15 predictions is also obtained for LHC energies using the HERAPDF2.0 set (Fig. 17). HERAPDF2.0 does not enter the PDF4LHC fit, but is given at the same central value of α s . However, these PDFs give a cross-section that is about 6% lower at Tevatron energies, and increase above the PDF4LHC15 predictions at higher center-of-mass energies. This value is the result of the ABM fit. As one can see from Fig. 18, the ABM12 set gives a prediction that is about 23% lower than the one from PDF4LHC15 at Tevatron energies, and 9 − 7% lower at LHC energies. The PDF+α s error is 1.2%, which does not account for this discrepancy. We note here that the variation range for α s used for the PDF+α s variation in the ABM12 set is determined by the fitting procedure and is slightly smaller than the range suggested by the PDF4LHC recommendation [109].
To understand how much of this difference comes from the choice of a different value of the strong coupling constant, we plot in Fig. 18 the prediction from CT14 at the same value of α s as the one obtained by ABM12. At α s = 0.118 the predictions from CT14 are in very good agreement with those from PDF4LHC15 (Fig. 16). At a lower value of α s , CT14 gives a cross-section that is about 10% smaller than the result at α s = 0.118 (12% at Tevatron energies). The dependence on the center-of-mass energy appears to be much milder than the one exhibited by ABM12. However, the PDF+α s uncertainty might improve the agreement between the two sets. Unfortunately, only one error set for CT14 at α s = 0.113 is available, and we cannot assess this uncertainty.

Recommendation for the LHC
In previous sections we have considered various effects that contribute to the gluon-fusion Higgs production cross-section at higher orders. In this section we combine all these effects, and as a result we are able to present the most precise prediction for the gluon-fusion crosssection available to date. In particular (for the Setup 1 of Tab. 1) for a Higgs boson with a mass m H = 125 GeV, the cross-section at the LHC with a center-of-mass energy of 13 TeV is σ = 48.58 pb Next, let us analyze the uncertainties quoted in our cross-section prediction. We present our result in eq. (8.1) with two uncertainties which we describe in the following. The first uncertainty in eq. (8.1) is the theory uncertainty related to missing corrections in the perturbative description of the cross-section. Just like for the central value, it is interesting to look at the breakdown of how the different effects build up the final number. Collecting all the uncertainties described in previous sections, we find the following components: In the previous table, δ(scale) and δ(trunc) denote the scale and truncation uncertainties on the rEFT cross-section, and δ(PDF-TH) denotes the uncertainty on the cross-section prediction due to our ignorance of N 3 LO parton densities, cf. Section 3. δ(EW), δ(t, b, c) and δ(1/m t ) denote the uncertainties on the cross-section due to missing quark-mass effects at NNLO and mixed QCD-EW corrections. The first uncertainty in eq. (8.1) is then obtained by adding linearly all these effects. The parametric uncertainty due to the mass values of the top, bottom and charm quarks is at the per mille level, and hence completely negligible. We note that including into our prediction resummation effects in the schemes that we have studied in Section 4 would lead to a very small scale variation, which we believe unrealistic and which we do not expect to capture the uncertainty due to missing higher-order corrections at N 4 LO and beyond. Based on this observation, as well as on the fact that the definition of the resummation scheme may suffer from large ambiguities, we prefer a prudent approach and we adopt to adhere to fixed-order perturbation theory as an estimator of remaining theoretical uncertainty from QCD.
The second uncertainty in eq. (8.1) is the PDF+α s uncertainty due to the determination of the parton distribution functions and the strong coupling constant, following the PDF4LHC recommendation. When studying the correlations with other uncertainties in Monte-Carlo simulations, it is often necessary to separate the PDF and α s uncertainties: Since the δ(α s ) error is asymmetric, in the combination presented in eq. (8.1) we conservatively add in quadrature the largest of the two errors to the PDF error.
As pointed out in Section 7, the PDF4LHC uncertainty estimate quoted above does not cover the cross-section value as predicted by the ABM12 set of parton distribution functions. For comparison we quote here the corresponding cross-section value and PDF+α s uncertainty with the ABM12 set 8 : σ ABM12 = 45.07 pb +2.00 pb (+4.43%) −2.88 pb (−6.39%) (theory) ± 0.52 pb (1.17%) (PDF+α s ) . (8. 3) The significantly lower central value is mostly due to the smaller value of α s , which however is also smaller than the world average.
It is also interesting to compare our prediction (8.1) to the value one would have obtained without the knowledge of the N 3 LO corrections in the rEFT. We find Here we adopt exactly the same prescription to estimate the uncertainty at NNLO and at N 3 LO, and we do not only rely on scale variation for the NNLO uncertainty estimate, as was often done in the past. Finally, we have also studied how our predictions change as we vary the center-of-mass energy and the value of the Higgs mass. Our predictions for different values of the protonproton collision energy and a Higgs mass of m H = 125 GeV are summarized in Tab. 10. In comparison to the official recommendation of the LHC Higgs Cross-section Working Group earlier than our work [48], our results have a larger central value by about 11%. The difference can be attributed to the choice of optimal renormalization and factorization scale, the effect of the N 3 LO corrections, the different sets of parton distribution functions and value of α s as well as smaller differences due to the treatment of finite quark-mass effects. In comparison to the earlier recommendation from some of the authors in ref. [120], our result has a central value which is higher by 3.5%. The difference can be attributed to the effect of the N 3 LO corrections, the different sets of parton distribution functions and value of α s as well as smaller differences due to the treatment of finite quark-mass effects.
Additional cross-section predictions for a variety of collider energies and Higgs boson masses can be found in Appendix E.

Conclusion
In this paper we have presented the most precise prediction for the Higgs boson gluonfusion cross-section at the LHC. In order to achieve this task, we have combined all known    Table 12: Earlier recommendation for the gluon-fusion Higgs cross-section at a proton-proton collider by some of the authors in ref. [120].
higher-order effects from QCD, EW and quark-mass corrections. The main component that made our computation possible was the recent computation of the N 3 LO correction to the cross-section in an effective field theory where the top quark was integrated out.
In an appendix we present analytic expressions for the partonic subchannels of the N 3 LO partonic cross-sections which have not been presented elsewhere in the literature, in the form of a series expansion around the threshold limit. The N 3 LO corrections moderately increase (∼ 3%) the cross-section for renormalization and factorization scales equal to m H /2. In addition, they notably stabilize the scale variation, reducing it almost by a factor of five compared to NNLO. The N 3 LO scalevariation band is included entirely within the NNLO scale-variation band for scales in the interval [m H /4, m H ]. Moreover, we have found good evidence that the N 3 LO scale variation captures the effects of missing higher perturbative orders in the EFT. We base this conclusion on the following observations: First, we observed that expanding in α s separately the Wilson coefficient and matrix-element factors in the cross-section gives results consistent with expanding directly their product through N 3 LO. Second, a traditional threshold resummation in Mellin space up to N 3 LL did not contribute significantly to the cross-section beyond N 3 LO in the range of scales µ ∈ [m H /4, m H ]. Although the effects of threshold resummation are in general sensitive to ambiguities due to subleading terms beyond the soft limit, we found that within our preferred range of scales, several variants of the exponentiation formula gave very similar phenomenological results, which are always consistent with fixed-order perturbation theory. Finally, a soft-gluon and π 2 -resummation using the SCET formalism also gave consistent results with fixed-order perturbation theory at N 3 LO. While ambiguities in subleading soft terms limit the use of soft-gluon resummation as an estimator of higher-order effects, and while it is of course possible that some variant of resummation may yield larger corrections, it is encouraging that this does not happen for the mainstream prescriptions studied here.
Besides studying the effect of QCD corrections in the EFT at high orders, we also investigated the cross-section in the EFT after inclusion of exact LO and NLO QCD corrections in the full Standard Model theory (with finite top, bottom and charm quark masses) and 1/m t corrections at NNLO. We also included known two-loop electroweak corrections and an estimate of three-loop mixed QCD-EW corrections into our final prediction.
No prediction for the cross-section would be complete without estimating the residual uncertainties that may affect our result. We have identified several sources of theoretical uncertainties, namely, the truncation of the threshold expansion, the QCD scale variation, missing higher-order corrections in the extraction of parton densities, missing finite quarkmass effects beyond NLO and missing mixed QCD-EW corrections. After adding all these uncertainties linearly, we obtain a residual theoretical uncertainty of about 5 − 6%. We have also studied the sensitivity of the cross-section on the choice of parton distribution functions. The CT14, MSTW and NNPDF sets are in good agreement among themselves, and have been combined together according to the PDF4LHC recommendation. They yield a combined uncertainty due to both α s and parton densitites of the order of ∼ 3.5%. The PDF4LHC sets give cross-section values that are in good agreement with the cross-section as computed with HERAPDF sets. However, the ABM12 set of parton densities yields results which are significantly lower and outside the quoted range of uncertainty.
We expect that further progress can be made in order to improve even more the precision of our computation. A forthcoming computation of the N 3 LO cross-section in the EFT in a closed analytic form will remove the truncation uncertainty. Future computations of the NNLO QCD cross-section in the full Standard Model (including finite top, bottom and charm masses) and a complete computation of three-loop mixed QCD-EW corrections will remove further significant sources of uncertainties. Progress in the determination of parton densities, with more precise LHC data and more precise computations of crosssections used in the extraction of parton densities, will be crucial to corroborate the PDF + α s uncertainty and to resolve discrepancies due to systematic effects.
To conclude, we have presented the predictions for the Higgs boson cross-section in gluon fusion, based on very high orders in perturbation theory. In this way, we have obtained the most precise prediction of the Higgs boson production cross-section at the LHC to date. We are looking forward to comparisons of our results with precise measurements of the Higgs boson cross-section at the LHC in the future. The analogous result in the on-shell scheme can be derived combining the OS decoupling constant ζ g and the OS Wilson coefficient with α s running in the full theory, that one can find in the literature (see, for example, ref. [17,18]

B. Numerical implementation of the Mellin inversion
In this section we describe our numerical implementation of the inverse Mellin transform, where f g (N ) are the Mellin moments of the gluon density andσ gg (N ) the (resummed) partonic cross-section. The integration contour is the straight vertical line Re(N ) = c, chosen according to the minimal prescription [27], i.e., all the poles inσ gg (N ) lie to the left of the contour, except for the Landau pole in Mellin space, which lies to the right of the contour. The position of the Landau pole in Mellin space is given by .

(B.2)
We parametrize the integration contour as N = c + it, and we obtain In order to evaluate the integral, we need to know the Mellin moments of the gluon density for complex values of the Mellin variable N . To our knowledge, the public PDF sets do not provide grids which allow one to immediately obtain the Mellin moments of the PDFs, and so we need to use our own method to perform the inverse Mellin transform in eq. (B.3). This method is described in the remainder of this appendix. We start by truncating the integral (B.3) at some large value t max , and we approximate the integral over the range [0, t max ] by a Gauss-Legendre quadrature of order l, where t (l) k are the zeroes of the l-th Legendre polynomial P l (x). The Gauss-Legendre weights are given by where a l denotes the coefficient of x l in P l (x). The advantage of eq. (B.4) is that we only need to know the Mellin moments of the gluon density on the finite set of points N (l) k . We can evaluate these Mellin moments numerically once and for all (for a given PDF set and factorization scale) and store them in a grid, Note that the integral (B.6) is numerically convergent for Re N (l) k > 0. Equation (B.4) is our master formula for the computation of inverse Mellin transforms. We have generated grids f g N (l) k for various choices of PDF sets and factorization scales, making it straightforward to compute eq. (B.4) for any of these choices. Let us make some comments about the master formula (B.4). First, we see that the right-hand side of eq. (B.4) depends on three free parameters: the real part c of the integration contour, the truncation t max and the order l of the Gauss-Legendre quadrature. While the inverse Mellin transform must obviously be independent of these parameters, they may introduce some systematic uncertainties. In our implementation we choose c = 2.5, and we checked that the value of the integral remains unchanged under small deformations of this value. Next, it is easy to check that the Mellin moments are highly suppressed for Im(N ) 1. In our implementation we choose t max = 125, and we checked that the contribution to the integral from the range [100,125] is completely negligible. Note that this implies that the bulk of the value of the integral comes from the region where Im(N ) is small. We therefore partition the range [0, t max ] into subregions of increasing length, and in every subregion we approximate the integral by a Gauss-Legendre quadrature of order l = 20. In this way we make sure that the sum in eq. (B.4) receives mostly contributions from points where Im(N ) is small.
Finally, let us briefly comment on the choice of the straight-line contour in the inverse Mellin transform. Indeed, we could deform the contour such as to maximize the convergence of the numerical integration. In particular, in ref. [123] it was argued that the inverse Mellin transform converges faster if the contour is chosen as N = c + t exp(iφ), π/2 < φ < π. In this case, however, we have Re(N ) = c + t cos φ, and so Re(N ) < 0 for large enough t, which contradicts the convergence criterion for the integral (B.6). Hence, as we need to perform the integral (B.6) numerically in our approach, we cannot choose the optimized integration contour of ref. [123]. We note, however, that since it is sufficient to generate the grids f g N (l) k for a sufficiently large number of points once and for all, speed is not an issue and we do not loose anything by choosing φ = π/2.

D. Color and flavor number dependence for the first coefficients of the threshold expansion
In this appendix we present analytic result for the first few coefficients in the threshold expansion for each partonic channel. We use the notation The non-zero coefficients C ij [m, a, b, n] for m = 0, 1, 2 and n ≤ 5 are given in the remainder of this appendix.

E. Cross sections in the HXSWG recommended mass range
We present here the gluon-fusion Higgs production cross-section at a proton-proton collider for center-of-mass energies of 2, 7, 8, 13 and 14 TeV and for a Higgs boson of mass from 120 GeV to 130 GeV. The choice of these parameters follows the indications of the Higgs Cross Section Working Group [86]. The components that enter (linearly) the theory uncertainty have been discussed in the text. To summarize the main points of that discussion, the scale variation uncertainty is assessed at each energy and Higgs mass through a scan over µ ∈ [m H /4, m H ]. The uncertainties due to truncation, unknown N 3 LO PDFs and unknown finite-mass effects are also evaluated every time, following the procedure described in the text. For the missing QCD-EW effects, we find that a reasonable estimate yields a 1% uncertainty. At 2 TeV, however, we adopt δ(EW) = 0.8%, which is the most conservative estimate we obtain over the mass range analyzed. Finally, we assign a 1% uncertainty to missing finite-top mass effects at NNLO [21,45]. We use the PDF set PDF4LHC15.      Table 17: Gluon-fusion Higgs production cross-section at a proton-proton collider for √ s = 14 TeV. Details on the calculation of the theory error are given at the beginning of this Appendix and in the main text.