Charged Current Drell-Yan Production at N3LO

We present the production cross section for a lepton-neutrino pair at the Large Hadron Collider computed at next-to-next-to-next-to leading order (N3LO) in QCD perturbation theory. We compute the partonic coefficient functions of a virtual $W^{\pm}$ boson at this order. We then use these analytic functions to study the progression of the perturbative series in different observables. In particular, we investigate the impact of the newly obtained corrections on the inclusive production cross section of $W^{\pm}$ bosons, as well as on the ratios of the production cross sections for $W^+$, $W^-$ and/or a virtual photon. Finally, we present N3LO predictions for the charge asymmetry at the LHC.


Introduction
Cross sections for the production of leptons are among the ultimate precision observables measurable at the Large Hadron Collider (LHC) (see for example refs. [1][2][3]). As a consequence, they provide a unique window into the inner workings of collision processes at very high energies. The insights gained by studying them improve our understanding of the collider experiment and the fundamental mechanisms of scattering processes alike. To derive meaningful conclusions from such observations we must put our theoretical prediction for such scattering processes at the highest possible level. Here, we take one significant step in this direction and compute next-to-next-to-next-to leading order (N 3 LO) predictions for the production cross section of a lepton-neutrino pair in QCD perturbation theory. We focus on the charged current Drell-Yan [4] (CCDY) process, where a W boson is produced from the annihilation of two quarks with opposite isospin. In nature, the produced W boson decays within the blink of an eye, but the inclusive production probability of the bosons is easily related to the probability of producing a pair of fermions consisting of a neutrino and a charged lepton with invariant mass Q 2 . These particles represent a stable final state and the charged lepton can be detected by the LHC experiments.
To derive theoretical predictions for the inclusive CCDY cross section we use the factorisation of hadronic cross sections into parton distribution functions (PDFs) and partonic cross sections. In order to achieve high precision predictions of the hadronic cross section, it is paramount to go beyond the Born approximation for the partonic cross section. In particular, we compute the desired partonic cross sections in the framework of perturbative QCD through N 3 LO in the perturbative expansion in the strong coupling constant. The analytic formulae for the partonic cross sections are some of the main results of this paper and are provided in electronic form together with its arXiv submission.
We then move on and study the impact of N 3 LO corrections on various inclusive cross sections involving vector bosons. We focus on the progression of the perturbative series in QCD for some of the cleanest hadron collider observables. When combined with the available N 3 LO results in the literature [5][6][7][8][9][10][11][12][13][14], our results are an important input to gauge the relevance and the impact of N 3 LO corrections on more differential observables, like fiducial cross sections. We also study cross section ratios for vector boson production, and observe a remarkable perturbative stability for these ratios.
This article is structured as follows: In Section 2 we briefly review the computation of the N 3 LO corrections to off-shell W production. In Section 3 we present our main result, namely phenomenological predictions for the W cross section at N 3 LO in QCD, and we discuss the main QCD uncertainties which affect the cross section at this order. In Section 4 we extend this analysis to ratios of vector boson cross sections and the charge asymmetry at the LHC. In Section 5 we draw our conclusions.

