Scale dependence and collinear subtraction terms for Higgs production in gluon fusion at N3LO

The full, explicit, scale dependence of the inclusive N3LO cross section for single Higgs hadroproduction is obtained by calculating the convolutions of collinear splitting kernels with lower-order partonic cross sections. We provide results for all convolutions of splitting kernels and lower-order partonic cross sections to the order in epsilon needed for the full N3LO computation, as well as their expansions around the soft limit. We also discuss the size of the total scale uncertainty at N3LO that can be anticipated with existing information.


Introduction
During the past year both multi-purpose experiments at CERN's Large hardon collider (LHC), CMS [1] and Atlas [2], have observed a new boson with a mass of about 125 GeV, which is strongly believed to be the long-sought Higgs boson. The couplings of the new boson to Standard Model (SM) particles are currently compatible to the SM predictions with a minimal Higgs sector. Nevertheless, effects from physics beyond the Standard Model (BSM) may reside in small deviations of the couplings from the SM values, effects that will be, to a certain extent, accessible with the increased statistics and energy reach of the LHC high energy run starting in 2015.
The dominant production mode for the Higgs boson at the LHC is gluon fusion, accounting for about 90% of the total production cross section at the observed mass of about 125 GeV. Indeed, the Higgs boson has, up to now, been observed in channels in which its production is gluon induced. Next-to-leading order (NLO) QCD corrections for gluon fusion, in the five-flavour heavy-quark effective theory (HQET) were computed at the beginning of the 1990s [3,4]. Since then NLO corrections in the full theory including top-mass, and top-bottom interference effects were calculated by [5][6][7][8][9], and next-to-next-to-leading order (NNLO) corrections in HQET by [10][11][12]. Electroweak corrections are also available at the NLO level [13][14][15][16], and so are mixed QCD-EW corrections [17], and EW corrections to higgs plus jet including top and bottom quark contributions [18]. Recently the NNLO cross section for gluon induced higgs production in association to a jet was calculated, in a way that also allows for differential distributions to be produced [19]. All available fixed order contributions to Higgs production via gluon fusion were recently included in the program iHixs [20], which, moreover, allows for the incorporation of BSM effects through modified Wilson coefficients within the effective theory approach. The latter has been explicitly shown to be an excellent approximation for Higgs masses below the top-antitop threshold [21][22][23], and even more so for the Higgs boson at 125 GeV. Despite these advances and due to the slow perturbative convergence of the gluon fusion cross section, the remaining uncertainty due to variation in renormalisation and factorisation scale still amounts to about ∼ 9% for a 125 GeV Higgs boson at the LHC with 8 TeV centre-of-mass energy.
Beyond fixed order, threshold resummation has been performed to NNLL accuracy by traditional resummation methods [24] leading to a ∼ 7.5% uncertainty [25], and within the SCET framework [26][27][28] leading to a ∼ 4% scale uncertainty. The latter is generally considered too optimistic.
Information from the LHC high energy and high luminosity data set is projected to allow the determination of the Higgs couplings with precision of ∼ 10% or better [29][30][31]. This uncertainty includes experimental systematics and statistics, but also errors from the determination of parton distribution functions and of the strong coupling, as well as theory systematics, the latter being the limiting factor in several cases. It is evident that a prerequisite to this goal is the reduction of the theory scale uncertainty to the ∼ 5% level or lower. The question arises then, whether computing the cross section to the next order in perturbation theory, N 3 LO, within the EFT approach, an admittedly formidable task, would achieve this goal.
Information about certain N 3 LO contributions has been available for several years. The three-loop, virtual contributions have been calculated and were part of the full N 3 LO Higgs decay to gluons in [32]. However, disentangling the pure virtual contributions from this computation is not possible. The quark and gluon form-factors are known up to threeloop order [33][34][35][36]. In [37] the soft 'plus'-contributions to the N 3 LO cross section were derived using mass factorisation constraints. This allowed the authors of [37] to derive a soft approximation of the N 3 LO cross section whose renormalisation scale dependence is rather mild, resulting in 4% renormalisation scale uncertainty (keeping the factorisation scale equal to the Higgs mass). Recently further attempts to modify the resummation procedure such that its prediction at fixed order better matches the threshold and high energy limits of the known fixed order results, were made [38], resulting in another soft approximant with a scale uncertainty of 7%. It still remains true that without the full N 3 LO expression, it is difficult to judge which of these prescriptions is closer to reality.
Recently, some new ingredients of the full N 3 LO cross section have appeared. In [39,40], the real-virtual and double-real master integrals of the NNLO cross section have been calculated to higher orders in ǫ. In [41], the convolutions of collinear splitting kernels with lower-order partonic cross sections have been computed, which is also an ingredient for our result and has been re-derived in this work. Very recently, the soft limits of all master integrals appearing in triple-real radiation corrections (i.e. the emission of three additional partons) have been worked out [42].
In this paper, we compute the full dependence of the N 3 LO cross section on the factorisation and renormalisation scales, which can be obtained from lower-order results. Furthermore, we provide the soft limits of all convolutions that we calculated, which may become useful when expanding the full N 3 LO corrections around threshold. In section 2 we review how the dependence on factorisation and renormalisation scales enters higher-order calculations. In sections 3 and 4 we list the splitting kernels and partonic cross sections needed for our results and present the method used to compute their convolutions, respectively. In section 5 we give results for the estimated scale uncertainty of the N 3 LO gluon fusion cross section and conclude in section 6.

