Resummed Drell-Yan cross-section at N$^3$LL

We present the resummed predictions for inclusive cross-section for Drell-Yan (DY) production as well as onshell $Z,W^\pm$ productions at next-to-next-to-next-to leading logarithmic (N$^{3}$LL) accuracy. Using the standard techniques, we derive the $N$-dependent coefficients in the Mellin-$N$ space as well as the $N$-independent constants and match the resummed result through the minimal prescription matching procedure with that of existing next-to next-to leading order (NNLO). In addition to the standard $\ln N$ exponentiation, we study the numerical impacts of exponentiating $N$-independent part of the soft function and the complete $\bar{g}_0$ that appears in the resummed predictions in $N$ space. All the analytical pieces needed in these different approaches are extracted from the soft-virtual part of the inclusive cross section known to next-to-next-to-next-to leading order (N$^3$LO). We perform a detailed analysis on the scale and parton distribution function (PDF) variations and present predictions for the 13 TeV LHC for the neutral Drell-Yan process as well as onshell charged and neutral vector boson productions.


Introduction
Standard Model(SM) has been very successful so far in describing the physics of elementary particles. Precision study has played an important role in establishing the SM through the latest discovery of Higgs boson at the Large Hadron Collider (LHC). The properties of the Higgs boson is being studied with higher accuracy. Recent observations at the LHC demonstrate that the systematic precision study is essential to look for any deviation from the SM in search of new physics beyond the SM (BSM). While there is no promising sign of new physics signature so far at the LHC, it is extremely important to know the SM predictions for the standard processes like Higgs and DY or Z, W ± productions to utmost accuracy. Not only this could help in BSM searches but also help to understand the perturbative structure of the underlying gauge theory. Drell-Yan production has been a standard candle at the hadron colliders and is extremely important for luminosity monitoring. This is one of the hadronic processes which is well understood theoretically. For example, next to next to leading order (NNLO) quantum chromodynamics (QCD) correction [1][2][3] to this process was computed three decades ago. DY is also an important process experimentally for several BSM searches. Experimentally, one has a very clean environment for precise measurements in terms of the kinematics of the final state lepton pairs. Higher order perturbative QCD corrections to DY provides ample opportunity to explore the structure of the perturbation series. Thus DY serves as an important process in collider experiments. At the LHC, the strong interaction dynamics dominates over the others and hence There have been attempts to go beyond NNLO accuracy in order to improve the precision from the theoretical side.
The calculation of complete N 3 LO cross-section is extremely difficult due to increasing number of subprocesses involved, however there have been significant progress to obtain third order contribution to this process in QCD. Very recently the first result at complete N 3 LO from only virtual photon mediator has been calculated in [4]. From the theory side, DY is seen to be extremely stable with respect to factorization and renormalisation scales already at NNLO. The scale uncertainty has been seen to be reduced to 2% for a canonical variation of factorization and renormalisation scales compared to NLO where uncertainty is about 9.2%, whereas the K-factor seem to improve marginally from 1.25 at NLO to 1.28 at NNLO. However keeping in mind the importance of this process, it is worth studying the results from next orders and devise methods to incorporate more and more higher order corrections. Since a complete calculation beyond NNLO level is difficult, the soft-virtual (SV) contributions is often computed as first step. In addition, the later constitutes a significant part of the cross-section in the region where the partonic scaling variable z → 1, called the threshold region. The SV cross-sections are known for many SM processes e.g. Higgs production [5][6][7][8][9][10][11], associated production [12], bottom quark annihilation [13], pseudo-scalar Higgs production [14].For DY production, using the three loop quark form factor [15], exploring the universal structure of the soft part [16] of SV cross section to Higgs production [5], the dominant soft-virtual (SV) corrections for DY at third order was obtained [17] and later it was confirmed in [18].
The SV contributions dominate at every order in the perturbation theory through large logarithms spoiling the reliability of the fixed order predictions. The resolution to this is to resum these large logarithms to all orders. Resummation of these large logarithms is thus very important to correctly describe the cross section in the threshold region. In [19][20][21], a systematic approach was proposed to resum these logarithms to all orders. The large logarithms arise in the hard partonic cross section when the total available center-of-mass energy (ŝ) becomes close to the invariant mass (Q) of the final state, in other words the partonic scaling variable z = q 2 /ŝ → 1. This results from the soft gluon emissions, as a consequence of which the cross-section is enhanced by the large logarithms that appear as distributions namely Dirac delta δ(1 − z) and + distributions: In Mellin space these singular terms are transformed into powers of logarithms of the Mellin variable N . In Mellin N space, these contributions can be systematically resummed to all orders and they display a nontrivial pattern of exponentiation. In the threshold region, the fixed order predictions often fail to describe the cross-section well and hence the resummation of these large logarithms becomes very important to correctly describe the region. Moreover, it has been very well established that the resummed contributions give a sizeable contributions to the cross-section. In fact many SM fixed order calculations have been improved with the corresponding resummed results, for example, the inclusive scalar Higgs production in gluon fusion [6,[22][23][24][25] (see also [26] for renormalisation group improved prediction) as well as in bottom quark annihilation [27], deep inelastic scattering [28,29], DY production [6,23,30], pseudoscalar Higgs production [31][32][33], spin-2 production [34,35] etc. Threshold resummation not only improves the inclusive fixed order results but also differential observables like rapidity [20,[36][37][38][39] and in the context of LHC precision measurements, it is important to include these corrections and they are shown to improve the fixed order results.
In the resummed predictions for the cross-sections, there is an intrinsic ambiguity on what is exponentiated and what is not. In the standard approach, one exponentiates only large-N pieces coming from the soft function which are enhanced in the threshold region. However one can also define large logarithms in terms of a new variable N = N exp(γ E ), γ E being the Euler-Mascheroni constant. Theoretically this is allowed, since γ E arises as a mathematical artifact due to dimensional regularization in d-space time dimensions. Moreover, this does not spoil the fact that the large-N pieces are exponentiated in the threshold region. In this terminology, one exponentiates N instead of N . Numerically, however this makes a difference already at the leading logarithmic accuracy. It has been already seen in [29] that the perturbative convergence is improved if one exponentiates the large-N terms. Apart from the standard threshold exponentiation, one can in fact exponentiate the complete soft function i.e., all the large-N terms as well as the δ(1 − z) terms arising from the soft function. We call this 'Soft exponentiation' which renders some part of the N -independent constant (g 0 ) for the standard N -threshold. In addition to these procedures, one can also exponentiate the complete form factor along with the soft function. This was studied in the context of the SM Higgs production [24,25] and was shown to improve the scale uncertainty better than the standard threshold approach. The form factor is process dependent and therefore is non-universal unlike the soft function. However, the form factor as well as the soft function both satisfy the similar Sudakov K+G type equation [8,9,[40][41][42][43]. Hence the solution to K+G equation for the form factor is an exponential N independent constant justifying the exponentiation. The numerical impact of this has already been studied in the past for the DY production in [44] where the authors show that both in DIS scheme and in MS scheme the complete form factor exponentiates to the orders currently known.
The goal of the present article is to study the effect of threshold logarithms at N 3 LL accuracy and match it to the known NNLO results. We perform this study for the neutral DY production as well as for onshell Z and W ± productions. The paper is organized as follows. In Sec. 2.1, we collect the useful formulae required for the invariant mass distribution for DY and total cross-section for Z, W ± productions at the LHC. Next we discuss the theoretical set up in the context of resummation. Here we describe in detail the factorization of soft-virtual coefficient in Sec. 2.2. Next we set up in Sec. 2.3 different resummation prescriptions as well as derive some useful formulas needed. Sec. 3 we study in detail the effect of threshold logarithms for different prescriptions and present our results along with the estimation of uncertainties. We finally conclude in Sec. 4.

Drell-Yan and Z, W ± production
The hadronic cross-section for DY or onshell Z, W ± production at the LHC can be written as where σ = dσ dQ τ, q 2 for DY production, with Q being the invariant mass of the di-lepton pair. Here f a (x 1 , µ 2 f ) and f b (x 2 , µ 2 f ) are the non-perturbative parton distribution functions (PDFs) of the partons a, b carrying momentum fractions x 1 , x 2 of the incoming protons at the factorization scale µ f . These PDFs are appropriately convoluted with perturbatively calculable partonic coefficients ∆ ab (z, q 2 , µ 2 f ). For the on-shell Z, W ± production, σ = σ V , V = Z/W ± and Q = M V , the mass of the vector boson. The partonic coefficients are obtain from the partonic cross section using perturbation theory. For the DY production V we include contributions from γ and Z as well as their interference.
The partonic cross-section can be decomposed as 2) The first term ∆ (sv) is called the SV partonic coefficient and it contains distributions such as δ(1 − z) and D + , whereas the second term ∆ (reg) contains those terms that are regular in the scaling variable z. The prefactors for DY and Z, W ± production are given as below: where S is the hadronic centre-of-mass energy and n c = 3 in QCD. For DY production, the factor F (0) is found to be, Here α is the fine structure constant, c w , s w are sine and cosine of Weinberg angle respectively. M Z and Γ Z are the mass and the decay width of the Z-boson.

5)
Q a being electric charge and T 3 a is the weak isospin of the electron or quarks. In the threshold region, the SV terms which consist of distributions contribute significantly at the hadron level. After mass factorization, the partonic coefficient in the threshold region experiences further factorization in terms of the form factor and soft-collinear function. In the next section we will discuss in detail on the structure of distributions in the SV coefficient which will form the basis for the resummation.

Soft-virtual cross-section
In the following, we briefly describe the theoretical set up that is required to study the impact of threshold corrections within the framework of resummation a la Sterman, Catani and Trentedue [19,20]. We do this in order to understand the role of various pieces that contribute to the resummed result. Exploiting the factorization of infrared sensitive contributions and gauge and renormalisation group invariances, inclusive cross section for the DY and on-shell Z/W ± productions in the threshold limit can be expressed in terms of form factor of the neutral/charged current, soft distribution function and Altarell-Parisi kernels (see [8,9]). The resulting expression expressed in z space is free of both ultraviolet and infrared divergences and captures the distributions D j with given logarithmic accuracy to all orders in perturbation theory. In the Mellin N space, we can achieve the same and in addition, one has the advantage to reorganize the series in such a way that order one contributions of the form a s β 0 log(N ) can be resummed systematically to all orders in the large N limit. Here, a s is defined by a s = g 2 s (µ 2 r )/16π 2 with g s begin the strong coupling constant and µ r the renormalisation scale and β 0 is the first coefficient of the coupling constant beta function. Note that in Mellin N space, the convolutions in z space become simple products. The z space result can be used to compute the soft-virtual contributions in power series expansion of strong coupling constant a s .
In d = 4 + space time dimensions, the threshold enhanced partonic soft-virtual crosssection to all orders in perturbation theory in z space can be written [8,9] as Here Ψ is a distribution function which is finite in the limit → 0. The symbol C denotes the Mellin convolution (denoted below as ⊗) which in the above expression should be treated as with f (z) being a function containing only δ(1 − z) and plus distributions. The finite exponent in the above is refactorized in the threshold limit and gets contribution from the form factor F (â s , Q 2 , µ 2 , ) with q 2 = −Q 2 , soft-collinear function Φ(â s , z, q 2 , µ 2 , ) (later called as soft function) as well as mass factorization kernels Γ(â s , z, µ 2 f , µ 2 , ) and takes the following form in dimensional regularization: µ keeps the strong coupling (â s ) dimensionless in the d = 4 + dimensions. Z(â s , µ 2 r , µ 2 , ) denotes the overall UV renormalization constant which for the processes considered here is unity due to conserved current.
The bare quark form factor satisfies the Sudakov K+G equation [8,9,[40][41][42][43] which follows as a consequence of the gauge invariance as well as renormalisation group invariance, The function K contains all the infrared poles in whereas the function G is finite in the limit → 0. The renormalisation group invariance leads to the following solutions of these functions in terms of cusp anomalous dimensions (A): The cusp anomalous dimensions are known to fourth order [45-47, 47-49, 49-56] and are collected in Appendix C. The µ r independent piece of the G can be written in perturbative series as where the coefficients G (j) ( ) can be decomposed as where The coefficients G ik are the finite coefficients found in terms of QCD color factors and can be extracted from explicit calculation of quark form factor. Note that up to the third order one also needs coefficients G 22 , G 31 and thereby needs the three-loop calculation of the form factor [15]. We have collected them in the Appendix C. Similar to the cusp anomalous dimension, the coefficients f i have been found to be maximally non-abelian to third order in strong coupling i.e. they satisfy (2.14) The initial state collinear singularities are removed using the Altarelli-Parisi (AP) splitting kernels Γ(â s , µ 2 f , µ 2 , z, ). They satisfy the well-known DGLAP evolution given as, where P (z, µ 2 f ) is the AP splitting functions. The perturbative expansion for these splitting functions has the following form: As already discussed, only the qq channel contributes to the SV cross-section and thus we find that, only the diagonal terms of the splitting functions contribute to the SV crosssection. The diagonal part of the splitting functions is known to contain the δ(1 − z) and distributions and can be written as, The splitting functions are known exactly to four loops [45,[57][58][59].
The finiteness of the soft-virtual cross-section demands that the soft-collinear function Φ will also satisfy similar Sudakov type equation like the form factor i.e. one can write where K(â s , z, µ 2 r µ 2 , ) contains all the poles and G(â s , z, q 2 µ 2 r , µ 2 r µ 2 , ) is finite in the dimensional regularization such that Ψ becomes finite as → 0. The solution to the above equation has been found [8,9] Φ (j) can be found from the solution of the form factor by the replacement as A → −A, G( ) → G( ). Notice that G( ) are now new finite z-independent coefficients coming from the soft function whereas the z dependence has been taken out in Eq. (2.19). This can be found by comparing the poles and non-pole terms inΦ (j) with those coming from the form factors, overall renormalisation constants, splitting kernel and the lower order SV terms. The coefficient G has same structure as the form factor in Eq. (2.12) after (2.21) The coefficients f i are same as those appear in the quark form factor in Eq. (2.12). The coefficientsG ij required up to three loops have been extracted in [60] and also collected in the Appendix C. Note that one has to perform the following expansion in Eq. (2.19) in order to get all the distributions and delta function coming from the soft function, It is worth noting that G as well as the complete soft function Φ I satisfy the maximally non-abelian property up to three loops. Moreover Φ I is also universal in the sense that it only depends on the initial legs and is completely unaware of the color neutral final state. Expanding ∆ (sv) in powers of a s as with the born contribution being ∆ (0) = δ(1 − z). The SV correction at the three loops are known [17] which we collect here for completeness in the Appendix A.
In the following, we will study the numerical impact of resummed result resulting from ∆ (sv) ab after performing the Mellin transformation in the large N limit. We start with Ψ which is finite while the individual contributions to it contain UV and IR singularities. Decomposing the later ones as sum of singular and finite parts as Substituting the above equations in Eq.(2.8), we can easily show that all the singular terms in the limit → 0 cancel among themselves. In addition, D 0 terms in finite part of C ln Γ go away when added to Φ fin D resulting in a finite distribution. Substituting the Ψ in Eq.(2.6), we obtain where (supressing dependence on µ f and µ r ), the N independent constant C 0 is given by So far, we showed how various collinear soft gluon emissions as well as the wide angle soft emissions can be systematically summed to all orders to obtain Eq.(2.25) in the z space when partonic variable z → 1. Note that C 0 is obtained by first collecting those terms that are proportional to δ(1 − z) terms of Ψ and then expanding the exponential of them in powers of a s . The remaining function G + contains only distributions D j . Hence, one can predict the following structure for G + : (2.28) where each G 1 sums certain terms of the a i s D i−1 to all orders, and G 2 sums a i s D i−2 terms to all orders, etc etc. The result ∆ sv expressed in terms of C 0 and the exponential of G + using Eq.(2.28) systematically sums the distributions D j to all orders and hence can predict these distributions to all orders provided A and D are known to desired accuracy in a s . For example, knowing A 1 , we can predict all the terms a i s D i with i = 1, 2, ..., ∞ in Φ, similarly given A 1 and D 1 , we can predict a i s D i−1 with i = 1, 2, ..., ∞ etc. Hence, expression given in Eq.(2.25) has the predictive power for ∆ sv to all orders in a s given the logarithmic accuracy in z space, quantified by terms of the form a i s D j . Note that when the exponential of Φ is expanded using convolution rules given in Eq.(??), we will get not only D j but also δ(1 − z). In other words, δ(1 − z) terms in ∆ sv can come from both exp(G + ) as well as C 0 .
Often in certain kinematic regions, these contributions can be enhanced when convoluted with the parton distribution functions spoiling the reliability of the perturbation theory. Hence we need to include these potentially large terms to all orders in perturbation theory for any sensible predictions. Such an exercise in the z space is technically challenging due to the complexity involved in computing the convolutions of D j . However, in the Mellin N space, the convolutions become simple products allowing us to study the impact of these large logarithms to all orders in a systematic fashion. In the following, we will describe how this can be done in Mellin N space.

Threshold Resummation
In the last sub-section, we showed that threshold effects for partonic coefficients can be obtained near threshold as a product of well-defined functions, each organizing a class of infrared and collinear enhancements as can be seen from Eq. (2.8). This refactorization is valid up to corrections which are nonsingular at threshold when partonic z → 1. While the z space result captures the entire underlying infrared dynamics in the threshold limit, it can be better described in the Mellin-N space where the threshold limit z → 1 translates into N → ∞. We found that the form in Eq. (2.8) was already suitable for all order study, however complications arise in performing the convolution. On the other hand any such convolution becomes simple product in Mellin space and all the distributions coming from the soft function are thus translated into large logarithms in Mellin N .
Following [20], the resummed partonic SV coefficient function can be organized as follows:σ where G N is obtained by computing the large N limit of Mellin moment of G + and then by decomposing as The N independent constant g 0 is given by G N is function of the universal coefficients A which are known to fourth order and D known to third order in a s . G N collects and resums all the large-N logarithms to all orders and it can be expressed as a resummed perturbative series which takes the following form: Following [22,23], we computed g i up to i = 4 (for g i up to i = 3, see [30]) and they are given in Appendix B.1. Note that g i coefficients are universal in the sense that it depends only on whether the born process is qq channel or gg channel. In the Mellin N space, the δ(1 − z) in z-space directly translates into N independent piece whereas the plus-distributions give rise to the ln(N) as well as N independent constants in the large N limit. Part of these constant pieces, namely G 0 , is absorbed into the coefficients g 0 in the standard resummation approach. Hence, g 0 contains only N independent pieces which come from the form factor, soft distribution function, AP kernels and N independent part of the Mellin moment of G + (z, q 2 ). The condition G N = 0 for N = 1 allows the constants g i to contain N independent terms. Note that the expressions for g 0 and g i obtained this way depend on the condition G N = 0 for N = 1. In other words, there is an ambiguity in treating the N independent terms in the resummed results. Exploiting this, in [20], the N independent constants were defined by demanding G N = 1 when N = 1. With this, g 0 has the following perturbative expansion: The successive terms in the resummed exponent Eq. (2.41) along with the corresponding terms in Eq. (2.38) define the accuracies leading logarithmic (LL), next to LL (NLL), NNLL and N 3 LL etc. Terms independent of N can be treated, in principle, by the same methods that resum terms enhanced by logarithms of N .
In summary, the resummed result will differ depending on how we treat the Nindependent constants. We define various schemes that differentiate how these constants are treated in our numerical implementation for the phenomenological studies. This allows us to investigate numerical impact of the various resummed results in detail.
• Standard N exponentiation: This is the case we have discussed so far where we define large logarithms are functions of The N dependent functions G N in this case can be computed by simply performing the Mellin moment of G + (z, q 2 ) in the large N limit and keeping only those terms that vanish when N = 1.
• Standard N exponentiation: This approach differs from the previous one in the definition of large-N variable. In this case the large logarithm is simply ln N and these terms are exponentiated to all orders through the resummed exponent. It is evident that this only accounts for reshuffling of γ E between g 0 and G N in Eq. (2.29) which now takes the following form: The resummed exponent G N also takes a different form compared to the standard N exponent, The resummed coefficients g i in the above equation which defines the resummed accuracy, differs from g i in Eq. (2.41). The present scheme is defined by demanding With this definition, the rest of the N independent terms from the Mellin moment of G + is combined with finite parts of form factor, soft distribution function and the AP kernels as The N independent constant g 0 is given by and the above result is expanded in powers of a s : Numerically this can make a difference and it was seen in the context of DIS previously. In case of DY also we find such differences which will be discussed in the next section. Up to N 3 LL accuracy, the resummed exponents g i , i = 1, .., 4 for both quark as well as for gluon initiated process in N exponentiation scheme can be found in [22,28] and we computed the results for the g 0i coefficients up to i = 3 which are listed in Appendix B.2 along with g i .
• Soft exponentiation: In the standard N (N ) exponentiation, one exponentiates lnN (ln N ) and certain N (N ) independent terms which arise from G + , subjected to the condition G N = 0 (G N = 0) when N = 1 (N = 1). The remaining N (N ) independent terms in the Mellin moment of G + along with C 0 give the coefficient g 0 (g 0 ). In principle, we can define a scheme wherein entire N (N ) independent terms of G + can be kept in the exponent. More specifically, we define the scheme (relaxing G N = 1 (G N = 0) for N = 1 (N = 1)) wherein we exponentiate all the terms coming from the finite part of soft distribution function and those from the AP kernels. That is, the exponential contains The remaining N independent terms define g Soft 0 that is obtained by expanding exp(L fin F δ ) in power series expansion in a s : • All exponentiation: The soft function and the form factor satisfy K+G type Sudakov integro-differential equations given in Eqs. (2.9), and (2.18) and the AP kernels satisfy renormalisation group equation Eq. (2.15) governed by AP splitting functions. Hence, their solutions given the boundary conditions demonstrate exponential. The z space solutions that we obtained carry all order information on the distribution D j in terms of universal cusp A, soft f and collinear B anomalous dimensions and certain process dependent constants resulting from the form factor. Hence it is natural to study the numerical impact of the entire contribution in the Mellin space without imposing any condition on the N dependent terms. This can be easily achieved and the result forσ N takes the following form The present scheme was already explored in [24,25] for studying inclusive cross section for the production of Higgs boson at the LHC. For similar study for the DY in DIS and MS schemes, see [44]. Here we will extend it to the N 3 LL accuracy. The relevant resummed exponent has been provided in Appendix B.4.
Note that a detailed comparison between the N -exponentiation and N -exponentiation has been done in [29] for the charge and neutral DIS processes. There, one finds that the Nexponentiation shows a faster convergence compared to the N -exponentiation. In fact, the convergence has already been achieved at NLO+NLL order in the threshold region in the case of N -exponentiation, whereas in N -exponentiation, this occurs after the NLO+NLL order. Notice that the leading logarithmic term also differs between these two approaches. In the case of N -exponentiation, all the γ E terms are exponentiated through the variable N = N exp(γ E ); but in the N -exponentiation these γ E terms are distributed among the exponent and the N independent term g 0 . As a result the deviation starts already at the LL accuracy. In the next section, we will discuss how various schemes discussed so far can affect the predictions. Note that they all give same result at the LL accuracy, however from NLL they differ. At NNLO level, we have the contributions from all the channels and at N 3 LO only SV contribution is known so far. Hence, our numerical predictions will be based on fixed order N 3 LO sv results for the parton coefficients and on parton distribution functions known to NNLO accuracy. Note that the resummed result has to be matched to the fixed order result in order to avoid any double counting of threshold logarithms. Hence, the matched result which is usually denoted by N n LL is computed by by taking the difference between the resummed result and the same truncated up to order a n s . Hence, it contains contributions from the threshold logarithms to all orders in perturbation theory starting from a n+1 s : The Mellin space PDF (f i,N ) can be evolved using QCD-PEGASUS [61]. Alternatively they can be related to the derivative of z-space PDF as prescribed in [20,22]. The contour c in the Mellin inverse integration can be chosen according to Minimal prescription [62] procedure. Notice that the second term in Eq. (2.46) represents the resummed result truncated to N n LO order, i.e. the same order to which singular SV results are available. In the next section we present the numerical results for the DY production as well as on-shell Z, W ± production for LHC where we match the existing N 3 LO fixed order SV results with the N 3 LL resummation derived in this article.

Numerical Results
In this section, we present the numerical impact of resummed threshold corrections for neutral DY production as well as on-shell Z/W ± production at the LHC. For neutral DY production we consider all the partonic channels at the FO up to NNLO with off-shell γ * , Z intermediate states. Detailed analysis is done for 13 TeV LHC, however it can be extended to other energies as well as to other colliders.

Soft-virtual correction for neutral DY invariant mass
We start our discussion by examining the SV corrections at N 3 LO. For our numerical study, we use the following electro-weak parameters for the vector boson masses and widths, Weinberg angle (θ w ) and the fine structure constant (α): In fig. (1), we present the invariant mass distribution (left panel) of the di-lepton production for the neutral case to N 3 LO sv in QCD for 13 TeV LHC as well as the corresponding K-factors (right panel). It is worth noting here that at O(α 3 s ) level the δ(1 − z) contribution is comparable but opposite in sign to the sum of logarithmic contributions as is mentioned in [17]. The 3-loop SV corrections are found to be positive up to around Q = 400 GeV and remain negative for 400 GeV < Q < 2200 GeV and become positive thereafter as threshold logarithms dominate in the high Q region. At around 3500 GeV, the 3-loop SV corrections contribute by about 2%. The observed values of Q where this change in the sign happens are not fixed but can change with the center of mass energy of the hadrons. While the perturbation series is asymptotic and the higher orders terms are very small, the reliability of the theory predictions depends somewhat on the uncertainties due to the unphysical factorization (µ F ) and renormalization (µ R ) scales as well as those due to choice of PDFs. To this end, we estimate the 7-point scale uncertainties in the invariant mass distribution at various orders in the perturbation theory by varying the scales µ = {µ f , µ r } in the range 1 2 ≤ µ Q ≤ 2. The scale uncertainties are conveniently presented in terms of the invariant mass distribution at higher orders normalized with respect to LO ones. In fig. (2) we present these normalized distributions up to N 3 LO sv as a function of τ = Q 2 /S. At LO, there is no dependence on µ r , hence the observation that these scale uncertainties are minimum around τ = 0.001 (corresponding to about Q = 400 GeV) can be directly related to the behavior of the corresponding quark fluxes. At higher orders, the dependence on µ r and µ f is known and the scale uncertainties are found to increase with Q in the region Q > 400 GeV. For Q = 1500 GeV, they are found to be 12.55%, 6.23%, 1.50% and 1.91% respectively at LO, NLO, NNLO and N 3 LO sv . For the 3-loop SV case, the scale uncertainties are expected to get further reduced only after including the regular terms that are yet to be computed in the fixed order perturbation theory. However, as we increase Q value, even N 3 LO sv show reasonable reduction in scale uncertainty as threshold logarithms dominate over the regular terms for larger Q values. For completeness, we note that the scale uncertainties for Q = 3000 GeV are found to be 21.39%, 10.95%, 3.04% and 2.16% at LO, NLO, NNLO and N 3 LO sv respectively.  We have studied the impact of different resummation schemes as described in the previous section. First we compare the resummed results between two approaches: the Standard N and Standard N prescriptions. We find that the perturbative convergence is better in the case of N exponentiation for the scale choice µ r = µ f = Q. This can be clearly seen from fig. (3) where the convergence is already achieved at NLO+NLL whereas in N exponentiation it happens only after NLO+NLL order. At Q = 2500 GeV, we see the corrections received in Standard N exponentiation is 21.6% at NLO+NLL, 2.2% at NNLO+NNLL whereas in the Standard N exponentiation these are 6.7% and 2.3% respectively. This observation is also true for different scale choices. This is expected since naively one can expect that as we exponentiate more and more terms the convergence becomes better. In the rest of the discussion we will mention 'Standard' only in the context of N exponentiation unless otherwise stated.

Resummed prediction for neutral DY invariant mass
We now investigate the differences resulting from two approaches viz. the Soft exponentiation and All exponentiation to study their perturbative behavior. To illustrate this, we show fig. (4) where we took the ratio with respect to the Standard N results at each order. Notice that LO+LL results are same for all these three approaches by construction. To this end one sees that at lower orders the resummed cross-sections are improved over N exponentiations. At NNLL the Soft exponentiation gets additional 0.12% corrections compared to the Standard N approach at Q = 100 GeV. However at N 3 LL level the Soft exponentiation does not improve over the StandardN results and both approaches provide almost same results. On the other hand, All exponentiation still gets some contribution from higher orders through the exponentiation of complete g 0 even at N 3 LL order. The increment is however very small giving only 0.12% corrections over the StandardN scenario.
We have quantified the impact of resummed results through K-factor. In fig. (5) we present the resummed K-factors (K N LO+N LL , K N N LO+N N LL , K N N LO+N 3LL ) up to order N 3 LL. We define the K-factor as dσ resum dQ / dσ LO dQ , where resum represents all the    resummed corrections up to NNLO+N 3 LL. One observes that the perturbative convergence is improved in the case of All exponentiation compared to others although marginally. The K factor defined this way will be useful to directly compare against the experimental results. For All exponentiation case, we find that the K-factor is 1.294 at Q = 100 at NNLL which changes to 1.286 at N 3 LL. The K-factor increases with Q. At higher Q = 2500 GeV the K-factors become 1.362 at NNLL and 1.350 at N 3 LL. Next we study the uncertainties resulting from unphysical scale in these approaches. We follow the canonical variation of µ f and µ r around the final state invariant mass Q within [1/2, 2]Q imposing additional constraint 1/2 ≤ µ r /µ f ≤ 2 as was done in the third order SV prediction in the previous section. We notice that different approaches for resummation provide a systematic scale reduction at lower invariant mass of the di-lepton pair. For example, in the Standard N case, the scale uncertainty reduces from 13.37% at NLO+NLL to 1.99% at NNLO+NNLL and 0.56% at N 3 LO sv + N 3 LL. Similar pattern is seen for the Soft and All exponentiation as well as seen in fig. (6). However, when we compare among themselves, we see that in the case of All exponentiation the scale uncertainty is reduced to 1.65% at NNLO+NNLL compared to 1.99% forN exponentiation and 2.09% for Soft exponentiation at the same order. At the N 3 LO sv + N 3 LL, however All exponentiation gives relatively larger scale uncertainty compared to the other two approaches. At some high invariant mass (say Q = 2500 GeV), we see a better scale estimate at order NNLO+NNLL where we observe that the scale uncertainty systematically improved from N exponentiation from 0.53% to 0.51% for Soft exponentiation and 0.43% for All exponentiation. However at N 3 LO sv + N 3 LL order we observe an over-estimation of scale uncertainty which gets larger for different approaches and can reach the size of NLO scale uncertainty. This shows that the sub-leading regular pieces are also important to capture the scale dependence properly. This behavior is unlike the Higgs case where one sees a certain scale improvement for exponentiation of complete g 0 . We will again come back on this discussion at the end of this section. We have also estimated in our resummed predictions the uncertainties from the nonperturbative PDFs. We convolute the resummed coefficient at N 3 LL level with n different sets of a given PDF group and estimate the uncertainty from the lhapdf routines. We use the PDFs provided by ABMP16 (n= 30) [65] , CT14 (n=57) [66], MMHT2014 (n=51) [64], NNPDF31 (n=101) [67] and PDF4LHC15 (n=31) [68] groups. These results are shown in fig. (7) in terms of δσ/σ where δσ is the difference between the extrema obtained from n different sets and σ is the one obtained from central set n = 0. These PDF uncertainties in general are found to increase with the invariant mass of the di-lepton pair and, for the range of Q considered here,we find that they are smallest in the low Q-region for AMP16 and are largest for CT14 case. These uncertainties for Q = 1500 GeV are found to be 6.14% (AMBP16), 16.99% (CT14), 6.17% (MMHT2014), 4.21% (NNPDF31) and 7.43% (PDF4LHC15).
Finally, we discuss the matching relation presented in eq. (2.46). We notice that the matching relation eq. (2.46) can be interpreted in two ways. One can match the N 3 LO sv fixed order results (with n = 3) with the resummed results subtracted up to O(a 3 s ) (with n = 3) in order to avoid any double counting from the fixed order. So far, we have followed this approach. Instead we can match the complete NNLO fixed order result with the resummed result subtracted up to O(a 2 s ), which also avoids double counting and retains the  threshold terms at O(a 3 s ) in N -space in the threshold limit N → ∞. The difference in these two approaches is sub-leading and has to be related with the fact that N -space threshold results when transformed back into distribution space produces sub-leading logarithms in addition to the plus distributions. In fig. (8) we compare these two approaches setting all the scales same as Q in Standard N approach. We see that the threshold terms defined in Mellin-N space provide much better perturbative convergence compared to the z-space definition. This is a well-known observation which shows that the sub-leading pieces are also important at this order. As far as scale uncertainty is concerned, this approach gives better estimate of scale uncertainty at N 3 LL level reducing in some cases by a factor of two, however the general behavior does not change much.

Resummed prediction for Z/W ± productions
In this section we present the resummed results for on-shell Z and W ± productions to N 3 LO sv +N 3 LL accuracy. We use 13 TeV as centre of mass energy at the LHC. We set all the parameters same as the previous section. For pdf, we chose the central value from MMHT2014 set at the corresponding order. At the LHC, the underlying parton fluxes for W + production are larger than for W − case, consequently the production cross sections for the former case are larger than the latter one. This is true also for higher centre of mass energies. In tab. (1), (2), (3), we present the central predictions for on-shell Z, W + and W − respectively with the corresponding percentage scale uncertainties. Note that the scale uncertainties are calculated again using the same procedure i.e. the sevenpoint scale variation around the central scale which is now the vector boson mass i.e. the central scale has been chosen as (µ r , µ f ) = (1, 1)M V , with V = Z for Z production and V = W ± for W -boson production. In all the cases we observed that the fixed order scale  uncertainties are systematically reduced while going to higher orders, however at N 3 LO sv , it again increases which is due to the fact that, at this order there is still missing pieces which are essential to the scale uncertainty. Similar observation is also seen to the matched resummed prediction. We see that compared to the fixed order, the resummed results provide better perturbative convergence. The scale uncertainty is also seen to improve starting from NNLO level compared to fixed order. The resummed K-factors as defined before, however increases from NNLO+NNLL to N 3 LO sv +N 3 LL for all the cases. The absolute size of the perturbative corrections however decreases at N 3 LO sv +N 3 LL compared to the previous orders confirming the reliability of perturbation theory.

Conclusions
We have studied the Drell-Yan production as well as on-shell Z, W ± productions in the context of threshold resummation. We have used all the necessary ingredients available to perform this task, in particular the threshold enhanced large-N as well as the N -independent constants. The standard threshold resummation heavily reuses the results of the SV crosssection at a particular order. In particular we showed how the the large N -independent constants can be found at N 3 LL level using the existing SV results. We also explore other possibilities of doing resummation where we exponentiate the complete soft pieces coming from the soft distribution function and also exponentiate the complete g 0 coefficients including the form factor. All these different approaches show a systematic behavior of the resummed perturbative series which gets better when more and more terms are being exponentiated in terms of perturbative convergence. We have matched our resummed N 3 LL results with the existing NNLO(N 3 LO sv ) cross-section and presented results for 13 TeV LHC. We observe a systematic decrease of the size of the corrections at the third order. At this accuracy however the missing regular pieces are also important to tame the scale un-certainty. The results for inclusive DY and Z, W ± production demonstrate the ambiguity on exponentiation of N -independent terms in the resummed results. A Soft-Virtual coefficient in N -space The SV coefficient up to three loops are presented here (denoting L = lnN), The coefficients g 0i are given in Eq. (B.1).

B Resummed coeficients
Here we collect N -dependent and N -independent coefficients for all different prescriptions for resummation.

B.1 Resummation ingredients for the Standard N exponentiation
For the standard N exponentiation we present here the N independent coefficients g 0 to three loops in Eq. (2.41) below g 01 = G 11 2 + G 11 2 + B 1 2 L qr − 2 L f r + A 1 5 ζ 2 , (B.1) g 02 = G 21 1 +G 12 2 β 0 +G 11 − 2 β 0 L qr +G 2 11 2 + G 21 1 The resummed exponent as in Eq. (2.30) is calculated to the N 3 LL accuracy and collected below: All the anomalous dimensions and constants can be found in Appendix ??

B.2 Resummation ingredients for the Standard N exponentiation
Below we present the resummed exponent for the Standard N -exponent as given in Eq. (2.40).
constants g 0 are given by

B.3 Resummation ingredients for the Soft exponentiation
In the case for 'Soft exponentiation', all the terms coming from the soft function are exponentiated and hence this means all the contribution to the finite (N-independent) piece from the soft function is also being exponentiated. This renders the g 0 coefficients of the Standard N threshold and changes also the resumed exponent. We write these changes below in terms of the Standard N threshold exponent and pre-factor, (B.14) The N -independent constants in the case can be put in the following form: where the coefficients ∆ Soft g 0i are given by,

B.4 Resummation ingredients for the All exponentiation
In the case for 'All exponentiation', the complete g 0 is being exponentiated along with the large-N pieces. This brings into modification only for the resummed exponent compared to the 'Standard N exponentiation'. We write the resummed exponent in this case in terms of N exponents as, g All 1 = g 1 , g All 2 = g 2 + a s ∆ All g 2 , g All 3 = g 3 + a 2 s ∆ All g 3 , g All 4 = g 4 + a 3 s ∆ All g 4 , (B. 19) where ∆ All g i terms are found from exponentiating also the complete g 0 prefactor and they are given as, ∆ All g 2 = g 01 , The universal D coefficients are given as,