Lepton-pair production at hadron colliders at N$^3$LO in QCD

We compute for the first time the complete corrections at N$^3$LO in the strong coupling constant to the inclusive neutral-current Drell-Yan process including contributions from both photon and $Z$-boson exchange. Our main result is the computation of the QCD corrections to the inclusive production cross section of an axial-vector boson to third order in the strong coupling in a variant of QCD with five massless quark flavours. Since the axial anomaly does not cancel for an odd number of flavours, we also consistently include non-decoupling effects in the top-quark mass through three loops. We perform a phenomenological study of our results, and we present for the first time predictions for the inclusive Drell-Yan process at the LHC at this order in QCD perturbation theory.


Introduction
The production of a pair of massless leptons at hadron colliders like the Large Hadron Collider (LHC) at CERN -the so-called Neutral-Current Drell-Yan (NCDY) process -is one of the most important and most studied hadron collider processes.From an experimental perspective, its clean final-state signature makes it an ideal candidate for luminosity measurements and detector calibration.The DY process also plays a key role in the measurement of parton distribution functions (PDFs) at the LHC and in many searches for physics beyond the Standard Model (SM).From a theoretical perspective, the invariantmass distribution of the produced lepton pair is arguably the simplest hadron collider observable, and it often serves as a template to understand the structure of higher-order QCD corrections at hadron colliders more generally.It is thus important to have a solid theoretical understanding for the NCDY process, including higher-order corrections in both the strong and electroweak coupling constants.
The uncertainty on the NNLO results due to missing higher orders in the QCD perturbative expansion was estimated to be at the percent level by varying the renormalisation and factorisation scales by a factor of two around a central scale related to the invariant mass of the lepton pair.Very recently, also next-to-next-to-next-to-leading order (N 3 LO) corrections in the strong coupling have been computed, but only for the contributions from an intermediate off-shell-photon [23].Unlike for N 3 LO corrections to inclusive Higgs production processes [24][25][26][27][28][29][30], it was found that for large ranges of invariant masses, the N 3 LO corrections are sizeable and lower the value of the cross section by a few percent, which is more than the missing higher order uncertainty estimated from varying the perturbative scales at NNLO.This result was recently confirmed by the independent computation from the double-differential computation for lepton-pair production via an off-shell photon at N 3 LO in ref. [31] and the computation of the fiducial cross setion at N 3 LO in ref. [32].A similar behaviour was observed for the production of a lepton-neutrino pair at N 3 LO [33].This shows that N 3 LO corrections are highly needed if we want to reach a precision at the percent level for vector-boson production at hadron colliders.In particular, a complete calculation of the NCDY process in the SM at N 3 LO, including also the contributions from the exchange of the Z-boson, is highly desired.
The computation of higher-order corrections including Z-bosons, however, is much more challenging.Unlike the photon, the Z-boson has both vector and axial couplings to fermions.Higher-order computations typically diverge, and they are conventionally regularised using dimensional regularisation, where the space-time dimension is analytically continued from D = 4 to D = 4 − 2 dimensions.Axial couplings involve a γ 5 matrix, which is an intrinsically four-dimensional object, and so its treatment in dimensional regularisation is ambiguous (see ref. [34] and references therein for a review).Moreover, it is well known that the axial current is anomalous in QCD, and its divergence is described by the famous Adler-Bell-Jackiw (ABJ) anomaly equation [35][36][37] (this anomaly is absent for the vector current, as a consequence of Yang's theorem [38]).In the complete SM, the axial anomaly cancels.Calculations involving massive top quarks, however, are technically extremely challenging, and therefore computations are usually done in an effective theory where the top quark is infinitely heavy and is integrated out.This procedure does no longer naively work in the presence of an axial current, because the resulting effective theory is anomalous [39].
The goal of this paper is to present for the first time complete results for the invariant mass distribution for the NCDY process in the SM at N 3 LO in the strong coupling constant, including contributions from both an intermediate photon and Z-boson.Our main contribution is the computation the inclusive N 3 LO cross section for an axial-vector state in QCD.We treat the γ 5 matrix by working in the Larin scheme [40][41][42][43][44] and we work in a theory with five massless active quark flavours where the top quark has been integrated out.Since this effective theory is anomalous, the top quark does not completely decouple.We include non-decoupling effects via a Wilson coefficient multiplying the axial current in the effective theory.As a result, we obtain the complete N 3 LO cross section in the SM, up to terms that are power-suppressed in the top quark mass.The finite top-mass effects are known to be small at NNLO, and we expect the size of the missing power-suppressed terms to be negligible at N 3 LO as well.As a result, we obtain for the first time complete phenomenological predictions for the NCDY process to third order in the strong coupling.
This paper is organised as follows: In section 2 we review the structure of the QCD corrections to the NCDY process.In section 3 we present the main ingredient of our computation, namely the cross section for axial-vector production at N 3 LO in QCD, and we discuss our treatment of the γ 5 matrix in dimensional regularisation and the non-decoupling top-mass contributions.In section 4 we discuss the phenomenological implications of our results.In section 5 we draw our conclusions.

The Neutral-Current Drell-Yan process
The Neutral-Current Drell-Yan (NCDY) process describes the production of a pair of (massless) leptons as the result of the collision of (anti-) protons: P(P 1 ) + P(P 2 ) → l l(Q) + X . (2.1) Here, P 1 and P 2 are the momenta of the scattering protons, Q is the invariant mass of the lepton pair l l, and X collectively denotes additional QCD radiation.The goal of this paper is to compute this process to third order in the strong coupling constant and to discuss the associated phenomenology.More precisely, we want to compute the cross section for the process in eq.(2.1) differentially in the invariant mass Q.The production cross section is mediated by a virtual photon or Z-boson which subsequently decay to the final-state leptons.Our formalism is accurate up to corrections that are suppressed in the electroweak coupling constant α EW : Figure 1 schematically shows the NCDY process for the production of a lepton pair via an intermediate vector boson.The cross section can be written in a factorised form, We now describe the factors in the right-hand side of eq.(2.3) in detail.
The propagation of the virtual gauge bosons is described by the Breit-Wigner distribution, The decay of the gauge bosons into leptons is captured by the interference contributions to the width: x 5 y n v x 3 r 2 P e e u K l 8 8 c w R 9 4 n z 9 V 6 4 7 1 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x X e P 4 f x 5 y n v x 3 r 2 P e e u K l 8 8 c w R 9 4 n z 9 V 6 4 7 1 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x X e P 4 f x 5 y n v x 3 r 2 P e e u K l 8 8 c w R 9 4 n z 9 V 6 4 7 1 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x X e P 4 f  We use the notation M 1 •M * 2 to indicate that the interferences are summed over the colours and spins of all external particles.Equation (2.5) describes the interferences of the decay amplitudes of gauge bosons into two final-state particles with masses m l and ml at Born level.It is chosen such that if the formula is evaluated on-shell, Q = m B and B = B , then it corresponds to the partial width Γ B→l l of the gauge boson B in its rest-frame.Since we are working to leading order in the electroweak coupling, we only need to consider the tree-level contributions to the decay.For massless leptons, we find: The charges in the above equation for a Z boson or photon are defined as : Here, m W and m Z are the masses of the W and Z bosons.In the SM we have the following charge assignments for up-type and down-type quarks and electrically-charged leptons: Finally, the first factor in eq.(2.3) describes the production cross section for the interference of two off-shell gauge bosons B and B : (2.9) In the previous equation we introduced the parton distribution functions f i (x) that are convoluted with the partonic coefficient functions η ij→B/B +X .We suppress the dependence of all quantities on the factorisation scale to improve the readability.The variables x 1,2 represent the fraction of the momenta of the protons carried by the initial-state partons.
We introduce the variables (2.10) The partonic coefficient functions are given by where the factor N ij represents a process-dependent averaging over initial-state spins and colours, dΦ B/B +X i is the phase space measure for a particular final state, and M ij→B+X i is the scattering matrix element for the production of this final state.In order to compute the partonic coefficient functions to a given order in perturbation theory, we consider a perturbative expansion of the product of matrix elements truncated at a fixed order in the coupling constant.The partonic coefficient functions are the main ingredients to our computation, as they incorporate the entirety of the higher-order QCD corrections.We will study their structure in more detail in the remainder of this section.
The coupling of the Z-boson to fermions involves a vector and an axial-vector part (see figure 2).Accordingly, we may split the partonic coefficient function into two terms Z < l a t e x i t s h a 1 _ b a s e 6 4 = " E 8 f Q K Y C 8 A W f q z P w k v p S S F / 7 e X P A = " > A A A C M H i c b Z D L S s N A F I Y n 9 V 5 v V Z d u B o u g Y E t S v G 0 E L w t d K t h W b N o w m U 7 a w Z k k z J w I J e S R 3 P g o u l F Q x K 1 P 4 f S y 0 O o P A x / / O Y c z 5 / d j w T X Y 9 q u V m 5 i c m p 6 Z n c v P L y w u L R d W V m s 6 S h R l V R q J S N 3 4 R D P B Q 1 Y F D o L d x I o R 6 Q t W 9 + / O + v X 6 P V O a R + E 1 9 G L W l K Q T 8 o B T A s b y C u d H J e w G i t C 0 g j m W X j 1 L 7 z O 3 Q 6 Q k L V c m r m A B b F 1 5 a W 0 n y F q 3 J U M n A 8 L D H m / P V b z T h W 2 v U L T L 9 k D 4 L z g j K K K R L r 3 C k 9 u O a C J Z C F Q Q r R u O H U M z J Q o 4 F S z L u 4 l m M a F 3 p M M a B k M i m W 6 m g 4 M z v G m c N g 4 i Z V 4 I e O D + n E i J 1 L o n f d M p C X T 1 e K 1 v / l d r J B A c N l M e x g m w k A 4 X B Y n A E O F + e r j N F a M g e g Y I V d z 8 F d M u M f G B y T h v Q n D G T / 4 L t U r Z 2 S / v X u 0 W j 0 9 H c c y i d b S B t p C D D t A x u k C X q I o o e k D P 6 A 2 9 W 4 / W i / V h f Q 5 b c 9 Z o Z g 3 9 k v X 1 D d p y q O g = < / l a t e x i t >  accounting for the coupling to the vector and axial-vector parts separately: (2.12) The inclusive cross section does not contain any non-zero partonic coefficient functions that depend on both a vector and an axial-vector coupling.We are interested in the perturbative coefficient functions computed up to a fixed order in the expansion in the strong coupling constant: where α S (µ) denotes the strong coupling constant in the MS scheme.At leading order in perturbation theory only the partonic coefficient functions with a quark -anti-quark pair in the initial state are non-vanishing: η Y, (0) The partonic coefficient functions for the production of an electroweak gauge boson can be split into 9 different contributions, corresponding to particular combinations of the charges associated with different ways of coupling the electroweak gauge bosons to the quarks.For the purpose of simplicity, we write η ij = η Y ij→B/B +X .We then have the following decomposition: Figure 3 shows representative diagrams that illustrate how the different charge structures arise.Figure 3a shows the tree-level diagram where the initial-state q q-pair annihilates and produces the gauge boson.The flavors of the initial-state quarks are connected by the gauge boson.This is made manifest by the Kronecker delta symbol δ ij multiplying the corresponding partonic coefficient function η 7 ij .In contrast, there are classes of diagrams with a quark and anti-quark of the same flavor, but the quark lines are not identical.In such diagrams the gauge boson may or may not connect the different quark lines, as is illustrated in fig.3b  , where the gauge boson is connected only to one quark line of an initial-state quark.Furthermore, it is possible that the gauge boson couples to quark lines that are disconnected from the initial state, as shown for example in fig.3f.These contributions are taken into account by η 1 ij and η 2 ij .The partonic coefficient functions η V ij→B/B +X and η A ij→B/B +X begin to differ starting from NNLO in QCD perturbation theory.At this order, it is for the first time possible to connect two different quark lines with the produced electroweak gauge boson.As long as the gauge boson is twice connected to the same quark line and all quarks are treated as massless, helicity conservation dictates that vector and axial-vector partonic coefficient functions will be identical.In particular, η 2 ij , η 5 ij , η 6 ij , η 7 ij are identical for vector and axial-vector production at all orders in perturbation theory.
The partonic coefficient functions for the production of a vector boson η V ij→B/B +X were computed at N 3 LO in QCD perturbation theory in refs.[2-9, 23, 45].The computation of the partonic coefficient functions for an axial-vector boson η A ij→B/B +X through N 3 LO in perturbative QCD is one of the main results of this article.We present analytic formulae for these partonic coefficient functions in terms of ancillary files attached to the arXiv submission of this article.In the next section we discuss in detail the computation and the structure of the partonic coefficient functions for axial-vector production in QCD.
3 Axial-vector production at N 3 LO in QCD Let us consider N 3 LO corrections to the production cross section for a color-singlet axialvector state Z A .In broad terms, our computational strategy is similar to the computation of the N 3 LO corrections to the inclusive Higgs boson production cross sections in gluon fusion [24,46] and bottom-quark fusion [26], and the charged-current and photon-only DY processes [23,33].More precisely, we generate Feynman diagrams using QGRAF [47] and perform the spinor and color algebra based on custom C++ algorithms using GiNaC [48] and FORM [49].We generate all required integrands for interferences of matrix elements for this cross section up to N 3 LO.We then reduce all real and virtual interference diagrams to a set of so-called master integrals using IBP identities [50][51][52] and the framework of reverse unitarity [13,[53][54][55][56].The master integrals were computed using the method of differential equations [57][58][59][60][61] in refs.[46,[62][63][64][65].We work in dimensional regularisation, and compute all matrix elements in D = 4 − 2 dimensions.After inserting the master integrals into our partonic coefficient functions, we renormalise UV divergences in the MS scheme.Furthermore, we remove collinear initial-state singularities via mass factorization counterterms comprised of DGLAP [66][67][68] splitting functions [69][70][71][72].Finally, we arrive at fully analytic formulae for our finite partonic coefficient functions for the production of an axial-vector boson through N 3 LO in perturbative QCD.The functions are expressed in terms of the class of iterated integrals introduced in ref. [46].While the general strategy is the same as for the N 3 LO cross section considered in ref. [24,26,33,46], there are several technical steps that distinguish axial-vector production from those other processes, due to the ambiguities in how to treat γ 5 in dimensional regularisation.In the remainder of this section we discuss our treatment of γ 5 in detail.
Consider the operator describing the coupling of an axial current in QCD to the axialvector Z A : where N f is the number of quark species.In the SM we have and N f = 6.Since the coupling of Z A to the axial current involves the γ 5 matrix, which is only well-defined in four dimensions, care is needed how this operator is analytically continued to arbitrary dimensions when working in dimensional regularisation.Moreover, it is well-known that the singlet axial current is anomalous in QCD and satisfies the Adler-Bell-Jackiw (ABJ) anomaly equation [35][36][37] (see also ref. [73]).In four dimensions γ µ and γ 5 anticommute, as one can easily verify by writing γ 5 as in refs.[40,41]: However, if the spacetime dimension is extended to D-dimension this no longer holds true.
Here we work in the Larin-scheme [42][43][44] and define the axial-vector current explicitly as For D = 4, this definition is identical to the one of eq. ( 3.1) as on can easily see by computing the anti-commutation relation In our computation the γ µ matrices are interpreted as D-dimensional objects, and the Dirac algebra is performed in D-dimensions.Any Feynman diagram in the Drell-Yan process mediated by an axial-vector current will involve two insertions of the axial-vector current of eq.(3.3).We use the identity which is valid in strictly D = 4 dimensions, to contract Lorentz indices of the Levi-Civita tensors, and we treat the metric tensors in the above equation as D-dimensional.The above extension to D space-time dimensions modifies the computed cross sections at sub-leading order in the dimensional regulator.In the presence of divergences, these modifications are propagated into finite and singular terms of the bare partonic cross sections and their renormalisation below.
It is well-known that the axial current J µ A,f is anomalous in QCD and develops a UVdivergence starting from one-loop order.This UV-divergence cancels in the complete SM, due to the relation In variants of the SM with an odd number N f of fermion species, however, eq.(3.6) is violated, and the UV-divergence does not cancel.The UV-divergence can be removed by renormalising the axial current by the equation (valid through at least three loops) where J µ A,f R is the renormalised axial current.We see that the renormalisation mixes the axial currents for different flavors.The non-singlet and singlet counterterms Z ns and Z s are not pure MS counterterms, but they also include finite terms whose purpose is to ensure that the renormalised singlet axial current computed in the Larin-scheme satisfies the all-order ABJ anomaly equation [35][36][37]: with F F = µνρσ Tr F µν F ρσ , and the renormalised singlet axial current is Both the non-singlet and singlet counterterms are known to three loops in QCD [43,[74][75][76]: We have computed the partonic coefficient functions in the Larin-scheme in a variant of the SM with only N f = 5 massless active quark flavors, and the top quark is considered infinitely-heavy and thus absent from the computation.If only the strong coupling constant is UV-renormalised, then the partonic coefficient functions still exhibit poles in the dimensional regulator even after mass factorisation.However, all these poles cancel once the axial current is renormalised as in eq.(3.7), and we obtain finite results for all partonic coefficient functions.Besides the explicit analytic cancellation of all ultraviolet and infrared poles, we performed several additional checks to validate our results.In particular, we find agreement with our previous computation of ref. [23] of partonic coefficient functions for vector boson production for the part of the partonic coefficient functions that are identical between vector and axial-vector production.Moreover, purely virtual corrections were computed in ref. [75,[77][78][79] and we find agreement.Finally, we also checked that our results have the expected dependence on the perturbative scale µ = µ F = µ R .This check, however, involves a subtle point, which we discuss in detail in the remainder of this section.
Our goal is to compute the coefficient functions η A ij , which correspond to the axialvector contribution from Z-boson exchange in the limit where we drop all power-suppressed terms in the top-quark mass.If we denote the partonic coefficient functions obtained from our procedure by η A,5 ij , we find that generically η A,5 ij = η A ij .To see this, consider the dependence of the coefficient functions on the perturbative scale µ := µ F = µ R .In the SM, this dependence is governed by the the well-known Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation [66][67][68]: where β = − ∞ n=0 β n a n+1 s denotes the QCD β-function and P kl is the DGLAP splitting kernel, and we introduced the convolution: (3.12) The coefficient functions that we computed, however, satisfy the evolution equation where we defined implicitly and we find it useful to introduce the quantity N A = N f f =1 A f .Note that N A = 0 in the complete SM with N f = 6 quark species.γ J is the anomalous dimension of the (renormalised) axial current [43,75]: Here we defined the matrix where I N f is the N f × N f unit matrix and E N f is the N f × N f matrix whose entries are all 1.We have Clearly, the term proportional to γ J in eq.(3.13) is absent in the complete SM, and so η A,5 ij cannot be identified with the large-m t limit of the SM.Note that in the complete SM we have N A = 0, and so γ J = 0.The interpretation and resolution of this puzzle is as follows: If we start from the complete SM, with N f = 6 and m t < ∞, the UV-divergences from the axial anomaly cancel.Naively, one would expect that by making the top quark infinitely heavy, it decouples from the theory and we land on an effective field theory with N f = 5 massless quarks, defined by simply ignoring operators involving the top quark.This naive expectation, however, is wrong: Since the UV divergences from the axial current do not cancel for N f = 5, our EFT is anomalous and the top quark does not naively decouple (see, e.g., ref. [39]).We include non-decoupling effects in the form of a Wilson coefficient The Wilson coefficient satisfies a renormalisation group equation which precisely cancels the contribution from the anomalous dimension γ J in eq.(3.13): We see that the partonic coefficient functions η A ij in the large-m t limit of the SM are obtained from η A,5 ij by the replacement A f → C A,f (µ) (cf.eq.(3.14)): The Wilson coefficient admits the perturbative expansion: It is easy to check that by combining the evolution equations (3.13) and (3.19), the function η A ij defined in eq.(3.20) satisfies the SM evolution equation (3.11).
The Wilson coefficient depends explicitly on the (on-shell) top-quark mass, making the non-decoupling of the top quark manifest.It was recently computed through three loops in ref. [76,80].We have independently obtained all the coefficients except for the non-logarithmic three-loop coefficient c (3,0) A , and we find full agreement.Indeed, the coefficients of the logarithms, i.e., the c (n,k) A with k = 0, are determined from the evolution equation (3.19).The non-logarithmic terms c (n,0) A are not constrained by the evolution equation and need to be computed explicitly.We have recalculated the two-loop coefficient c (2,0) A by taking the large-m t limit on the known NNLO cross section in the full SM, including all finite quark-mass effects [10][11][12].The two-loop coefficients are: At three-loops, we have [76,80]:

Phenomenological Results
In this section we present for the first time phenomenological results for the invariant-mass distribution of a massless lepton pair at N 3 LO in QCD, We work in the five-flavour scheme with N f = 5 active massless quark flavors.The topquark is considered infinitely heavy, though the non-decoupling effects are included through the Wilson coefficient C A,f .We mention that finite quark-mass effects have been calculated through NNLO and are known to be very small [10][11][12].We therefore expect our computation to give a very good estimate of the QCD corrections in the full SM where the full finite top-mass effects are retained.Our choice for the numerical values of the SM input parameters is summarized in table 1.The strong coupling constant is evolved from α S (m Z ) to the renormalisation scale µ R using the four-loop QCD beta function [81][82][83][84] in the MS-scheme using N f = 5 active massless quark flavors.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 [85], and bands correspond to varying the perturbative scales µ F and µ R by a factor of two around the central scale Commonly, this is referred to as 7-point variation.

Scale dependence and perturbative convergence
In table 2 we show values for a selection of representative invariant masses Q 2 for the choice of central scale µ F = µ R = Q, together with the corresponding QCD K-factors: We include an estimate of the residual perturbative uncertainty based on seven point variation of the factorization and renormalisation scale, as well as the uncertainties from PDFs and the value of the strong coupling constant (see below).In figure 4 we show the inclusive cross section for the production of a massless lepton pair as a function of Q 2 .In figures 5 and 6 we show the NCDY cross section normalized to its value computed at N 3 LO as a function of the invariant mass Q.We observe qualitatively the same features as for the photon-only [23,31] and the charged-current DY processes [33].Specifically, we observe that in the range of invariant masses between ∼ 40 GeV and ∼ 400 GeV, the scale variation bands from NNLO and N 3 LO do not overlap, indicating that conventional scale variation at NNLO underestimated the true size of the N 3 LO corrections.We note, however, that the size of the bands at NNLO was particularly small for the NCDY process, often at the sub-percent level depending on the invariant masses considered.

LO
In figure 7 we show the dependence of the cross section for Q = 100 GeV on one of the two perturbative scales with the other held fixed at some value in the interval [Q/2, 2Q].We observe a very good reduction of the scale dependence as we increase the perturbative order, with only a very mild scale dependence at N 3 LO.Just like for the photon-only and W cases, the bands from NNLO and N 3 LO do not overlap.

Ratios of K-factors
In the previous section we have seen that the NCDY process shows relatively good perturbative stability as we increase the order.This behavior is very reminiscent of the photon-only and charged-current DY processes considered in refs.[23,33].To investigate if and to what degree QCD corrections differ among different DY-type processes, we study ratios of K-factors between the NCDY process and the photon-only and charged-current DY processes.In particular, in this section we refer to the ratio of the N n LO cross section to its LO counter part as the N n LO K-factor, K (n) .In contrast, ratios of cross sections rather than K-factors involving the W and γ * cross sections were considered in ref. [33].Before we discuss these ratios, we need to make a comment.Whenever one studies ratios of cross sections, there is an ambiguity in how to choose the perturbative scales.For example, if one believes that QCD corrections should be similar between W , γ * and Z production (as motivated for example by the universality of certain limits, like the threshold limit), it is natural to vary the scales in the numerator and the denominator in a correlated way.Alternatively, one may choose the scales in an uncorrelated way (e.g., because different partonic channels are weighted differently by the PDFs for these processes, breaking the universality of the QCD corrections), typically leading to larger scale variation bands.In ref. [33] it was shown that for the W and γ * cross sections the correlated prescription leads to a vanishingly small scale dependence at the sub-permille level at N 3 LO, while the uncorrelated prescription produces excessively large scale variation bands at N 3 LO (much larger than the absolute shift from NNLO to N 3 LO).As a consequence, neither the correlated nor the uncorrelated prescriptions are expected to give reliable estimates for missing higher-order terms for these ratios at N 3 LO.Therefore, in ref. [33] a new prescription was considered, which uses the relative size of the last considered order compared to the previous one as an estimator of the perturbative uncertainty: In the following we use this prescription to obtain uncertainty bands for the ratios.
In figures 8 -11 we show the ratios of the K-factors as a function of the invariant mass Q.In all cases we observe a remarkable similarity of the K-factors, which agree among themselves within at worst 2% for all ratios considered.In particular, we see from figure 11 that the K-factors agree within 1% for the NCDY process computed with or without including the contributions from the Z boson.At the same time, we observe that there is a dependence of the shape of the QCD corrections on the invariant mass, and different invariant-mass regions may receive slightly different QCD corrections.This shows that, if we want to reach a level of precision of 1%, care is needed when using K-factors obtained from one process in one region of phase space to reweight other processes or other regions of phase space.However, our results also demonstrate that shape differences mainly introduced at lower orders in perturbation theory and higher order corrections are remarkably similar.

Uncertainties related to PDFs
In order to assess the dependence of our predictions on the methodology of how the PDFs are extracted, we follow the prescription of ref. [85] for the computation of PDF uncer-  Q tainties δ(PDF) using the Monte Carlo method.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 α down S = 0.1165.These sets allow us to deduce an uncertainty δ(α S ) on our cross section following the prescription of ref. [85].We combine the PDF and strong coupling constant uncertainties in quadrature to give

Q [GeV]
Currently there is no available PDF set extracted from data with N 3 LO accuracy, and so we are bound to use NNLO PDFs in our predictions.We estimate the potential impact of this mismatch on our results using the prescription introduced in ref. [25].The PDF theory (PDF-TH) uncertainty is then obtained by studying the variation of the NNLO cross section as NNLO-or NLO-PDFs are used: Here, the factor 1 2 is introduced as it is expected that this effect becomes smaller at N 3 LO compared to NNLO.

Q [GeV]
Relative Uncertainty [%] LHC @ 13 TeV PDF4LHC15 P P  e + e -+X Figure 12: Relative uncertainty of the NCDY process at N 3 LO due to incomplete knowledge of parton distribution functions and the strong coupling constant as a function of Q. δ(PDF), δ(PDF + α s ) and the sum of δ(PDF + α s ) and δ(PDF-TH) are shown in red, brown and green respectively.
In figure 12 we show the combined uncertainty from PDFs, the value of the strong coupling constant α S and the missing N 3 LO PDFs.The size of these uncertainties is comparable to the uncertainties obtained in refs.[23,33] for the photon-only and chargedcurrent DY processes.ferent PDF sets.PDF4LHC15 is a combination of the CT14 [86], MMHT14 [87] and NNPDF3.0[88] PDF sets and we show in fig.13 predictions based on these individual sets relative to the prediction based on the PDF4LHC15 nnlo mc set.The red band in fig.13 reflects the δ(PDF) uncertainty of PDF4LHC15 and we observe that the predictions based on the individual PDF sets are contained within this band and that their spread is comparable in size to this band.Since the publication of the PDF4LHC15 combination a plethora of developments and the inclusion of LHC data into global PDF fits has led to updated PDF.
In fig.14 we study in particular the PDF sets ABMP16 [89], CT18 [90], MSHT20 [91], NNPDF3.1 [92] and NNPDF4.0[93].We observe, that the spread among the newer sets is less than the the spread of their predecessors.In particular, in the range where Q is comparable to the Z and W -boson masses the different PDF sets seem to agree nicely.We observe furthermore that NNPDF4.0 leads to a significant enhancement of the NCDY cross section at large values of Q.It would be interesting to study the impact of new PDF sets on high precision processes in further detail in the future and we are looking forward to an updated version of PDF4LHC [94].

Contributions from interferences
We conclude this phenomenological study by investigating some properties of the NCDY cross section at N 3 LO and how it receives contributions from the photon, the Z-boson and their interference.
In figures 15 and 16 we show the relative size of the three contributions.We see that, as expected, the photon-only contribution dominates below the Z-threshold.Above threshold all contributions are of similar size.Close to the threshold the cross section is dominated entirely by the Z boson, and the interference changes sign when the threshold is crossed, as expected.The observed pattern is hardly modified by perturbative corrections.
Computations with γ 5 matrices are technically more involved than for purely vectorial couplings.For this reason, the axial-vector contribution is frequently approximated by simply using the same partonic coefficient functions as for the vector contribution and  simply changing coupling of the quarks in an appropriate fashion.While this is correct for contributions from Dirac traces with an even number of γ 5 matrices, a mismatch is introduced starting from NNLO for contributions involving traces with an odd number of γ 5 matrices.Moreover, effects from the axial-anomaly and the resulting non-decoupling of the top-quark are neglected in this way.

Conclusions
In this paper we have computed for the first time the complete inclusive N 3 LO corrections in the strong coupling constant for the production of a massless lepton pair.This result has been made possible by combining our computation for the N 3 LO corrections to vector production from ref. [23] with the computation of the axial-vector production cross section.A main ingredient in our computation is the use of the Larin scheme to treat the γ 5 matrix in dimensional regularisation.Our computation has passed several non-trivial checks: Apart from the explicit cancellation of all ultraviolet and infrared poles, we have checked that our implementation of the Larin scheme reproduces the contributions to the axial-vector cross section where both γ 5 matrices appear in the same Dirac trace, and that our final result satisfies the expected DGLAP renormalisation group equation once the nondecoupling top-mass effects are included.We include all partonic coefficient functions for the Drell-Yan production cross section as described by eq.(2.15) as ancillary files together with the arXiv submission of this article.We have studied the phenomenological impact of our computation.We find that the complete N 3 LO corrections to the NCDY process, including the axial-vector contributions, has the same features as for the photon-only and W cases of refs.[23,33].We find that the dependence of the complete NCDY process on the perturbative scales is similar to the production of a γ * or W .In particular, the bands obtained by varying the perturbative scales by a factor of two do not overlap when going from NNLO to N 3 LO, showing once more the importance of considering N 3 LO corrections for precision LHC observables.We have also considered ratios of K-factors between γ * , W and the complete NCDY process, and we find that, while the K-factors are very similar, the ratio depends on the invariant mass of the lepton pair, reaching a few percent depending on the invariant masses considered.This shows that care is needed when taking K-factors from one process and to apply it to another process, especially when aiming for precision results.Finally, we have also studied the impact of the choice of the PDF on our results.In particular, we have analysed how the missing N 3 LO impact our predictions, which is one of the main bottlenecks to further improve theoretical predictions for LHC observables.Since the NCDY process is one of the main measurements used to constrain PDFs from LHC data, we expect that our computation will play an important role in the future to determine precisely the proton structure.

2 < l a t e x i t s h a 1 _ b a s e 6 4 =
0 n j v B I G l f D + o l y 9 z u M o w D G c w B m E c A l V u I M a 1 I H A I z z D K 7 x 5 y n v x 3 r 2 P e e u K l 8 8 c w R 9 4 n z 9 V 6 4 7 1 < / l a t e x i t > dˆ dQ " d v + 9 O t w 1 8 R i 2 s x x O G / 6 + M p 2 D j I A = " > A A A C D 3 i c b V D L S s N A F J 3 4 r P U V d e l m s C

Figure 1 :
Figure 1: Schematic depiction of the Drell-Yan process and its factorisation into the production probability of a virtual gauge boson and the subsequent decay to final-state leptons.
< l a t e x i t s h a 1 _ b a s e 6 4 = " S n Y Q S 6 1 M 9 Q i h K l q m d y E 6 x Q 3 C + L w = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N U Y 9 E L x 4 h k U e E D Z k d e m F k d n Y z M 2 t C C F / g x Y P G e P W T v P k 3 D r A H B S v p p F L V n e 6 u I B F c G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q e N U M W y w W M S q H V C N g k t s G G 4 E t h O F N A o E t o L R 7 c x v P a H S P J b 3 Z p y g H 9 G B 5 C F n 1 F i p / t A r l t y y O w d Z J V 5 G S p C h 1 i t + d f s x S y O U h g m q d c d z E + N P q D K c C Z w W u q n G h L I R H W D H U k k j 1 P 5 k f u i U n F m l T 8 J Y 2 Z K G z N X f E x M a a T 2 O A t s Z U T P U y 9 5 M / M / r p C a 8 9 i d c J q l B y R a L w l Q Q E 5 P Z 1 6 T P F T I j x p Z Q p r i 9 l b A h V Z Q Z m 0 3 B h u A t v 7 x K m h d l 7 7 J c q V d K 1 Z s s j j y c w C m c g w d X U I U 7 q E E D G C A 8 w y u 8 O Y / O i / P u f C x a c 0 4 2 c w x / 4 H z + A L u v j O c = < / l a t e x i t >

< l a t e x i t s h a 1 _
b a s e 6 4 = " m X T T q J L V G z Q H k Y p t y o i k s c q K U l E = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W Ol Z t g v V 9 y q O w d Z J V 5 O K p C j 0 S 9 / 9 Q Y x S y O U h g m q d d d z E + N n V B n O B E 5 L v V R j Q t m Y D r F r q a Q R a j + b H z o l Z 1 Y Z k D B W t q Q h c / X 3 R E Y j r S d R Y D s j a k Z 6 2 Z u J / 3 n d 1 I T X f s Z l k h q U b L E o T A U x M Z l 9 T Q Z c I T N i Y g l l i t t b C R t R R Z m x 2 Z R s C N 7 y y 6 u k f V H 1 L q u 1 Z q 1 S v8 n j K M I J n M I 5 e H A F d b i D B r S A A c I z v M K b 8 + i 8 O O / O x 6 K 1 4 O Q z x / A H z u c P z d + M 8 w = = < / l a t e x i t > f < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 O 7 9 j k R F Q 0 3 0 b G v n j g Y g i b a y A U w = " > A A A B 7 X i c b V D L T g I x F L 3 F F + I L d e m m k Z i 4 I j O G q E u i G 5 e Y y C O B C e m U D l Q 6 7 a T t m J A J / + D G h c a 4 9 X / c + T c W m I W C Z 3 V y z r 2 5 5 5 4 w E d x Y z / t G h b X 1 j c 2 t 4 n Z p Z 3 d v / 6 B 8 e N Q y K t W U N a k S S n d C Y p j g k j U t t 4 J 1 E s 1 I H A r W D s e 3 M 7 / 9 x L T h S j 7 Y S c K C m A w l j z g l 1 k m t X k g 0 j v r l i l f 1 5 s C r x M 9 J B X I 0

Figure 2 :
Figure 2: Feynman rule coupling a Z boson to fermions.

Figure 3 :
Figure 3: Representative diagrams for different charge structures appearing in the NCDY cross section.The diagrams show interference diagrams where the vertical dashed line represents the final state phase space cut.Solid, wavy and curly lines represent propagating quarks, gauge bosons and gluons respectively.

Figure 4 :
Figure 4: The invariant-mass distribution Σ(Q 2 ) at the LHC with √ S = 13 TeV at different orders in perturbation theory.

Figure 5 :Figure 6 :1
Figure 5: The K-factors Σ N k LO /Σ N 3 LO as a function of invariant masses 10 GeV≤ Q ≤150 GeV for k ≤ 3. The bands are obtained by varying the perturbative scales by a factor of two around the central µ cent.= Q.

Figure 7 :
Figure 7: Dependence of the invariant-mass distribution through N 3 LO on one of the two perturbative scales with the other held fixed.The bands are obtained by varying the other scale by a factor of two around the central scale Q = 100 GeV.

Figure 8 :
Figure 8: Ratio of the K-factors for W + and W − production as a function of Q 2 at different orders in perturbation theory.

Figure 9 :
Figure 9: Ratio of the K-factors for NCDY and W − production as a function of Q 2 at different orders in perturbation theory.

Figure 10 :
Figure 10: Ratio of the K-factors for NCDY and W + production as a function of Q 2 at different orders in perturbation theory.

Figure 11 :
Figure 11: Ratio of the K-factors for NCDY and γ * production as a function of Q 2 at different orders in perturbation theory.

Figure 13 :Figure 14 :
Figure 13: Dependence of the NCDY process at N 3 LO on the choice of the PDF set relative to the combined PDF4LHC15 nnlo mc set.The red band corresponds to the δ(PDF) uncertainty.

Figure 15 :
Figure 15: Decomposition of the invariant-mass distribution into the absolute value of contributions from photon and Z exchange and their interference as a function of Q at different order in perturbations theory.

Figure 16 :
Figure 16: Relative contributions to the invariant-mass distribution from photon and Z exchange and their interference as a function of Q at different order in perturbations theory.

Figure 18 :
Figure 18: Relative contribution from the non-decoupling top-mass effects as a function of Q at different order in perturbations theory.

Table 1 :
Numerical values of the SM input parameters used for our phenomenological predictions.

Table 2 :
Cross section values for a selection of representative invariant mass Q for the choice of central scale µ F = µ R = Q, together with the corresponding QCD K-factors and uncertainties (see main text for a detailed description).