Setup of the computation
In this paper we compute higher-order corrections in the strong coupling constant to the charged-current Drell-Yan (CCDY) cross section, i.e., the inclusive cross section for the production of a lepton-neutrino pair of invariant mass Q 2 at a proton-collider with centerof-mass energy √ S. Since we are only interested in QCD corrections, the lepton-neutrino pair can only be produced via the intermediate of an (off-shell) W boson. We can then factorise the production of the W boson from its subsequent decay and cast the cross section in the following form: where v is the vacuum expectation value and m W and Γ W are the mass and width of the W boson, and σ W ± denotes the inclusive cross section for the production of an off-shell W ± boson with virtuality Q 2 . Using QCD factorisation, this cross section can be written in the form where µ F and µ R denote the factorisation and renormalisation scales, and the f i (x, µ 2 F ) denote the parton density functions (PDFs) to find a parton species i with momentum fraction x inside the proton. Furthermore, τ = Q 2 S and η ± ij is the partonic coefficient for the production of an off-shell W boson from the parton species i and j. In the above equation we made use of convolutions defined by and we introduced the normalisation factor where n c corresponds to the number of colours in QCD.
The partonic coefficients can be expanded into a perturbative series in the renormalised strong coupling constant a S = α S (µ R )/π Above we have suppressed arguments of the functions indicating the dependence of the partonic coefficients on the perturbative scales. At leading order (LO) in α S , it is only possible to produce a W boson from the annihilation of two (massless) quarks with opposite isospin: η where V u i d j denotes the Cabibbo-Kobayashi-Maskawa quark-mixing-matrix. Beyond LO [4], also other partonic channels open up. Perturbative corrections to the CCDY cross sections have been computed at next-to-leading order (NLO) in refs. [15][16][17][18] and at next-to-nextto-leading order (NNLO) in refs. [19][20][21][22]. The main result of this paper is to present for the first time phenomenological results for the production of an off-shell W boson at next-tonext-to-next-to-leading oder (N 3 LO) in perturbative QCD. Before we discuss our results in the next sections, we review in this section the main steps of the computation of the third-order corrections to the partonic coefficients. The partonic coefficients are computable from Feynman diagrams, with z = Q 2 /ŝ and s the partonic center-of-mass energy. We have followed the same strategy as that for the computation of the inclusive cross section for Higgs boson production through gluon fusion and bottom-quark fusion [6,8,9,23] and the inclusive Drell-Yan cross section [10]. In particular, the results were obtained using the framework of reverse unitarity [24][25][26][27][28] in order to compute all required interferences of real and virtual amplitudes contributing to the N 3 LO cross section. The required phase-space and loop integrals were carried out implicitly using integration-by-parts (IBP) identities [29][30][31], together with the method of differential equations [32][33][34][35][36]. This method allows one to represent the required integrated and interfered amplitudes in terms of linear combinations of master integrals. The purely virtual corrections are essentially identical to the case of the production of an off-shell photon, apart from the colour structure involving a cubic Casimir operator, which is absent here because of the non-diagonal flavour-structure of the charged-current interactions. The three-loop corrections for virtual photon production were first computed in refs. [37][38][39][40][41][42][43], and recomputed and confirmed in ref. [8]. Contributions with one real parton in the final state were considered in refs. [44][45][46][47][48][49] and the master integrals we used for our calculation were documented in refs. [44,48]. Master integrals with two and three real partons were obtained for the purpose of ref. [6] and are based on results from refs. [23,[50][51][52][53][54].
We work in the MS-scheme in conventional dimensional regularisation. The only interaction vertex that involves an axial coupling is the vertex involving the W boson. The non-diagonal flavour-structure of the charged-current interactions forces the W boson to be coupled twice to same connected fermion line in interference diagrams. As a consequence, vector and axial vector contributions to the hadronic cross section are identical and we only work with a vector current in the generation of our partonic coefficient functions. The ultraviolet (UV) counterterm for the strong coupling constant has been determined through five loops in refs. [55][56][57][58][59]. Infrared (IR) divergences are absorbed into the definition of the PDFs using mass factorisation at N 3 LO [60][61][62]. The mass factorisation involves convoluting lower-order partonic cross sections with the three-loop splitting functions of refs. [63][64][65]. We have computed all the convolutions analytically in z-space using the PolyLogTools package [66]. After combining our interfered matrix elements with the UV and PDF-IR counterterms we send the dimensional regulator to zero and obtain our final results. Our partonic coefficients are expressed in terms of the iterated integrals in z defined in ref. [6]. While many of these iterated integrals can be expressed in terms of harmonic or multiple polylogarithms [67,68], some integration kernels involve elliptic integrals. It is currently unknown how to express them in terms of known functions. In ref. [6] it was shown how to obtain fast converging series representations which allow one to achieve a relative numerical precision of (at least) 10 −10 in the whole range z ∈ [0, 1] for all partonic coefficients.
The analytic results for the partonic coefficients are provided as ancillary material with the arXiv submission. Besides the explicit analytic cancellation of the UV and IR poles, we have performed various checks to establish the correctness of our computation. First, we have checked that all logarithmic terms in the renormalisation and factorisation scales produced from the cancellation of the UV and IR poles satisfy the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [69][70][71]. Second, the limit of soft-gluon emission, which corresponds to z → 1, is independent of the form of the hard interaction process and only depends on the quantum numbers of the initial-state partons. Hence, the soft-virtual cross sections must be identical to the ones for neutralcurrent Drell-Yan production (again, apart from the contribution from the cubic Casimir operator). We have reproduced the soft-virtual N 3 LO cross section of refs. [51,[72][73][74][75], and also the physical kernel constraints of refs. [76][77][78] for the next-to-soft term of the quark-initiated cross section. The high-energy limit of the cross section, which corresponds to z → 0, must also be identical between charged-and neutral-current production, and we have checked that our partonic cross sections have the expected behaviour in the highenergy limit [79,80]. Finally, we have checked that we reproduce the numerical results for the ratios of the inclusive W + and W − cross sections up to NNLO given in ref. [81].

W -production at N LO
In this section we discuss our phenomenological results for inclusive (off-shell) W -production at N 3 LO. All results were obtained in a theory with N f = 5 massless quark flavours (two up-type and three down-type flavours). The values of the mass and the width of the W boson are m W = 80.379 GeV and Γ W = 2.085 GeV, and the vacuum expectation value is v = 246.221 GeV. The elements of the CKM matrix relevant to our computation are [82]: Note that only the absolute values of the entries of the CKM matrix enter our computation. In our phenomenological results the top quark is absent, which is equivalent to having V td i = 0 and only considering N f = 5 massless degrees of freedom in loops. This approximation is motivated because off-diagonal CKM matrix elements are small and diagrams without a coupling of the top quark to the electroweak gauge boson decouple in the limit of infinite top quark mass. Corrections to this approximation, which are expected to be very small, can be computed separately and are beyond the scope of this article. The strong coupling constant is evolved to the renormalisation scale µ R using the four-loop QCD beta function in the MS-scheme assuming N f = 5 active, massless quark flavours. Unless stated otherwise, all results are obtained for a proton-proton collider with √ S = 13 TeV using the zeroth member of the combined PDF4LHC15 nnlo mc set [83].    2 show the dependence of the fixed-order cross sections on the factorisation scale µ F and renormalisation scale µ R , which are introduced by the truncation of the perturbative series. We show the variation of the cross section for Q = 100 GeV on one of the two scales with the other held fixed at Q. We observe that the dependence on the perturbative scales is substantially reduced as we increase the perturbative order. The dependence on the scales looks very similar to the case of the N 3 LO cross section for the neutral-current process studied in ref. [10].
Next, in order to quantify the size of the N 3 LO corrections, we investigate the Kfactors: is the hadronic cross section including perturbative corrections up to n th order evaluated for µ F = µ R = Q and δ scale (σ (n) W ± ) is the absolute uncertainty on the cross section from varying µ F and µ R independently by a factor of two up and down around the central scale µ cent = Q such that 1 2 ≤ µ F µ R ≤ 2 (7-point variation). Table 1: The QCD K-factor at N 3 LO for charged-current and neutral-current Drell-Yan production. All central values and uncertainties are computed according to eq. (3.2). The results for neutral-current Drell-Yan production are taken from ref. [10]. Table 1 shows the results for the K-factors for both charged-and neutral-current Drell-Yan production for various values of Q (the results for the neutral-current process are taken from ref. [10]). We observe that the K-factor reaches up to 5% for small values of Q, and the impact of the N 3 LO corrections gets smaller as Q increases. Moreover, we find that the K-factors at N 3 LO are identical between the charged-and neutral-current processes, confirming and extending results at lower orders, and hinting towards a universality of QCD corrections to vector-boson production in proton collisions. Figure 3 shows the values of the cross section normalised to the N 3 LO cross section as a function of the virtuality Q. The uncertainty bands are obtained by varying the renormalisation and factorisation scales independently up and down as described above around the central scale µ cent = Q. We observe that for Q 50 GeV the scale-variation bands at NNLO and N 3 LO do not overlap. A similar feature was already observed for  virtual photon production in ref. [10], hinting once more towards a universality of the QCD corrections to these processes.   The fact that the scale variation bands do not overlap puts some doubt on whether it gives a reliable estimate of the missing higher orders in perturbation theory, or whether other approaches should be explored (cf., e.g., refs. [84,85]). In ref. [10] it was noted that for virtual photon production there is a particularly large cancellation between different initial state configurations. We observe here the same in the case of W boson production. This cancellation may contribute to the particularly small NNLO corrections and scale variation bands, and it may be a consequence of the somewhat arbitrary split of the content of the proton into quarks and gluons. If these cancellations play a role in the observed perturbative convergence pattern, then it implies that one cannot decouple the study of the perturbative convergence from the structure of the proton encoded in the PDFs. We will return to this point below, when we discuss the effect of PDFs on our cross section predictions.    ����� ������� � -��� ����� ���������_����_�� μ�����=� Figure 6: The cross sections for producing a lepton-neutrino pair via an off-shell W boson as a function of the invariant mass of the final state, or equivalently the virtuality of the W boson, cf. eq. (2.1). Figure 6 shows the nominal production cross section of a lepton-neutrino pair at the LHC at 13 TeV centre of mass energy, as defined in eq. (2.1). Figure 7 shows the variation of K-factors as a function of the energy of the hadron collider for Q = 100 GeV. The orange, blue and red bands correspond to predictions with the perturbative cross section truncated at NLO, NNLO and N 3 LO, and the size of the band is obtained by performing a 7-point variation of (µ F , µ R ) around the central scale µ cent = Q. We observe that the NLO, NNLO and N 3 LO K-factors are relatively independent of the centre of mass energy.

���
Parton distribution functions are extracted from a large set of measurements and are consequently subject to an uncertainty related to the input as well as to the methodology used to extract the PDFs. Here, we follow the prescription of ref. [83] for the computation of PDF uncertainties δ(PDF) using the Monte Carlo method. Furthermore, also  the strong coupling constant is an input parameter for our computation. The PDF set PDF4LHC15 nnlo mc uses α S = 0.118 as a central value and two additional PDF sets are available that allow for the correlated variation of the strong coupling constant in the partonic cross section and the PDF sets to α up S = 0.1195 and α up S = 0.1165. This sets allow us to deduce an uncertainty δ(α S ) on our cross section following the prescription of ref. [83]. We combine the PDF and strong coupling constant uncertainties in quadrature to give In our computation we use NNLO-PDFs, because currently there is no available PDF set extracted from data with N 3 LO accuracy. It is tantalising to speculate if the observed convergence pattern is related to the mismatch in perturbative order used for the PDFs and the partonic cross section. We estimate the potential impact of this mismatch on our cross section predictions using a prescription introduced in ref. [5] that studies the variation of the NNLO cross section as NNLO-or NLO-PDFs are used. This defines the PDF theory uncertainty δ(PDF-TH) = 1 2 Here, the factor 1 2 is introduced as it is expected that this effect becomes smaller at N 3 LO compared to NNLO. Figure 8 displays the uncertainties δ(PDF), δ(PDF+α S ) and δ(PDF-TH) as a function of Q in orange, red and green respectively. In particular, the green band indicates the sum δ(PDF + α S ) + δ(PDF-TH). We observe that the estimate for δ(PDF-TH) plays a significant role especially for low values of Q. The traditional PDF uncertainty has a stronger impact for larger values of Q. Overall, we observe that the relative size of δ(PDF) and δ(PDF-TH) is large in comparison to the effect of varying the scales. We conclude that future improvements in the precision of the prediction of this observable will have to tackle the problem of the uncertainties discussed here. In particular, we emphasize

Predictions for cross section ratios
In the previous section, we have seen that the conventional variation of the perturbative scales by a factor of 2 does not give a reliable estimate of the size of the missing higher orders. This motivates us to study the ratios of cross sections for the production of gauge bosons with virtuality Q: Indeed, since the charged-and neutral-current Drell-Yan processes show very similar Kfactors and dependences on the perturbative scales, it is conceivable that some of the uncertainties (e.g., PDF effects) cancel in the ratio, and the ratios may exhibit an enhanced perturbative stability. In the remainder of this section we analyse and compare different prescriptions to estimate the missing higher orders in the perturbative expansion of the cross section ratios for √ S = 13 TeV. We focus here on the following prescriptions: • Prescription A: We take the ratio of the perturbative expansion of the numerator and the denominator computed at a given order in perturbation theory. We choose the renormalisation and factorisation scales in the numerator and denominator in a correlated way, i.e., we always choose the same values for the scales in both the numerator and the denominator.
• Prescription B: Similar to Prescription A, but we do not correlate the renormalisation and factorisation scales between the numerator and the denominator. In other words, we perform independently a 7-point variation of the scales in the numerator and the denominator, and we take the envelope of the values obtained.
• Prescription A : We choose the scales in a correlated way, but we expand the ratio in perturbation theory and only retain terms through a given order in the strong coupling.
• Prescription B : Similar to Prescription A , but we choose the scales in an uncorrelated way.
• Prescription C: We take the relative size of the last considered order compared to the previous one as an estimator of the perturbative uncertainty: The superscript n on R (n) XY (Q) indicates the order at which we truncate the perturbative expansion. The values of µ F = µ R = µ cent are chosen in a correlated way in the numerator and the denominator of R XY . This estimator is based on the expectation that in a well-behaved perturbative expansion the subleading terms should be smaller than the last known correction. By construction, it leads to an estimate of the missing higher orders that is symmetric around the central value.
Note that Prescriptions A, A , B and B are fully equivalent to all orders in perturbation theory. The truncation of the perturbative series can, however, introduce differences in the results obtained from these four prescriptions, especially at low orders in perturbation theory. For example, the bands obtained by varying the scales will always be larger if the scales are varied in an uncorrelated way, because in that case the size of the band is obtained by taking the envelope of a strictly larger set of values.   Table 2 shows the prediction for R W + W − for Q = m W computed at the first few orders in perturbation theory. First, we see that the central value is extremely stable in perturbation theory, changing only at the permille level as we go from NNLO to N 3 LO. The central value is pretty much independent of the choice of the central scale µ cent and whether the ratio is expanded in perturbation theory or not (primed vs. unprimed prescriptions). While in general the different prescriptions lead to vastly different estimates of the missing higher orders, the predictions are similar between the primed vs. unprimed prescriptions, especially as we increase the perturbative order. This is to be expected: If the perturbative order is increased, the differences stemming from expanding or not the denominator of the ratio should decrease, which is indeed what we observe. We therefore only discuss the unprimed prescriptions from now on. Second, we observe that Prescription B leads to an estimate that is more than an order of magnitude larger than for Prescriptions A and C. In particular, this makes one wonder if correlated scales (Prescription A & C) tend to underestimate the size of the missing higher-order terms beyond N 3 LO. We believe that results obtained from uncorrelated scales (Prescription B) lead to estimates that are too conservative. Indeed, since the central value of the ratio only receives permillelevel corrections from NNLO to N 3 LO and exhibits extremely good perturbative stability, one expects higher-order corrections to be at the sub-permille level, which is indeed the size of the band obtained by varying the scales in a correlated way (Prescription A & C). It would be unreasonable to expect that the missing higher-order corrections shift the central value by 1% or even more, which is the size of the bands obtained from the uncorrelated prescription (Prescription B). A correlated prescription is also motivated by the fact that the neutral-and charged-current processes are expected to receive very similar QCD corrections, a fact which is corroborated by the results from the previous section. Finally, we observe that Prescription C leads to an estimate that is always slightly larger (by a factor ∼ 3 for Q = m W ) than the one obtained from Prescription A at N 3 LO. We have observed that the size of the higher-order terms estimated from Prescription A always encompasses the next order in perturbation theory. In order words, Prescription A is a reliable estimator of missing higher-orders for on-shell W production at the LHC.  Figure 9 shows the ratio R W + W − as a function of the virtuality Q. The bands were created using Prescription A (left), B (middle) and C (right). Just like for on-shell production, we observe that correlated and uncorrelated scales lead to vastly different estimates for the size of the bands. In particular, Prescription B gives an extremely conservative estimate of the bands over the whole range of Q considered. Unlike for on-shell production, the bands obtained from Prescription A do not overlap for Q > 200 GeV, which indicates that Prescription A does not correctly capture the size of higher-order corrections in this range of virtualities. Instead, Prescription C seems to give the most reliable estimator of the size of the residual perturbative corrections for ratios of cross sections over the whole range of virtualities considered.  Figure 10: The ratios R W + W − (left), R W + γ * (middle) and R W + γ * (right) as a function of the virtuablity Q. The uncertainty bands are obtained with Prescription C for the central scale µ cent = Q.

���
In fig. 10 we extend our analysis from R W + W − (left) to R W + γ * (middle) and R W + γ * (right). In all cases the bands are estimated using Prescription C. Similar to our discussion in the previous paragraph, we find that Prescription C delivers reliable estimates also for the latter two ratios. In all cases we find that the residual perturbative uncertainty is very small, making ratios of production cross sections of electroweak gauge bosons very stable under perturbative corrections, and therefore ideal precision observables.
While Prescription C seems to give reliable estimates for all cross section ratios considered over a wide range of kinematics, we have to point out that Prescription C has the obvious shortcoming that it gives a vanishing result whenever two consecutive perturbative orders give identical numerical predictions. In order to estimate perturbative uncertainties we therefore suggest that multiple different perspectives and estimators should be considered, in order to test the quality of the estimates obtained. For the future, it would be interesting to investigate other prescriptions to estimate the impact of missing perturbative orders, including prescriptions based on statistical methods [84,85]. A more detailed study of these effects, however, would go beyond the scope of this paper.
We finish by presenting results for an observable which is closely related to the ratio R W + W − . The lepton-charge asymmetry is defined as: Figure 11 shows our predictions for the lepton-charge asymmetry as a function of Q at different orders in perturbation theory. All uncertainty bands are obtained from Prescription C. Just like the cross section ratios studied earlier, we observe a good perturbative stability and a very small residual dependence on the perturbative scales at N 3 LO. In particular, for Q = m W , we find  Figure 11: The lepton-charge asymmetry A W as a function of the virtuality Q. The uncertainty bands are obtained with Prescription C.
So far we have only discussed the uncertainties on cross section ratios from the truncation of perturbative orders. We therefore conclude by commenting on uncertainties on ratios related to PDFs, similar to those considered in Section 3. Just like in the case of the perturbative uncertainty, a choice has to be made whether or not to treat PDF and strong coupling constant variation as correlated in numerator and denominator of the ratio. The fact that PDFs and α S are universal quantities suggests indeed a correlated treatment. In this case, ratios of cross sections could provide a remarkable tool to reduce some of the largest theoretical uncertainties afflicting the observables in question.

Conclusion
In this paper we have computed for the first time the N 3 LO corrections to the inclusive production cross section of a lepton-neutrino pair at a proton-proton collider in QCD perturbation theory. One of the main results of this article are analytic formulae for the partonic coefficient functions for this cross section, which we make available as ancillary material with the arXiv submission of this paper.
We have studied the phenomenological impact of our results by providing numerical results for (off-shell) W -production at the LHC. All our results are differential in the virtuality of the W boson, or equivalently in the invariant mass of the lepton-neutrino pair. We find that the N 3 LO corrections are at the level of a few percent and they stabilise the perturbative progression of the series. We have also studied ratios of cross sections for W + , W − and γ * , as well as the charge asymmetry at the LHC, i.e., the amount of produced W + bosons relative to the amount W − bosons. We find that these ratios feature a remarkable perturbative stability and that they can be predicted with very high precision. For the future, it would be interesting to study, in addition to the QCD corrections discussed here, how the inclusive production probability of weak bosons is impacted by electroweak and QED corrections, for example the electroweak and mixed QCD-electroweak corrections to Drell-Yan processes of refs. [86][87][88][89][90][91][92][93][94][95][96]. While such a study is beyond the scope of this article, we would like to stress their importance here.
Combined with the results for virtual photon production of ref. [10], our results have also allowed us to investigate the progression of the perturbative series through N 3 LO in QCD for one of the simplest classes of hadron collider observables, namely fully inclusive vector boson production cross sections. Understanding the perturbative convergence of this class of processes is an important proxy to understand the precision that can be reached more generally for LHC observables at this order in perturbation theory. Here we summarise our main observations and conclusions, and the possible implications for a precision physics program at the LHC. First, we observe in all cases that the N 3 LO corrections shift the central value of the predictions by a few percent. This implies that, very likely, also for more differential observables percent-level precision can only be achieved after the inclusion of N 3 LO corrections. It would be interesting to develop and extend techniques to perform differential calculations at N 3 LO (see, e.g., refs. [12][13][14][97][98][99] for first steps in this direction). Second, we observe that the uncertainties related to PDFs generically dominate over the residual perturbative uncertainty at N 3 LO. This motivates to push for determining PDFs at this order in perturbation theory, by extracting them from experimental measurements confronted to theory calculations performed at the same order, and to evolve them using the (yet unknown) DGLAP evolution kernels at N 4 LO. Finally, we observe that the conventional method to estimate the missing higher-order terms by varying the factorisation and renormalisation scales by a factor of two around a hard scale does not give reliable results at N 3 LO. This calls for an improved method to estimate the missing higher-order terms, e.g., by studying the progression of the perturbative series as done here, or by considering statistically-motivated techniques such as those of refs. [84,85]. However, we have also observed that the determination of the residual perturbative uncertainty may not be completely decoupled from the study of PDF effects: Indeed, we observe substantial cancellations between different partonic channels. It is tantalising to speculate if these cancellations are responsible for the non-overlapping scale variation bands at N 3 LO, and if they persist once a complete set of N 3 LO PDFs is available. A detailed study of these effects, however, goes beyond the scope of this paper and is left for future work.