Sources of explicit scale dependence
Predictions for observable quantities in quantum field theory are independent of arbitrary scales, when calculated at all perturbative orders. The scale dependence of all predictions is an artefact of the truncation of the perturbative series, and is usually considered a measure of the effect of missing higher orders in any given computation. This dependence occurs explicitly, through terms in the final result that depend on logarithms involving the scale, and implicitly, through the running of α s and the evolution of the parton distribution functions. In this section we describe the occurrence of the explicit scale dependence.
Let us, for the moment, introduce only one scale, In dimensional regularisation the scale µ appears during renormalisation, when the bare coupling is replaced by the renormalised one, where we have chosen the MS-scheme. Z α is the renormalisation constant of the strong coupling, and the factor of µ 2ǫ ensures that the coupling and thus the action remain dimensionless in D = 4 − 2ǫ dimensions as well. We define a(µ) ≡ αs(µ) π throughout the paper. Divergences, of UV or IR nature, manifest themselves as poles in the regularisation parameter ǫ. The leading divergences, ǫ −2n , . . . , ǫ −n−1 for the n-th order correction, vanish, among real and virtual contributions, after renormalisation counter terms are included. The remaining poles of the UV-renormalized partonic cross section, starting from ǫ −n , vanish only after subtraction of collinear counterterms.
Specifically, let us denote byσ ij the partonic cross section after renormalisation 1 (which still contains divergences of infrared (IR) origin), spin,col The expansion ofσ ij can be written aŝ where we have written explicitly the pole coefficients at every order in a and the associated logarithms L f ≡ log(µ 2 /s). The relation ofσ ij to the total, inclusive cross section is given by convolution with the parton distribution functions, f i (x), and the collinear counter terms, Γ ij (x), by where summation of repeated indices is implied and τ = m 2 h /S with S the total centre-ofmass energy of the collision.
The convolution is defined by 1 Note the 1/z factor in the definition ofσij that is necessary to make eq. 2.6 work. and the collinear counter term reads The P

