Do we need N$^3$LO Parton Distributions?

We discuss the uncertainty on processes computed using next-to-next-to leading (NNLO) parton distributions (PDFs) due to the neglect of higher order perturbative corrections in the PDF determination, in the specific case of Higgs production in gluon fusion. By studying the behaviour of the perturbative series for this process, we show that this uncertainty is negligible in comparison to the theoretical uncertainty on the matrix element. We then take this as a case study for the use of the Cacciari-Houdeau method for the estimate of theoretical uncertainties, and show that the method provides an effective way of treating theoretical uncertainties on the matrtix element and the PDF on the same footing.

Gluon fusion, the dominant Higgs production channel at the LHC, has a slowly convergent expansion in perturbative QCD: the inclusive cross section is currently known up to next-to-next-to-leading order (NNLO) [1][2][3], and a recent approximate determination of the N 3 LO result has been presented [4], while rapid progress on the exact computation has been reported [5].
With N 3 LO results around the corner, it is natural to ask whether these will be of any use, given that fully consistent N 3 LO parton distributions (PDFs) are not likely to be available any time soon, essentially because the determination of N 3 LO anomalous dimensions would require a fourth-order computation, for instance of deep-inelastic structure functions, or Wilson coefficients. Clearly, this question is related to the more general issue of theoretical uncertainties on PDFs: current PDF uncertainties [6] only reproduce the uncertainty in the underlying data, and of the procedure used to propagate it onto PDFs, but not that related to missing higher-order corrections in the theory used for PDF determination. Henceforth in this paper we will call 'theoretical uncertainty' the uncertainty due to the fixed-order truncation of the perturbative expansion, sometimes [7] also called missing higher-order uncertainty, or MHOU.
Here we address this set of issues in the specific context of Higgs production in gluon fusion. We use the dependence on the perturbative order of the prediction for this process as either the PDF or the matrix element are taken at different orders as an estimate the theoretical uncertainty on either. We then address the more general issue of how one may estimate theoretical uncertainties on PDFs and matrix elements, specifically by using the approach of Cacciari and Houdeau [8].
We first compute the cross-section using the ggHiggs code [4,9], with default settings 1 .
Results are shown in Figure 1, where we show the cross section evaluated at increasingly high perturbative order (henceforth loosely referred to as the "order of the matrix element"), also including the approximate N 3 LO from Ref. [4], using in each case LO, NLO or NNLO PDFs (henceforth referred to as the "order of the PDF"). We use NNPDF2.3 PDFs [12] (with NNPDF2.1 LO [13] as LO set [14]). What is shown here is the total cross-section at the hadronic level, obtained summing over all parton subchannels, except at N 3 LO, where only the gluon-gluon channel is included in the estimate of Ref. [4]. We assume α S (M Z ) = 0.119 in all cases, as we are interested in studying the behaviour of the perturbative series for a fixed value of the coupling constant. The uncertainty bars in Figure 1 are all obtained by varying the renormalization scale in the range m H /2 ≤ µ R ≤ 2m H (the choice µ R = m H /2 as central scale is sometimes advocated instead [10], as it leads to faster convergence of the perturbative expansion: this choice would not change our conclusions). The variation of the renormalization scale should provide an estimate of the missing higher-order corrections to the matrix element when the PDF is kept fixed 2 , though, as well known, for this process it substantially underestimates them.
Be that as it may, it is clear that the dependence of the result on the order of the matrix element is much stronger than the dependence on the order of the PDF: on the scale of the variation of the matrix element, results obtained when the order of the PDF is varied are almost identical (especially beyond LO). Hence we could conclude here our brief investigation, having answered in the negative the question which is asked in the title: at least as far as Higgs in gluon fusion is concerned, based on the behaviour of the perturbative expansion at known orders, it is very likely that using NNLO PDFs in the N 3 LO computation would lead to results which are essentially indistinguishable from those consistently obtained using N 3 LO PDFs at N 3 LO.
However, it is worth elaborating a little more on our result. Specifically, it would be desirable to be able to provide a quantitive estimate of the theoretical uncertainty on the PDF, as well as of the combined theoretical uncertainty on the hadron-level process due to both the matrix element and the PDF. In principle, it is possible to use scale variation in order to determine the uncertainty on the PDF, too: it is, however, quite cumbersome in practice as it requires keeping track of the scale variation during the PDF fitting [15]. Indeed, to the best of our knowledge, it has never been done for any of the available PDF sets. Also, correlations between the behaviour of different processes upon scale variation would have to be kept into account: correlations between processes used for PDF determination among each other would be needed in order to determine the uncertainty on the PDF, and correlations between them and the processes for which a prediction is sought would be required in order to be able to combine uncertainties on the PDF and the matrix element.
On the other hand, it is clear by looking at Figure 1 that a good deal of information is contained in the behaviour of the perturbative expansion itself: it is then natural to try to systematically provide a determination of the uncertainty bar based on previous orders, rather than on scale variation. This has the further advantage that the theoretical versions of this code instead disagreed with ggHiggs, because of bugs affecting the top mass dependence at NLO and the factorization scale. Minor differences persist in the top mass dependence, but ggHiggs fully agrees with the non-public code of Ref. [11]. The cross section for Higgs production in gluon fusion computed with NNLO PDFs at increasing perturbative orders. At each order the uncertainty is shown as determined (from left to right) using scale variation (red circles, same as Fig. 1), the Cacciari-Houdeau method (blue crosses), and the same method but with rescaled parameter (see text, green squares); at N 3 LO the Passarino-David uncertainty is also shown (see text, purple diamonds).
uncertainty due to the matrix element and the PDF could then be treated on the same footing, and easily combined. A methodology to do so has been suggested Ref. [8] (Cacciari-Houdeau method, henceforth), based on assuming a prior distribution for the coefficients of the perturbative expansion, and then using Bayesian arguments to determine a confidence interval for the unknown coefficients based on the behaviour of the known ones. In Figure 2 we compare the theoretical uncertainty on the matrix element computed by scale variation, already shown in Fig. 1 (corresponding to points shown in the plot as red circles) to that obtained using the Cacciari-Houdeau method, i.e. essentially Eq. (85) of Ref. [8] (shown as blue crosses). In all cases, the PDF is kept fixed to the NNLO set. It should be noticed that this is a very simple-minded application of the ideas of Ref. [8]: in particular, we do not distinguish between different partonic subchannels, which could be in principle characterized by different perturbative expansion, and we study the perturbative behaviour of the total hadronic cross section rather than, for instance, the perturbative behaviour of the differential partonic cross-section.
The result of Fig. 2 is clearly not very satisfactory: at each order, the uncertainty bar, rather than being of the same order as the shift when going to the next order, is much smaller than it, and also rather smaller than the scale uncertainty, which already underestimated this shift. This could of course be due to the crude approximations we are making, as discussed above.
However, there is a more fundamental consideration. Namely, the approach of Ref. [8] is based on the assumption that the perturbative expansion coefficients whose behavior is being studied are roughly all of the same order -at least at low orders, well below the point where the perturbative series, which is at best asymptotic, starts diverging. The result shown in Fig. 2 has been obtained by writing the cross-section as a series and identifying the expansion coefficients σ i with the coefficients c i as given in Eq. (85) of Ref. [8]. However, it turns out that the coefficients σ n thus defined rapidly grow with the perturbative order. On the other hand, it is clear that (as already pointed out in Ref. [8]) another choice of expansion parameter would generally lead to a different behaviour. Indeed, the natural expansion parameter, even in the simplest cases, usually differs by the strong coupling by a (possibly large) factor: it might be given, for instance, by α S 4π or C A α S .
Lacking an analytic knowledge which may motivate a choice of the expansion parameter (especially in view of our very simple-minded approach), we rewrite Eq. (1) by rescaling α S by a real parameter λ: The approach of Ref. [8]) is then applicable if there exists a value of λ such that the rescaled coefficients c λ i are all of comparable order. We then simultaneusly test for the applicability of the method, and determine the optimal value of λ, by letting c λ n = κ and then performing a two-parameter fit of λ and κ to the three known coefficients (including the approximate N 3 LO result of Ref. [4]). We get an almost perfect fit (χ 2 below 1%), and a best-fit value of λ = 5.6, with NNLO PDFs and µ R = m H . The good quality of the fit means that the rescaled coefficients are indeed all of the same order, thus justifying the use of the method, but the large rescaling which is required explains the failure of the method before rescaling. Note that the rescaled expansion parameter is large, but still smaller than one, as one expects for a slowly convergent series.
We have checked that the best-fit λ is quite stable upon variations of the procedure. In particular it varies by a few percent if we change the order of the PDF, or if we decide to also include the leading-order coefficient in the fit (i.e. if we fit directly the coefficients σ n of Eq. (1) as σ n = κλ n ): this latter choice leads to a significantly worse χ 2 , but with essentially the same λ. If we change the renormalization scale to µ R = m H /2 the optimal λ decreases by about 20%, to λ = 4.3, while the fit quality deteriorates to χ 2 ∼ 1.1, still justifying the use of the method, given that the equality of the coefficients is only expected to be approximate. 3 The fact that the rescaling is somewhat smaller is an interesting feature of the method: it shows that the perturbative expansion converges somewhat faster with this choice of renormalization scale, as it is known to be in fact the case.
Armed with the knowledge of the necessary rescaling λ = 5.6, we recompute the Cacciari-Houdeau uncertainty using the rescaled parameterᾱ s . The result is also shown in Fig. 2 (green squares). It is clear that now the result provides a rather reasonable estimate of the theoretical uncertainty, which, up to NNLO, turns out to be of the same order as the observed perturbative shift at each order, and thus in particular it reflects the theoretical uncertainty better than scale variation. At N 3 LO, scale variation and Cacciari-Houdeau lead to similar answers. For comparison, in Fig. 2 we also show (purple diamonds) the uncertainty on the N 3 LO result estimated according to the method of Ref. [7]. In this reference, the theoretical uncertainty is determined by assuming that the perturbative series is an asymptotic series which is summed using various techniques (such as Borel summation). For a same-sign series the uncertainty band is taken to be at any given order as the interval between the known result up to that order, and the upper (more in general, the extreme) all-order asymptotic sum -so the lower edge of the band coincides with the N 3 LO central value, by construction. Interestingly, the size of the uncertainty band on the N 3 LO result found using the method of Ref. [7] is very close to that from Cacciari-Houdeau (which, at this order, is also similar to scale variation as already mentioned). We now finally turn to a determination of the theoretical uncertainty due to either the PDF, or the matrix element, or both. Results are shown in Figure 3, for the hadronic cross-section computed at each order using consistently PDFs at the corresponding order (LO matrix element with LO PDFs and so on). All uncertainties are now computed using the Cacciari-Houdeau method. In order to determine the uncertainty on the matrix element as the PDF is kept fixed, we have used the rescaled method as discussed above, with the PDF kept fixed either to its LO, NLO, or NNLO value: we then show at LO the uncertainty on the matrix element when the PDF is kept fixed at LO, at NLO the uncertainty when it is kept fixed at NLO and so on. In order to determine the uncertainty due to the PDF, no rescaling turns out to be necessary, so we use the method with α S taken as expansion parameter. Finally, the combined uncertainty is determined applying the rescaled Cacciari-Houdeau method to the series at the hadronic level in which the order of the PDF and the matrix element are varied simultaneously (with the NNLO PDF used also at N 3 LO); the same rescaling is used as for the uncertainty on the matrix element only (indeed, inclusion of the PDF changes the best-fit rescaling by an amount which is essentially irrelevant). Comparing Figure 3 with Figure 1 shows again that the Cacciari-Houdeau method, with rescaling when necessary, provides an estimate of the theoretical uncertainty which is in reasonable agreement with the behaviour of the perturbative expansion at the known orders. This supports its use in order to estimate theoretical uncertainties at the highest order at which they are known exactly (NNLO) or approximately (N 3 LO). Figure 3 confirms the conclusion we already reached by inspection of Figure 1, namely, that the dependence of results for Higgs production in gluon fusion on the perturbative order of the PDF is much weaker than that on the perturbative order of the matrix elementwhich, as well known, is unusually large. We conclude that an exact determination of the N 3 LO perturbative correction to the matrix element will lead to a substantial reduction of the theoretical uncertainty on the cross section for Higgs production in gluon fusion, even without knowledge of N 3 LO parton distributions.
While in the specific case of Higgs in gluon fusion the negligible impact of N 3 LO corrections to PDFs follows almost trivially from the huge hierarchy between the uncertainty on the PDF and that on the matrix element, one may ask whether this is true in general.
To see this, we have also considered the case of top pair production. The analogue of the plot of Fig. 1 for the total top pair production cross-section is shown in Fig. 4. Results are obtained using TOP++2.0 [16], including the recent full NNLO result of Ref. [17]. Here too we take α s (M Z ) = 0.119, and we use NNPDF2.3 PDFs (in the version with maximum number of flavors N f = 5, as this is what TOP++2.0 requires). Uncertainty bars are now obtained by varying both the renormalization and factorization scales by a factor two about the central value µ R = µ F = m t , with the ratio of the two scales constrained not to exceed two [18]. It is clear that, while the dependence on the order of the matrix element is still somewhat stronger than that on the order of the PDF, now the two are comparable, so that the dependence of the cross section on the perturbative order with fixed PDF differs by a non-negligible amount from that found when the order of the PDF is consistently varied along with that of the matrix element. In fact, because the dependence on the order of the matrix element and that on the order of the PDF are anti-correlated, the dependence of the physical (hadron-level) cross-section on the perturbative order is somewhat weaker than found if the order of the matrix element is varied while the PDF is kept fixed.
This example is sufficient to conclude that what is true for Higgs in gluon fusion is not true in general: for other processes N 3 LO corrections to PDFs might well be relevant. Also, the example raises two interesting questions. The first is the reason for this difference between Higgs and top. The question can be answered by studying the perturbative behaviour of the gluon luminosity, from which the dominant contribution to both processes originates, and comparing it to the correlation (defined as in Sect. 4 of [19]) between the gluon luminosity itself and the cross-sections which are being computed, see Fig. 5. It is clear that the correlation of the Higgs cross section to the gluon luminosity is rather strongly peaked in a region in which the gluon luminosity depends very weakly on the perturbative order, while the correlation of the top cross section is large in a significantly broader kinematical range (because even at LO the invariant mass of the final state is not fixed), including a region in which the perturbative dependence of the luminosity is sizable. This implies that, whereas the general behaviour of the cross-section may be easily understood in terms of the features of the relevant physical processes and of the parton luminosity (which in turn depends on the processes used for PDF determination), whether the perturbative dependence of the PDF is or not important has to be determined by a dedicated analysis of each process. In particular, this requires a systematic correlation analysis, such as that performed for several Higgs signal and background processes in Sect. 3.2 of Ref. [20]; along with a study of the perturbative dependence of parton luminosities.
The second question is how to best determine and use theoretical uncertainties on PDFs. The perturbative behaviour of the top cross section of Fig. 4 suggests that the theoretical uncertainty on the top cross-section would be overestimated if the uncertainty on the PDF were not included, i.e. that the latter actually reduces the uncertainty of the physical cross section. Hence, in this case, in order to properly include the theoretical PDF uncertainty one must keep into account its (anti)correlation with the theoretical uncertainty on the matrix element. It appears that this would be very difficult to do if theoretical PDF uncertainties were determined by scale variation. In this respect, a method such as Cacciari-Houdeau, based on the analysis of the perturbative behviour appears rather more promising.
A systematic investigation of both these issues will be left for further studies.