(n)
ij are the Altarelli-Parisi splitting kernels which govern the emission of collinear partons (see section 3.2).
Within the renormalized N n LO partonic cross section,σ ij , the logarithmic dependence on the scale µ arises when residual poles of order up to ǫ −n are multiplied with the expansion of the factor µ 2ǫ /s ǫ (the s −ǫ originating from the d-dimensional phase-space measure), These poles are required to cancel against the poles from the collinear counter terms convoluted with lower order partonic cross sections, see eq. (2.6) and (2.8). This requirement fixes the coefficientsσ (n,r) ij for −n − 1 < r < 0 which are also the coefficients of the logarithmic terms. In summary, all contributions to the N n LO cross section that are proportional to a power of log(µ 2 /s) can be obtained from calculating the convolutions of splitting kernels and lower-order, N m<n LO partonic cross sections.
The computation of all combinations of splitting kernels and partonic cross sections relevant for the N 3 LO corrections to the gluon fusion process is the main work of this publication. We calculate also all pieces of higher orders in ǫ that will add to the finite part of the N 3 LO corrections (but not necessarily to the scale dependent parts).
With the pole cancelation achieved, let us now define the finite, mass-renormalised partonic cross section withσ km now explicitly dependent on log(µ 2 /s). Alternatively, the relation above can be inverted,σ and solved for the highest order of σ one is interested in. For example, at NLO, this yields (using Γ This step-by-step procedure provides an additional test on the lower-order cross section, since a mistake in their mass-renormalisation will result in uncancelled poles at a higher order. We provide results in this framework, i.e. our convolutions involve splitting kernels and the finite partonic cross section σ ij . We can, then, set µ = µ f and use the renormalisation group equation for the strong coupling constant, to change the scale at which α s is evaluated in σ ij (µ r , µ f , z).
The third-order expansion of a(µ f ) in terms of a(µ r ) reads with L ≡ log(µ 2 r /µ 2 f ). The explicit logarithm L essentially counters the running of the coupling constant up to the order considered, such that the effect of varying the unphysical scales is weakened for higher orders in perturbation theory, i.e. the perturbative series is converging to its all-order value, as can be seen by taking the total derivative of the partonic cross section, If there are more scale dependent quantities entering the cross section, such as running MS-masses, their scale translations have to be included as well. The full dependence on the scale µ will, by construction, be of the next order in a once again.

Ingredients
From the previous section, we conclude that we need the following ingredients to obtain all convolutions required for the N 3 LO gluon fusion cross section: • The LO partonic cross section through O(ǫ 3 ).
• The LO splitting kernels P qq .
• The NLO splitting kernels P • The NNLO splitting kernels P (2) gg and P (2) gq (owing to the fact that at LO, only the gg-channel is nonzero).

Partonic cross section
We work in the effective five-flavour theory with the top quark integrated out. This approximation has been shown to be very good (less than 5%) for light Higgs masses, as can be seen by comparing the NLO results in effective and full six-flavour theory and by studying the importance of 1/m t corrections of the effective NNLO cross section [22,23]. We expect this behaviour to persist at N 3 LO.
The effective Lagrangian describing the interaction between gluons and the Higgs boson is given by where G a µν denotes the gluonic field-strength tensor. The Wilson coefficient C 1 which starts at O(a) has been computed perturbatively to four-loop accuracy [43,44] in the SM as well as to three-loop accuracy for some BSM models [45][46][47]. Through O(a 4 ), the Wilson coefficient in the SM reads with L t = log(µ 2 /m 2 t ). N F denotes the number of light flavours set to 5. The above expression denotes the renormalised Wilson coefficient, which is related to the bare one through the renormalisation constant, where we have suppressed the scale dependence of the strong coupling constant. The partonic cross section for the production of a Higgs boson through gluon fusion can then be cast in the form where we kept the squared Wilson coefficient factorised and pulled out all dimensionful prefactors, such that the ǫ 0 -piece of the leading order cross section becomes just All convolutions calculated in this work are done using theσ (n,m) ij and from here on, the term "cross section" will refer to these objects.
The sole dependence of the LO cross section on ǫ is an overall factor of (1 − ǫ) −1 = ∞ n=0 ǫ n from averaging over the D-dimensional polarisations of the initial gluons. Thus, the LO partonic cross section through O(ǫ 3 ) is trivially found to bẽ for all m = 0, . . . , 3.
At NLO, the dependence on ǫ is still fairly simple. There are only two master integrals and they are easily computed to all orders in ǫ.
The NNLO cross section through O(ǫ) necessitated the knowledge of the 29 master integrals to sufficiently high order in ǫ. The double-virtual master integrals can be found in work on the two-loop gluon form factor [48][49][50]. The real-virtual and double-real master integrals were computed by two groups independently during the last year [39,40]. The expression for the bare NNLO cross section in terms of master integrals was kindly provided to us by an author of [11].
In general, the partonic cross sections consist of three types of terms, delta-, plus-and regular terms.σ is defined via its action on a test function f (x) with a finite value at x = 1, The full expressions for the partonic cross sections through NNLO can be found in the ancillary files accompanying this arXiv publication. They agree with the ones given in [41] after compensating for the factor of 1/z that was not included in that publication.

Splitting kernels
The splitting kernel describes the probability of a parton j emitting a collinear parton i carrying a fraction x of the momentum of the initial parton. The splitting kernels are known up to three loops (P ij ) and may all be found in [51,52]. Note some different conventions that we use, though. Since we chose to expand all our results in a = αs π as opposed to αs 4π as in [51,52], our kernels P Also, since by P (n) qg we mean the emission of a single quark of a given flavour, we differ from the expression in [51,52] which parametrises the emission of any quark, by a factor of 1 2N F . Furthermore, there is also a conventional difference to the splitting kernels used in [41]. The authors of that publication use the quark-quark splitting kernel as defined in eq. (2.4) of [51]. This kernel, which we shall denote byP qq is used in the DGLAP evolution of pdfs. To compute all contributions to the N 3 LO gluon fusion cross section, we have to distinguish different initial-state channels such as qq (quark-antiquark), qq (identical quarks) and qQ (quarks of different flavour) which are convoluted with different combinations of pdfs. Thus, for channel-by-channel collinear factorisation, we require the three distinct, flavourdependent quark-quark kernels P qq , P qq , P qQ , (3.11) which describe the emission of an identical quark, the emission of an antiquark of the same flavour, and the emission of a quark or antiquark of a different flavour, respectively. The latter two kernels vanish at the one-loop order, P qQ but are nonzero for higher orders. In the notation of [52], this corresponds to the kernels P q i q j and P q iqj . The relation betweenP qq and our kernels is given bỹ We are not aware of results in [41] that involve the flavour-dependent quark-quark kernels. We close this section by giving the expressions for the four LO splitting kernels. For the lengthy higher-order kernels, we again refer to the machine-readable files accompanying this publication. The two NNLO kernels were taken from [53] in Form format and then translated to Maple input. Their regular parts were tested against the Fortran routine in [53], and their δ(1 − z) and D n (1 − z) parts were checked against [51]. (3.14)

Computation of the convolutions
In this section we will describe the method we used to compute the convolutions of splitting kernels and partonic cross sections that are needed to cancel collinear divergences at N 3 LO. Let us remark that our method is different from the technique used in [41], where the convolutions were calculated in Mellin space (where convolutions turn into ordinary multiplications) and the problem was essentially the calculation of the inverse Mellin-transform.
In the following, we restrict ourselves to a single convolution. Since the convolution product is associative, any multiple convolution appearing in the N 3 LO cross section can be obtained by repeating the steps using the result of the first convolution and the next convolutant. As already mentioned, both the splitting kernels and the partonic cross sections consist of three types of terms, delta-, plus-and regular terms.
where the regular pieces c where Li 1 (z) = − log(1 − z). HPLs can be defined recursively via the integral H(a 1 , a 2 , . . . , a n ; z) = z 0 dtf a 1 (t) H(a 2 , . . . , a n ; t) , a i ∈ {−1, 0, 1} , and in the special case where all a i = 0, the HPL is defined as H( 0 n ; z) = 1 n! log n (z) . For more comprehensive information about harmonic polylogarithms, we refer to [54][55][56][57][58][59][60]. Any convolution involving a delta function trivially returns the other convolutant (whether it be another delta function, a plus distribution or a regular function), (4.7) Convolutions involving two plus-distributions are more involved, yet no integral actually has to be solved. We comment on their calculation and list results for the required plus-plus convolutions in appendix A. For the remaining two types of convolutions, we end up with an actual integral that we need to compute, Specifically, a MPL may be a function of multiple variables that appear anywhere in the index vector (x 1 , . . . , x n ). The relation to HPLs reads H(a 1 , . . . , a n ; x) = (−1) k G(a 1 , . . . , a n ; x) , where k is the number of +1 indices in (a 1 , . . . , a n ). This sign difference is due to the fact that HPLs historically use 1 1−t as the weight function when adding a +1 to the index vector.
For more detailed information on multiple polylogarithms, see references [61][62][63] and references therein. Note that the order of the MPLs indices is often reversed. We follow the convention of [63]. The subsequent steps to solve the integrals are as follows: 1. We first remap the integral by x → 1 − x, such that the integration region becomes (0, 1 − z).

2.
HPLs with argument 1 − x and z 1−x have to be written as a combination of MPLs with the integration variable x as their argument, or no x-dependence at all. For example, where we've used that G(a; b) = log 1 − b a for a = 0. For MPLs of higher weights, one can find these translations by using the recursive definition of MPLs and changing variables in the integration. This becomes very tedious, though, so it proved to be more practical to use the symbol formalism developed in recent years [59,[63][64][65] and to follow the method presented in appendix D of [42]. For technical details, we refer the reader to said appendix.

At this stage, all integrations have been performed. The result still contains MPLs
where the variable z appears multiple times in the argument vector. Using the techniques from appendix D of [42] again, we can rewrite all the expressions in HPLs. 5. The final numerical check on the result consists of the comparison of the original integral using Mathematica's numerical integration (using the package HPL [57,58] to evaluate the HPLs numerically) and our final expression, using Ginsh, the interactive frontend of the computer algebra system GiNaC [56], for a random value of z.
The full set of convolutions can be found in machine-readable form (both Maple and Mathematica) in the ancillary files accompanying this arXiv publication. They were all compared analytically in Mathematica to the expressions given in [66], and complete agreement was found for all convolutions. For convolutions involving the two-loop quarkquark splitting kernels P qq and P (1) qQ , the results had to be combined according to eq. (3.12) to find equality.

Soft expansion of the convolutions
While the full N 3 LO corrections to the gluon fusion cross section may still be out of reach for the time being, a description in the soft limit could be feasible already in the close future. Note that this was the sucession at NNLO, as well, where the expansion of the cross section up to O((1−z) 16 ) [10,67] was published before the full computation [11,12]. The numerical agreement between the two computations proved to be excellent, so, anticipating the same behaviour at N 3 LO, the soft expansion of the N 3 LO corrections would be a very important result to obtain. The first pure N 3 LO piece of the third order soft expansion, the soft triple-real emission contribution, has recently been published [42].
In the limit z → 1, the partonic cross section (and all convolutions contributing to it) can be cast in the following form (suppressing partonic indices) We thus need to expand the regular part as a polynomial in (1 − z), times log(1 − z) terms. We proceed as follows: 1. We define z ′ ≡ 1 − z. Thus, our expressions now consist of HPLs with argument 1 − z ′ times powers of z ′ . The desired limit is z ′ → 0.
2. We want to rewrite the HPLs with argument 1 − z ′ as MPLs with argument z ′ , which results in changing the array of indices from {−1, 0, 1} to {0, 1, 2}, as can be easily seen by taking the integral definition eq. (4.10) for x 1 ∈ {−1, 0, 1} and changing variables t → 1 − t. The rewriting is achieved once again with the techniques from appendix D of [42].
3. The expansion of any MPL in its argument is straightforward, since there is a connection between MPLs and multiple nested sums [61], where the translation from MPLs to nested sums is given by The specific form of the MPL on the left-hand side of the equation above can be obtained via the scaling property, G(x 1 , . . . , x n ; z) = G(λx 1 , . . . , λx n , λz), where λ = 0 = x n . MPLs with a rightmost index of 0 must be rewritten using the shuffle product, e.g.
until all rightmost zeroes have been turned into explicit logarithms. The remaining MPLs can then be safely translated to nested sums.
The crucial point is that the variable x k only appears in the outermost sum in eq. (4.14), while the inner nested sums only depend on the x i<k , which in our case are the indices a i ∈ {0, 1, 2}. We thus easily obtain the desired expansion when we just truncate the sum over n k in eq. (4.14) at the highest power of z ′ we are interested in.
4. The validity of the soft expansions of the convolutions was checked numerically for some small values of z ′ . All soft expansions of the convolutions up to O (1 − z) 12 can be found in the ancillary files accompanying this arXiv publication.

Numerical results for the gluon fusion scale variation at N 3 LO
The total cross section for Higgs production through gluon fusion at N 3 LO depends on the factorisation and renormalisation scales explicitly, through logarithmic terms that have been derived in this work, and implicitly through the µ r dependence of α s and the µ f dependence of the parton distribution functions. In principle N 3 LO parton distribution functions should be used, but in practice, not only are they not available (nor will they be in the near future), but also their deviation with respect to the NNLO pdfs available is expected to be very small. On the other hand, the full, implicit, µ r dependence through α s , can only be estimated once the N 3 LO matrix elements are known, and in particular the coefficients 2 a are known, from mass factorisation constraints [37]. The question, then, arising is whether we can anticipate the scale uncertainty at N 3 LO with the information currently available.
To this end we parametrize the unknown delta and regular coefficients by a scaling factor K times the corresponding NNLO coefficients: There is no a priori reason why the scaling factor for the delta and the regular terms should be the same. However, it turns out that the numerical impact of the delta coefficient a is negligible (for scaling coefficients that do not break by orders of magnitude the pattern observed from lower orders), in contrast with the coefficient of the regular part, so we adopt here a common scaling factor to keep the parametrisation simple. For the same reason we use the same K scaling coefficient for all initial state channels. A loose argument about the size of K can be derived if one assumes a good perturbative behaviour at µ r = µ f = m H where all other terms of order a 5 vanish. Since a(m H ) ∼ 1/30 one expects K not to be much larger than 30. For comparison, the corresponding rescaling factors between NNLO and NLO are for m h = 125 GeV and µ f = µ r = m h . In what follows we study the inclusive cross section as a function of the scales, in the HQET approach, rescaled with the exact leading order cross section. We use the framework of the iHixs program [20] a Fortran code which contains the complete NNLO cross section for gluon fusion in HQET. The coupling α s was run to four-loop order according to eq. (2.13), while for the parton distributions, the MSTW08 NNLO set was used. Furthermore, to cross-check our results, a second implementation was programmed in C++, where the convolutions of splitting kernels and partonic cross sections were performed numerically. For both codes, the numerical evaluation of HPLs was performed using the library Chaplin [60]. The two implementations agreed for all parameter configurations that were tested.
In figure 1, the different orders of the hadronic gluon fusion cross section for the 8 TeV LHC and a Higgs mass of 125 GeV, along with several N 3 LO approximants for various numerical values of K are plotted as a function of the renormalisation scale µ r , while the factorisation scale is fixed to µ f = m h . Note that the convolutions of splitting kernels and partonic cross sections do not enter in this plot, since they are proportional to log(µ 2 f /m 2 h ). The µ r scale variation for LHC with 14 TeV centre-of-mass energy is shown in fig. 3. The µ f scale dependence, shown in figure 5 for 8 TeV centre-of-mass energy, is, as expected, extremely mild, in accordance with what is observed at NNLO.  Figures 2 and 4 display the overall scale dependence, with both scales set to be equal and varied simultaneously. We note that the curves for the approximate N 3 LO cross section with various Ks spread widely in the low scale region, i.e. for µ < 30 GeV. This is not unreasonable, though, as in this regime, the unknown N 3 LO contributions that are neglected become much more important due to the running of α s . Indeed, at the lowest renormalisation scale considered, µ = m h /16 ≈ 7 GeV, the coupling becomes as big as which is supposed to cancel the implicit logarithms in the running of α s , and which becomes large and negative, thus pulls the curve down for small scales, and is canceled by the currently unknown contributions whose magnitude is small at µ r = m h but is greatly enhanced due to α s at small µ r . It can hardly be overemphasised that the above prescription does not represent a proper calculation of the N 3 LO matrix elements, but just a way of parametrising their unknown numerical importance. Once the height of the N 3 LO curve at (µ r , µ f ) = (1, 1) is set, the shape of the full curve only depends on lower order cross sections (which we know exactly), the running of α s and the parton distribution functions, respectively. As mentioned above, the unknown, numerically important coefficient functions c (3,0) gg (z) contain logarithmic contributions that are singular at threshold, log(1 − z), contributions that are regular and contributions that are singular at the opposite, high energy limit log(z). The leading and several, but not all, subleading threshold contributions are associated with multiple soft emissions and can be recovered by resummation techniques. The Indeed, by comparing our results for the µ r -dependence of the N 3 LO cross section for the dominant gluon gluon initial state, with the numbers obtained via the recently released numerical program gghiggs [38], we find agreement between the two curves when setting K to 25, as is displayed in figure 6.
While it is plausible that the leading logarithmic contributions, being threshold enhanced, capture the bulk of the cross section, it is unclear whether the unknown subleading contributions, as well as the non-logarithmic terms, are really negligible. Their importance certainly rises for the LHC at 14 TeV, as the luminosity function suppresses the region away from threshold less, resulting in more phase space for real radiation. One might, therefore, want to be conservative about their magnitude, and hence on the size of the scale uncertainty to be anticipated before the full N 3 LO result is available. Table 1 shows the estimates for various values of the rescaling factor K, covering the range from relatively mild to extremely strong N 3 LO corrections, resulting in scale uncertainties varying from 2% to as large as 8% or more. The scale uncertainties cited here are evaluated by varying the scales in the interval [m h /4, m h ]. µ f is fixed to m h and only µ r is varied. K is varied from 0 to 30. Only the gg channel is plotted, and compared to the results obtained with [38].
The choice of the central scale around which the variation is performed has been an issue of debate lately, since different choices result in slightly different scale uncertainty estimates but also in different central values for the cross section. The choice is largely arbitrary, but various indications (like improved perturbative convergence, typical transverse momentum scales for radiated gluons, average Higgs transverse momentum etc.) point to a central scale choice that is lower than the traditional one at m h , closer to m h /2. An alternative indication comes from the considerations of [68], where it is argued, looking at examples from jet physics, that a reasonable indication would be the position of the saddle point in a contour plot of the cross section as a function of µ r and µ f . In figs. 7 and 8 we show such contour plots for Higgs production at LO, NLO, NNLO and N 3 LO (for three values of the parameter K). In the cases where a saddle point exists, its position points indeed to lower scale choices, and in the cases without a saddle point the plateau region is also located in lower scales. Given the extremely mild factorisation scale dependence, the saddle point or plateau region is largely determined by the µ r plateau in all previous figures.  The value on the contours is the cross section in picobarns. The x-axis is log 2 (µ f /m h ), the y-axis log 2 (µ r /m h ). Our preferred central scale choice is located at (-1,-1).

Conclusions
In this work we have presented all convolutions of lower-order partonic cross sections and splitting kernels that contribute to order a 5 to Higgs production in gluon fusion. The results agree with the ones previously published in [41]. Apart from the full expressions, we also provide all convolutions expanded around threshold, as the full N 3 LO corrections in this limit seem to be feasible in the near future.
We have also anticipated the scale dependence of the N 3 LO gluon fusion cross section, The value on the contours is the cross section in picobarns. The x-axis is log 2 (µ f /m h ), the y-axis log 2 (µ r /m h ). Our preferred central scale choice is located at (-1,-1).
into which the calculated convolutions enter. As is the case at NNLO, the factorisation scale dependence is extremely mild, at the per mille level or below. The overall scale uncertainty is driven by the renormalisation scale variation. The definite uncertainties depend on the size of the missing pure N 3 LO contributions. Scanning over a reasonable range for these contributions, we find that in the residual scale uncertainty can vary from 2% − 8% depending on the magnitude of the hard real corrections, whose computation is, to our view, a prerequisite for a solid estimate of the N 3 LO scale uncertainty. to thank Franz Herzog for pointing out the way of obtaining plus-plus convolutions by expanding a hypergeometric function. This work was supported by the Swiss National Foundation under contract SNF 200020-126632.

A. Convolutions of two plus-distributions
In the convolutions needed for collinear counterterms we face the problem of convolutions involving one or more plus-distributions. Here, we demonstrate how to obtain all convolutions containing two plus-distributions. In the second to last step, we mapped λ → 1 − λ and in the last step, the Euler definition of the hypergeometric function was used, where B(x, y) denotes the Euler Beta-function. On the other hand, we may also directly expand the integrands in I ab in terms of a delta function and a tower of plus-distributions, The above expressions agree with the ones given in [41] (eq. 22) and [70] (eq. C.28 -C.31). For the cases D 0 ⊗ D n , the combination of harmonic polylogarithms given in the references collapses to the single term − log(z) log n (1 − z)/(1 − z).