Next-to-soft-virtual resummed rapidity distribution for Drell-Yan process to $\rm{\textbf{NNLO}+\overline{\textbf{NNLL}}}$

We present the differential predictions for the rapidity distribution of a pair of leptons through the Drell-Yan (DY) process at the LHC taking into account the soft-virtual (SV) as well as next-to-soft virtual (NSV) resummation effects in QCD perturbation theory to next-to-next-to-leading-order plus next-to-next-to-leading-logarithmic ($\rm{{NNLO+\overline{NNLL}}}$) accuracy. We perform the resummation in two dimensional Mellin space using our recent formalism \cite{Ajjath:2020lwb} by limiting ourselves to contributions only from quark anti-quark ($q \bar q$) initiated channels. The resummed corrections to the fixed order results are computed through a matched formula using the minimal prescription procedure. We find that the resummation at next-to-leading-logarithmic (next-to-next-to-leading-logarithmic) level brings about 3.98\% (1.24\%) corrections respectively to the NLO (NNLO) results at the central scale value of $q=M_Z$ for 13 TeV LHC. We also observe that the sensitivity to the renormalisation scale gets improved substantially by the inclusion of NSV resummed predictions at $\rm \overline{NNLL}$ accuracy. Further, the lack of quark gluon ($qg$) initiated contributions to NSV part in the $\rm \overline{NNLL}$ resummed predictions leaves large factorisation scale dependence indicating their importance at NSV level as we go to higher orders in perturbation theory.


I. INTRODUCTION
The production of a pair of leptons, known as the Drell-Yan (DY) production [2], is one of the well studied processes at TeV colliders such as Tevatron and the Large Hadron Collider (LHC). This was possible due to wealth of precise theoretical predictions both in Standard Model (SM) and beyond SM taking into account corrections from various sources. Being a least contaminated process, DY production is used as luminosity monitor [3] at the LHC. In addition, precise data on the rates allows one to fit the parton distribution functions of hadrons [4][5][6]. Any deviation from the precise predictions can be used to set bounds on the parameters of models of new physics.
While the DY process at leading order (LO) is purely electroweak, the radiative corrections are dominated by Quantum Chromodynamics (QCD) and it is an active area of interest for several decades, see [7,8] for first nextto-leading order (NLO) results in perturbative QCD for DY process, and for invariant mass distribution of a pair of leptons at NNLO see [9,10]. For the same observable at N 3 LO level, the dominant soft-virtual (SV) contributions were obtained in [11,12] prior to complete result [13] at N 3 LO become available recently. Electroweak correction beyond LO can be found in [14,15]. In addition to invariant mass distribution, other differential distributions namely rapidity, transverse momentum are known to N 3 LO in QCD, see [16][17][18][19][20][21][22][23][24][25]. For, mixed QCD and EW corrections see [26] and for parton showers matched with NLO QCD results, see MC@NLO [27], POWHEG [27,28] and aMC@NLO [29] frameworks.
The fixed order predictions are improved by resumming large threshold logarithms resulting from soft glu-ons, see [30][31][32][33][34][35][36]. For the transverse momentum distribution, at small p T , the resulting large logarithms exponentiate in the impact parameter space [37,38]. In the soft-collinear effective theory (SCET) one performs resummation in momentum space, see [39] for inclusive production and [40] for transverse momentum distribution. In the Mellin space, resummation of large logarithms of the Feynman variable x F which describes the longitudinal momentum of the final state was achieved in [31] and it was found that there were two thresholds that could be resummed to all orders, also see [41] for a different scheme. In [42] the resummation of rapidity of W ± in the Mellin-Fourier (M-F) space was performed following a conjecture (see [43]) and later on it was applied for Drell-Yan production in [44,45]. A similar approach in SCET can be found in [46,47].
Following [21,22,31], in [48,49], we studied the soft gluon resummation for the rapidity distributions of Higgs boson and also of a pair of leptons produced in hadron colliders. In the threshold limit, i.e., when the scaling variables z 1 → 1 and z 2 → 1, the soft gluons contribute through delta functions and plus distributions in the partonic cross sections. These contributions can be resummed to all orders both in z 1 , z 2 space and in N 1 , N 2 space. The resummed results known to desired logarithmic accuracy can be used to predict certain highest logarithms in the fixed order, see [22,24,50]. The threshold limit denoted by (z 1 → 1 , z 2 → 1) corresponds to (N 1 → ∞, N 2 → ∞) in the Mellin-Mellin (M-M) space, giving large logarithms of the form ln n (N i ), where n = 1, · · · and i = 1, 2 and the resummation in M-M space resums terms of the form ω = a s β 0 ln(N 1 N 2 ) through a process independent function g(ω) and a process dependent but N i independent function g 0 . Here, arXiv:2112.14094v3 [hep-ph] 6 Jun 2023 a s = g 2 s (µ 2 R )/16π 2 with g s being the strong coupling constant and µ R , the renormalisation scale. The constant β 0 is the leading coefficient of the beta function in QCD.
Contrary to naive expectation, in certain inclusive [13,[51][52][53] and differential [54] observables, one finds that the contributions from subleading threshold logarithms, called next to soft-virtual (NSV) terms contribute significantly at every order in perturbation theory. They are found to be present in most of the partonic channels unlike the leading logarithms. Thanks to the availability of these fixed order results to unprecedented accuracy, there are enormous developments in the understanding of these subleading terms. In particular, questions related to their structure to all orders are still open, see, [55][56][57][58][59][60][61][62][63][64][65][66][67] for more details. Recently, in a series of articles [68][69][70][71], we studied a variety of inclusive observables to understand these subleading logarithms. Remarkably, we found that there exists a systematic way to sum them up to all orders in z as well as in the Mellin N spaces, exactly the way one sums up leading threshold ones. This was achieved only for the diagonal channels. One finds that resummation of both SV and NSV terms can be achieved N space. Later on, this was extended to study NSV terms present in rapidity distributions [1] of a pair of leptons in DY and a Higgs boson in gluon fusion as well as in bottom quark annihilation. For a generic case of n-colorless particles in the final state, see [72] for. Like the inclusive one, these subleading logarithms can be resummed to all order in multi-dimensional space (spanned by z l or N l ) along with the leading threshold logarithms In the present paper, we discuss the phenomenological aspects of resummed NSV terms for the production of a pair of leptons at the LHC. In the subsequent sections, we briefly recapitulate the relevant theoretical results followed by a detailed study on the numerical impact of NSV contributions and finally we conclude our findings.

II. THEORETICAL OVERVIEW
In the QCD improved parton model, the double differential distribution of a pair of leptons in DY process with respect to their invariant mass q 2 and rapidity y can be written as , where σ q B (x 0 1 , x 0 2 , q 2 ) is the Born cross section, The dimensionless scaling variable τ is given by where q being the momentum of the pair of leptons and S = (p 1 + p 2 ) 2 in center of mass energy of incoming hadron with momenta p i respectively. The rapidity y of the lepton pair is given by y = 1 2 ln p2.q p1.q z2 , µ 2 F where x 0 1 /z 1 , x 0 2 /z 2 are their momentum fractions respectively and are renormalised at the factorisation scale µ F . ∆ q d,ab (a s , z 1 , z 2 , q 2 , µ 2 F ) are the Drell-Yan coefficient functions (CFs) obtained from the partonic subprocesses after mass factorisation at the scale µ F and are calculable order by order in QCD perturbation theory in powers of a s . These CFs beyond leading order in perturbation theory contain distributions such as δ(1 − z i ) and ln m−1 (1−zi) with m ≤ 2n, n being the order of perturbation and regular functions of z i . Distributions show up only in the diagonal CFs, called SV terms and are denoted by ∆ SV d,q while the regular part, also called the hard part is given by ∆ q,H d,ab . The leading terms in the hard part in the threshold expansion are nothing but the NSV terms. Unlike the SV terms, the NSV terms get contributions from both diagonal as well as non-diagonal channels. Each term belonging to NSV contribution contains a pair of either and regular term log k (1 − z j ), j = 1, 2; k ≥ 0. Following the work by Catani and Trentadue [31], in [21,22,24], the resummation of SV terms for the rapidity distribution to all orders was achieved in the scaling variables z i in z 1 , z 2 space and later extended it to N 1 , N 2 space in [48,49] by performing two dimensional Mellin transformations in the large N 1 , N 2 limit to obtain resummed result in the M-M space. In the Mellin N l space, when N l are large, the NSV terms take the form ln k N j /N l with (j, l = 1, 2), (k = 0, 1 · · · ) up to 1/N α 1 /N β 2 , α, β ≥ 1. In [1], restricting to diagonal channels, they were systematically resummed to all orders along with SV terms.
The diagonal CF, taking into account both SV and NSV, denoted by ∆ SV+NSV d,q was shown to exponentiate in [1] as where the symbol "C" refers to convolution which acts on any exponential of a function f (z) takes the following expansion: Here, we keep only SV distributions, namely, δ(z l )δ(z j ), δ(z l )D i (z j ), D i (z l )D k (z j ), and NSV terms D i (z l ) ln k (z j ) and δ(z l ) ln k (z j ) with (l, j = 1, 2) (i, k = 0, 1, · · · ) resulting from the convolutions. In (2), ϵ is the complex valued parameter in the dimensional regularisation scheme. The function Ψ q d in the above equation has the following integral representation in z-space wherez l = 1 − z l , q 2 l = q 2 (1 − z l ), q 2 12 = q 2 z 1 z 2 and the subscript + indicates standard plus distribution.
In (4), P q (a s , z l ) = P q (a s , z l ) − 2B q (a s )δ(z l ) with P q (a s , z l ) being the splitting function in QCD which takes the following form, with A q and B q being the cusp and collinear anomalous dimensions, L q (a s , z l ) ≡ C q (a s ) ln(z l )+D q (a s ). The cusp (A q ), the collinear (B q ) anomalous dimensions and the constants C q and D q to third order are available in [73][74][75]. The anomalous dimensions A q , B q , C q and D q can also be found in the appendix of [71]. In (5), we drop those terms which contribute to non-diagonal NSV and beyond NSV terms throughout. The function Q q d in (4) is given as In the above equation, the SV coefficient D q d are known to third order [48] in QCD and the NSV coefficient φ f d,q is parametrized in the following way, The upper limit on the sum over k is controlled by the dimensionally regularised Feynman integrals that contribute to order a i s . The constant g q d,0 in (4) results from finite part of the virtual contributions and pure δ(z l ) terms of real emission contributions. The NSV coefficients φ q,(k) d,i in (7) are known to third order [1] and are listed below: where the constants C A = N c and C F = (N 2 c − 1)/2N c are the Casimirs of SU(N c ) gauge group, n f is the number of active flavours and ζ i are the Riemann zeta functions.
In [1], we systematically computed the analytic expression of resummed CFs in the M-M space after performing the double Mellin transformation on ∆ q d in (2) as Here, the resummed result for Ψ q d,N1,N2 takes the following form: where Here ω = a s β 0 ln N 1 N 2 , ω l = a s β 0 ln N l for l = 1, 2 and h q d,01 (ω) = 0. The expressions in (10) and (11) are slightly different from those in Eq. (12) and Eq.(13) of [1], due to the explicit w l dependence in the diagonal termsh q d,ii , which needs an explanation. We notice that the diagonal terms h q d,ii for i ≥ 2 (or in general h c d,ii for c=q,g,b which represent the Drell-Yan process, Higgs production in gluon fusion and bottom quark annihilation respectively) involve only the previous order information and hence can be included in the earlier order. Taking this fact into account, we redefine the diagonal terms at each order in the following way: The results of h q d,ii remain the same as in [1]. Due to this rearrangement, the expression in (10) is better predictive as compared to the Eq.(12) of [1]. This is clear when we compare Table II in this paper with Table 1 of [1], where we depict the predictions for the NSV logarithms at higher order logarithmic accuracy using the lower order information. We emphasize that, both the definitions correctly predict the higher order resummed terms at every logarithmic accuracy and the only difference comes in how much lower order information are required for higher order resummed predictions.
The SV resummation coefficients, which comprises of g q d,0 and g q d,i are extensively discussed in [48,76,77]. The NSV coefficients g q d,i , h q d,ik andh q d,ii are listed in the appendices A and B. Our next task is to include these resummed contributions consistently in the fixed order predictions to understand the phenomenological relevance of the NSV resummed results in the context of rapidity distribution for the lepton pair production in Drell-Yan process.

III. PHENOMENOLOGY
In this section, we study the impact of resummed soft virtual plus next-to-soft virtual (SV+NSV) results for the rapidity distribution of a pair of leptons produced through Z and γ * in the collision of two hadrons at the LHC with centre of mass energy 13 TeV. We take n f = 5 flavors, the MMHT2014(68cl) PDF set [78] and the corresponding a s (M Z ) through the Les Houches Accord PDF (LHAPDF) interface [79] interface at each order in perturbation theory. For the fixed order rapidity distribution, we use the publicly available code Vrap-0.9 [17,18]. The resummed contribution is obtained from ∆ q d,N1,N2 in (9) after performing Mellin inversions which are done using an in-house FORTRAN based code. The resummed results are matched to the fixed order result in order to avoid any double counting of threshold logarithms. The matched result is given below in (14). Here e q is the charge of electron. The numerical values for the various parameters are taken from the Particle Data Group 2020 [80] and are listed below : where σ q B is given by The first term in (14), d 2 σ q,N n LO /dq 2 dy , corresponds to contributions resulting from fixed order results up to N n LO. The second term on the other hand contains only SV and NSV logarithms but to all orders in perturbation theory. The subscript "tr" in ∆ q d,N1,N2 indicates that it is truncated at the same order as the fixed order after expanding in powers of a s . Hence, at a given order a n s , the non-zero contribution from the second term starts at order a n+1 s and includes SV and NSV terms from higher orders. The Mellin space PDF (f i,N ) can be evolved using QCD-PEGASUS [81]. Alternatively, we use the technique described in [31,82] to directly deal with PDFs in the z space. The contour c i in the Mellin inversion in (14) can be chosen according to Minimal prescription [83] procedure. Note that the same Minimal prescription procedure for the SV terms [83] will go through for the Mellin inversion of NSV terms as well. This is because the inclusion of NSV logarithms does not introduce any new Landau pole. Further, the Landau pole problem is directly related to the fact that the coupling a s enters into the non-perturbative region which results due to the integration of the argument of running coupling. And the position of the Landau pole is decided by the argument of a s in the integral representation (4) which is used to resum the large SV and NSV logarithms. Since the inclusion of NSV terms does not alter the argument of a s in (4) from the SV case, the Landau pole will remain same as that of SV resummation. In addition, as already noted in [83], there will not be any power corrections induced for the SV resummation and this observation continues to be valid for the NSV case as well due to the Minimal prescription formula in (14).
Our numerical results for the fixed order are based on NNLO computation of the CFs in which parton distribution functions are taken up to NNLO accuracy. Although the results of fixed order rapidity distribution at N 3 LO are presented in [25] for the Drell-Yan process, the corresponding numerical code is not publicly available, therefore we could not include N 3 LO results in our analysis. The resummed SV and NSV results are computed up to NNLL accuracy. To go beyond the NNLL accuracy, we need the collinear anomalous dimensions C q and D q to fourth order as well as the four loop NSV coefficient φ q,4 d,4 which are currently not available. To distinguish between SV and SV+NSV resummation, all along the paper, we denote the former by N n LL and the latter by N n LL for the n th level logarithmic accuracy.
In Tables I and II, we list the resummed exponents which are required to predict the tower of SV and NSV logarithms respectively in ∆ q d,N1,N2 at a given logarithmic accuracy. Let us first discuss the predictions for the SV logarithms. As it is shown in Table I, using the first set of resummed exponents {g q d.0,0 , g q d,1 } which constitute to the SV-LL resummation, we get to predict the leading SV logarithms of ∆ q,(i) d,N1,N2 to all orders in perturbation theory, i.e, ln l N 1 ln k N 2 with l + k = 2i (l, k ≥ 0) at the order a i s for all i > 1. Further, using the second set of resummed exponents {g q d,0,1 , g q d,2 } along with the first set, one can predict extra the next-to-leading SV logarithms i.e the towers {ln l N 1 ln k N 2 } with l +k = 2i−1, 2i−2 for ∆ q,(i) d,N1,N2 with i > 2. These towers of logarithms belong to the SV-NLL resummation. In general, using the n-th set {g q d,0,n , g q d,n+1 } in addition to the previous sets, we can predict the highest (2n + 1) towers of SV logarithms in N l with l = 1, 2, which constitutes to the SV-N n LL resummation, at every order in a i s for all i > n + 1, where n = 0, 1, 2 · · · .
Next, we discuss the predictions for the NSV logarithms present in ∆ q d,N1,N2 at a given logarithmic accuracy. As it is presented in Table II, using the first set of resummed exponents {g q d,0,0 , g q d,1 , g q d,1 , h q d,0 } which constitute to the LL resummation, we can predict the highest NSV logarithms of ∆ q,(i) d,N1,N2 to all orders in perturbation theory, i.e, { ln l N1 N1 ln k N 2 , ln l N2 N2 ln k N 1 } with l+k = 2i−1 at the order a i s for all i > 1. Similarly using the second set of resummation exponents {g q d,0,1 , g q d,2 , g q d,2 , h q d,1 } along with the first set, one can predict the next-to-highest NSV logarithms to all orders, which includes the tow- These towers of logarithms contribute to the NLL resummation. In general, using the n-th set {g q d,0,n , g q d,n+1 , g q d,n+1 , h q d,n } in addition to the previous sets, we get to predict the highest (n + 1) towers GIVEN

PREDICTIONS -SV Logarithms Logarithmic Accuracy
Resummed exponents The set of resummed exponents g q d,0,n , g q d,n which is required to predict the tower of SV logarithms in ∆ q,(n) d,N1,N2 at a given logarithmic accuracy in the Mellin N -space. Here, i, j ≥ 0 and L i l = ln i N l with l = 1, 2.

GIVEN PREDICTIONS -NSV Logarithms Logarithmic Accuracy
Resummed exponents of NSV logarithms in N l with l = 1, 2, which constitute to the N n LL resummation, at every order in a i s for all i > n + 1.
Below we present the resummed result given in (2) at various logarithmic accuracy and discuss the result-ing predictions for the NSV logarithms in (N 1 , N 2 )-space as displayed in Table II till a 4 s (N 4 LO). Note that we set µ 2 R = µ 2 F = q 2 in the expressions of the predictions throughout. We begin with the resummed result in the LL approximation given by Now we expand the above expression up to a 4 s (N 4 LO) and compare the predictions for leading NSV logarithms against those from fixed-order results. Note that, as can be seen from Table II, the LL resummation which comprises of only one loop anomalous dimensions and SV+NSV coefficients from fixed-order NLO results, predicts the leading loga- , a 4 s (N 4 LO) and so on respectively as shown below. At a 2 s (NNLO), we have The prediction at a 4 Our predictions for the leading NSV logarithms are compared against the fixed-order results in [9,10,84] up to third order.
Let us now turn to the resummed result at NLL accuracy which reads as Note that at NLL accuracy, we require anomalous dimensions up to two loops and second order SV+NSV coefficients obtained from NNLO results. After expanding the above result up to a 4 s , we obtain the predictions for the next-to leading NSV logarithms etc at a 3 s , a 4 s and so on respectively and are given by where γ E is the Euler-Mascheroni constant. The above predictions for the next-to leading NSV logarithms are in agreement with results given in [9,10,84] up to third order. Furthermore, we have compared our full third order results for ∆ q d (z 1 , z 2 ) in z space (see [1] for the third order results in z space) with the results obtained using the generalized threshold factorization approach presented in [84]. From that comparison, we have found that our third order prediction is in complete agreement with results given in [84] for terms of the type D i (z l ) ln j (z m ), i, j ≥ 0, l, m = 1, 2 in the z space. However, we could not compare the remaining δ(z l ) ln j (z m ) terms in our result because they were not available in [84].
Finally, using the NNLL resummation, which further embeds the three loop anomalous dimensions and third order SV+NSV coefficients obtained from N 3 LO results, we predict the next-to-next-to leading logarithms at a 5 s (N 5 LO) and so on. The resummed expression at NNLL accuracy is given by The prediction for the next-to-next-to leading NSV logarithm at a 4 s is provided below The above predictions for the NSV logarithms in rapidity distribution ∆ q,(i) d,N1,N2 are found to reproduce the corresponding predictions in the inclusive cross section ∆ q,(i) N computed in [71] in the limit N 1 = N 2 = N for i ≤4.
In [71], it has been shown, in the context of DY inclusive cross section up to N 3 LO in the Mellin N space that although the total contribution of NSV logarithms is smaller as compared to the SV counterpart, it is still numerically sizeable to not to be neglected. Despite being formally subleading, the considerable contribution of the NSV terms is due to their large coefficients and this trend was also observed in the case of Higgs boson production through gloun fusion [70,85]. Here, even in the case of rapidity distribution, we expect the same trend to be observed. Keeping this in mind, we ask the following questions that can shed light on the relevance of NSV terms in the context of rapidity distribution in di-lepton pair production in DY process at the LHC.
• In comparison to the fixed order corrections, how much is the effect of SV+NSV resummed results on the rapidity distribution?
• How sensitive is the SV+NSV resummed rapidity distribution to the choices of factorisation (µ F ) and renormalisation (µ R ) scales?
• How do the resummed NSV terms alter the predictions of SV resummed result?
In the following sections, we address the above questions in detail. Let us begin with analysing the impact of SV+NSV resummed results in comparison to the fixed order results, which is the topic of the next section. We present our results for the doubly-differential distribution with respect to invariant mass q and rapidity y by plotting it as a function of y for fixed values of q.

A. Fixed order vs Resummed rapidity distribution
In this section, we study the effects of SV+NSV resummation on the fixed order predictions for rapidity distribution of di-lepton pair production in DY process for 13 TeV LHC. Through (14), we get the resummed predictions at LL, NLL and NNLL matched with the corresponding fixed order results. The numerical impact of higher order contributions can be quantified through the K-factors defined below, where we have set renormalisation (µ R ) and factorisation (µ F ) scales at q. Below, in Table (III), we show the K-factor values of both fixed order and resummed results for benchmark rapidity values at the central scale µ R = µ F = M Z . We observe that there is an enhancement of 32.9% and 36.9% when we go from LO to NLO and NNLO respectively, at the central rapidity region. Furthermore, the fixed order values at LO, NLO and NNLO get incremented by the inclusion of SV + NSV resummed predictions at LL, NLL and NNLL by 4.9%, 3.98% and 1.24% respectively, around the central rapidity region. This can be seen from the right panel of fig.(1), where the resummed curves are found to lie above their respective fixed order ones implying the enhancement resulting from the resummed corrections. It should be noted that the K-factor curves of the resummed results at NLO + NLL and NNLO + NNLL overlap for a wide range of rapidity values which was not observed for the case of fixed order predictions. This indicates that the perturbative convergence is improved among the resummed results, thereby leading to the reliability of perturbative predictions by the inclusion of resummed corrections. We also notice that the K-factor values are closer for NNLO and NNLO + NNLL as compared to NLO and NLO + NLL over the full rapidity region. This suggests that the resummed contributions to the fixed order rapidity distribution decrease as we go to higher orders in perturbation theory.
From the above analysis of K-factors, we have observed  that the resummed predictions not only bring in considerable enhancement in the fixed order results, but also improve the perturbative convergence till NNLO+NNLL accuracy. However, both fixed order and resummed predictions suffer from the presence of unphysical scales, namely, the renormalisation µ R and the factorisation µ F scales. Therefore, a careful study of perturbative uncertainties of these predictions is needed by studying their sensitivity to the choices of µ F and µ R scales, which will be discussed in the following subsection.

7-point scale uncertainties of the resummed results
The uncertainty associated with the choice of renormalisation µ R and the factorisation µ F scales due to the truncation of the perturbative series can be assessed using the standard canonical 7-point variation, where µ = {µ F , µ R } is varied in the range 1 2 ≤ µ q ≤ 2, keeping the ratio µ R /µ F not larger than 2 and smaller than 1/2. In fig.(2), we compare the 7-point scale uncertainties of the SV + NSV resummed results (right panel) against fixed order ones (left panel) around the central scale choice (µ R , µ F ) = (M Z , M Z ) for 13 Tev LHC at various perturbative orders. Here, we find that the central scale lines of resummed predictions are shifted up with respect to that of corresponding fixed order results. This indeed suggests that there is a systematic enhancement in the rapidity distribution when we add the resummed corrections to the fixed order results as shown in Table (IV). This was also observed from the analysis of K-factor values given earlier. However, we notice that the uncertainty bands of the resummed predictions are wider than that of the corresponding fixed order ones over the entire rapidity range at every order of perturbation. Nevertheless, the uncertainty band decreases as we go to higher logarithmic accuracy from LO + LL to NNLO + NNLL. In addition, the error band of NNLO + NNLL is fully contained within the band of NLO + NLL over most of the rapidity region, unlike the fixed order case. In Table (IV), we present both fixed order and resummed predictions at various perturbative orders along with their asymmetric errors resulting from 7-point scale variation for benchmark rapidity values. Here, we notice an increment of 31.7% while going from LO + LL to NLO + NLL accuracy, which further improves by 0.3% at NNLO + NNLL for y = 0. Besides this, the scale uncertainty gets reduced significantly while going from LO + LL to NNLO + NNLL over the full range of rapidity.
For instance, the uncertainty ranges between (−16.19%, +15.36%) for LO + LL, (−7.50%, +7.06%) for NLO + NLL and it is (−2.18%, +3.30%) at NNLO + NNLL for central rapidity region. However, there is no improvement in the scale uncertainty of the resummed corrections when we compare it against the fixed order counterpart. This can be explained due to following possibilities: • The resummed logarithms are not dominant in this region to show their numerical relevance.
• The lack of NSV resummed corrections from offdiagonal channels can give rise to large scale uncertainties in the resummed predictions.
Recall that the resummation is inevitable to cure the perturbative series which suffers from certain large logarithmic terms at every order, in the kinematic threshold region, where the invariant mass q approaches the hadronic centre of mass energy which is 13 TeV in our case. Therefore, in order to see the impact of the resummed contributions, we need to extend our analysis to higher invariant mass region which is of the order of TeV. We will discuss the off-diagonal channel contribution later in detail towards the end of this section. Now, we move on to the analysis of 7-point scale uncertainty of the SV+NSV resummed predictions in comparison to the fixed order results for high invariant mass q = 2 TeV.    From the earlier discussions on the 7-point scale uncertainty for q = M Z , we found that the uncertainty bands of resummed predictions were wider than that of fixed order at every order of perturbation. Nevertheless, the width of uncertainty bands were found to decrease as we moved from LO + LL to NNLO + NNLL accuracy. In addition, we also observed appreciable amount of increment in the rapidity distribution by the inclusion of SV+NSV resummed effects. Now, here in fig.(3), we show the 7-point scale variation of the rapidity distribution for q = 2 TeV. The fixed order results are depicted in the left panel up to NNLO accuracy and resummed predictions are given in the right panel up to NNLO+NNLL accuracy. In general, we note that the width of uncertainty bands corresponding to both fixed order as well as resummed predictions are significantly reduced as compared to the uncertainty bands for q = M Z . Interestingly, the NLO + NLL uncertainty band is better as compared to NLO fixed order band over the entire rapidity region. Also, the NNLO uncertainty gets improved by the inclusion of resummed NNLL corrections around the central rapidity region. This indicates the relevance of resummed contributions at this invariant mass region. This was not observed for the case of q = M Z where the resummed contributions were not prominent.
In Table (V), we quote the central scale values of both fixed order and resummed rapidity distributions at q =2 TeV along with the 7-point scale uncertainties for benchmark rapidity values. Here, we observe that the percentage uncertainties of fixed order as well as resummed results get reduced substantially at each perturbative order when we compare them against the values given in Table(IV). For instance, the uncertainty at NNLO + NNLL is reduced from (-2.18%,+3.3%) to (-0.31%,+0.53%) as we go from q = M Z to q = 2 TeV around the central rapidity region. In addition, the uncertainty at NNLO + NNLL is evidently small as compared to the uncertainty of (-0.89%,+0.60%) at NNLO for y = 0. Similarly the uncertainty at NLO is (-2.93%,+2.84%) which comes down to (-1.25%,+1.39%) at NLO + NLL for the same value of y. As for the case of q = M z , there is a systematic reduction in the uncertainties while going from LO + LL to NNLO + NNLL over the entire rapidity region which can be seen from table (V). We also find that the resummed contribution at NNLL brings in 0.86% correction to NNLO whereas it was 1.24% for the case of q = M Z . This suggests that the correction resulting from resummation at NNLL accuracy decreases as we go to higher q values leading to better reliability of resummed results.
To summarise, we found that the uncertainties of the rapidity distribution decrease by the inclusion of the resummed corrections at q = 2 TeV over the full rapidity region. Furthermore, the reliability of the perturbative results due to resummed corrections is improved at this invariant mass value. Thus, it can be inferred that the relevance of resummation effects become evidently visible while going from q = M z to q = 2 TeV. To understand these observations in a better way, we now turn to study the effect of each scale individually on the SV+NSV resummed result.
Uncertainties of the resummed results with respect to µR and µF In the following, we examine the effect of µ R and µ F scales individually on the resummed result. We begin with plotting the dependence of the rapidity distribution on µ F as a function of the rapidity y while fixing the scale µ R at the invariant mass q, for q = 2 TeV as shown in fig.(4). The bands are obtained by varying the scale µ F by a factor of two up and down around the central scale µ R = µ F = 2 TeV. Here, the resummed band depicted in the right panel at NNLO+NNLL looks similar to that of the 7-point variation band shown in fig.(3) (right panel). This indicates that the contribution to the width of NNLO+NNLL band in fig.(3) mainly comes from the uncertainties arising from variations in the µ F scale. Note that the uncertainties at NNLO+NNLL arising from µ F and 7-point variation are identical and they lie between (-0.31%, +0.53%) for y = 0. Moreover, the µ F scale uncertainties decrease as we go to higher logarithmic accuracy in the resummed results. Now, we move on to compare the µ F scale uncertainty of the resummed predictions with respect to the fixed order results. We observe that the µ F scale uncertainty of NLO gets improved by the inclusion of NLL resummed predictions whereas the NNLO band increases when the NNLL corrections are added. Let us try to understand why the SV+NSV resummed result at NNLO+NNLL is more sensitive to the µ F scale variation as compared to fixed order NNLO result. As mentioned earlier, we perform the resummation of SV distributions and NSV logarithms present in the diagonal partonic channel. Unlike the SV distributions that get contribution only from diagonal quark anti-quark (qq) initiated channel, the NSV terms can originate from off-diagonal channels like quark gluon (qg), gluon gluon (gg) etc as well. Under the µ F scale variation, these various partonic channels get mixed due to the DGLAP evolution of the PDFs. Hence, it becomes essential to keep all the contributing partonic channels at a particular perturbative order as there can be compensations among those channels, thereby reducing the scale uncertainty at that order. The fixed order results used for our numerical analysis contain all the par- tonic channels while the resummed contributions are only from qq initiated channels. Thus, the scale dependence of fixed order result is expected to go down in comparison to the corresponding resummed prediction. However, as mentioned above, the inclusion of resummed corrections at NLL accuracy improves the NLO error band. This suggests that the contribution of qg channel is not prominent at NLO. We find that the oneloop correction from the qq-channel is 23.6% while the correction from the qg-channel is only -2.5% of the NLO rapidity distribution at the central rapidity value. Therefore, there is an improvement in the µ F scale uncertainty when we sum up the collinear logarithms resulting from the dominant qq-channel at NLL. On the other hand, at NNLO level, the a 2 s corrections from qq and qg channels are 4.5% and −1.25% respectively to the NNLO rapidity distribution. As a result, the magnitude of the NNLO result is determined by a significant cancellation between qq and qg channels which was not the case for NLO. Now, due to the unavailability of the qg resummed collinear logarithms in our analysis, the aforementioned cancellation at NNLO+NNLL is not balanced. Thus, the µ F variation band of resummed prediction at NNLO+NNLL in fig.( 4) displays that the qg resummation is required to improve the results.
Next we try to understand the behaviour of resummed rapidity distribution in comparison to the fixed order counterpart under µ R scale variation. Fig.(5) shows the dependence of the rapidity distribution on µ R keeping µ F fixed at 2 TeV. The bands are obtained by varying the scale µ R by a factor of two up and down around the central scale µ R = µ F = 2 TeV. Here, the LO rapidity distribution being independent of the scale µ R , does not have a band associated with it. On the other hand, there is a band when we add the resummed corrections at LL accuracy to the LO rapidity distribution. This is because the resummed corrections at LL capture the leading logarithmic contributions from all orders in perturbation theory, thereby giving rise to µ R scale uncertainty. Moreover, the inclusion of resummed corrections at both NLL and NNLL improves the µ R scale uncertainties of NLO and NNLO respectively. This is in contrast to the case of µ F scale variation discussed earlier. Although the improvement is minuscule at NLO, it is substantial at NNLO due to NNLO+NNLL which is indeed the highlight of this plot here as compared to µ F scale variations shown in fig.(4). For instance, the µ R scale uncertainty at NNLO is reduced from (-0.56%, +0.5%) to (-0.16%, 0%) for y = 0 by the inclusion of NNLL. As we know, each partonic channel is invariant under µ R scale variation when taken to all orders. Hence, there is an improvement when we include more higher order corrections within a channel which is qq in this case, by keeping the scale µ F fixed.
In conclusion, we observed that the uncertainties due to both µ F and µ R scales decrease as we go to higher logarithmic accuracy. As far as the µ F scale variation is concerned, the resummation of collinear logarithms resulting from qg channel also plays an important role. We notice that having the qg resummed contribution is more significant at NNLO level than at NLO due to relatively larger contribution from qg channel at NNLO. As a result, the 7-point scale uncertainty of the SV+NSV resummed predictions at NNLO+NNLL is mostly driven by the µ F scale variation. Note that the inclusion of SV+NSV resummed predictions reduces the the µ R scale sensitivity remarkably at NNLO+NNLL accuracy. So far we have discussed the effects of resummation on the fixed order results taking into account SV distributions and NSV logarithms together in the analysis. Now let us turn to understand which part of the SV+NSV resummation, i.e., whether it is the resummation of the distributions or  of the NSV logarithms, plays the main role in any kind of improvement observed so far.

B. SV vs SV+NSV Resummed results
In the previous section, we have studied the effects of SV+NSV resummation on the fixed order rapidity distribution in detail. We observed that there is a considerable amount of enhancement in the rapidity distribution by the inclusion of SV+NSV resummed predictions and more importantly the µ R scale uncertainty gets re-duced substantially at NNLO+NNLL accuracy. On the other hand, the µ F scale uncertainty shows improvement at NLO+NLL for higher values of q but not at NNLO+NNLL. In the following, we perform an analysis on the inclusion of resummed NSV logarithms by comparing it with the SV resummed results.
We begin with the analysis of K-factor values for SV+NSV resummed results in comparison to the SV counterpart till NNLO+NNLL level at the central scale µ R = µ F = q for q = M Z . In Table (VI), we compare the K-factor values of SV and SV+NSV resummed predictions at various orders for benchmark rapidity values. We find that there is an increment of 3.15%, 2.75% and 0.625% in the rapidity distribution when going from LL to LL, NLL to NLL and NNLL to NNLL, respectively, at the central rapidity region. Fig.(6) demonstrates this trend for a wider range of rapidity values. In addition, the K-factors curves of NLL and NNLL almost overlap with each other for a wide range of rapidity values which is not observed for the case of NLL and NNLL curves. This suggests that there is better perturbative convergence when the NSV logarithms are taken into account.
We now turn to study the scale uncertainties arising from SV+NSV resmmation in comparison to the SV resummation. We first analyse the behaviour of both SV and SV+NSV resummed rapidity distributions as a function of y under the 7-point scale scale variation as depicted in fig.(7) for q = 2 TeV. We observe that inclusion of SV as well as SV+NSV resummed corrections reduces the uncertainty of fixed order results at both NLO and NNLO accuracy. This reductions in the uncertainty is prominent for lower rapidity values |y| ≤ 0.5 as shown in the insets in fig.(7). As can be seen from table (VII), the uncertainty at NLO+NLL is comparable to that of NLO+NLL around the central rapidity region. However, the uncertainty at NNLO+NNLL gets worse when we   add the resummed NSV contributions at that accuracy. For instance, the uncertainty at NNLO+NNLL lies in the range (-0.34%,+0.23%) whereas it is increased to (-0.31%,+0.53%) at NNLO+NNLL for y = 0. This hint towards our earlier findings in the previous section that the sensitivity of the SV+NSV resummed results to the unphysical scales increases due to the lack of resummed NSV predictions from off-diagonal qg channel. Next we move on to compare the SV and SV+NSV resummed predictions under the variation of each of these scales separately.
We first consider the behaviour of both SV and SV+NSV resummed rapidity distributions as a function of y under the µ F scale variation with the scale µ R fixed at q = 2 TeV as depicted in fig.(8). In general, the bands corresponding to SV+NSV resummed predictions are wider than that of SV predictions over the entire rapidity region. We also find that the width of the bands corresponding to fixed order rapidity distributions gets reduced with the inclusion of both SV (NLL) and SV+NSV (NLL) resummed corrections at NLO. For instance, the uncertainty is (-1.36%,+1.7%) at NLO whereas it is reduced to (-0%,0.46%) and (-0.2%,+1.07%) at NLO+NLL and NLO+NLL respectively for the cen-tral rapidity value. This can be associated with the earlier observation of qq and qg contributions at NLO. We have already seen that qq is the dominating channel at NLO and hence the uncertainty is expected to get better as we include the resummed corrections coming from that channel. On the other hand, though the uncertainty at NNLO gets improved by the addition of NNLL SV resummed corrections, it gets worse when we include the NSV corrections through NNLL. These observations can be seen from the insets in fig.(8) for lower rapidity values |y| ≤ 0.5 . As we know, the SV resummed terms come only from diagonal qq channel, therefore they do not need any compensating factor to reduce its uncertainty. In contrary to this, the NSV resummed predictions which we have included here is incomplete due to missing contributions from off-diagonal qg channel. Consequently the NSV included results will show the residual µ F uncertainty due to mixing of various partonic channels. However, the scenario will be different if we keep the scale µ F fixed and vary the renormalisation scale µ R . Fig.(9) shows the comparison of SV and SV+NSV resummed results under the µ R scale variation for q = 2 TeV. Here, for the case of µ R variation as well, the un-   certainties of fixed order results get improved by the inclusion of both SV and NSV resummed corrections. Interestingly, at NNLO, the width of the bands gets reduced substantially when the NSV resummed correction at NNLL accuracy is added in comparison to its SV counterpart. This improvement by the inclusion of NNLL NSV resummed corrections is notable for lower rapidity values |y| ≤ 0.5 as shown in the insets in fig.(9). The uncertainty at NNLO around the central rapidity region lies between (-0.57%,+0.5%) which gets reduced to (-0.33%,+0.22%) when the SV (NNLL) corrections are added. And, it gets further improved to (-0.16%,0%) with the inclusion of NSV corrections (NNLL). This emphasizes that the resummed NSV contributions play vital role in bringing down the µ R scale uncertainty as we go to higher logarithmic corrections.
In summary, we found that the uncertainty becomes better with the inclusion of both SV (NLL) and NSV (NLL) resummed corrections at NLO under µ F as well as µ R scale variations. But at NNLO, under µ F scale variation, the inclusion of NSV NNLL corrections increases the uncertainty whereas the SV NNLL corrections brings it down significantly. This indicates that the NSV resummed corrections here require the resummed contributions from qg channel as a compensating factor to improve the uncertainty. Note that in all these analysis, we studied the impact of fixed order and resummed CFs using same PDF sets to desired logarithmic accuracy for both of them. For studies related to µ F variations, it is worthwhile to consider resummed PDFs if they are available. However, as far as the µ R uncertainty is concerned, the NSV corrections show nice behaviour especially at NNLO+NNLL accuracy with notable reduction in the uncertainty. This suggests that the resummed NSV terms plays substantial role in improving the µ R scale uncertainty in comparison to its SV counterpart.

IV. DISCUSSION AND CONCLUSION
Through this article, we provide for the first time, the numerical predictions for resummed next-to soft virtual contributions up to NNLO + NNLL accuracy to the rapidity distribution of pair of leptons in the Drell-Yan process at the LHC. By restricting ourselves to the mechanism where only neutral gauge bosons like photon and Z boson produce leptons, we have used our recent formalism [1] to systematically resum NSV logarithms to all orders. In our previous work on Drell-Yan inclusive crosssection, we have quantified the significant contribution of the NSV logarithms in the fixed order predictions [71]. This serves as the motivating factor to study the numerical significance of these collinear logarithms in the case of rapidity distribution as well.
We have quantified the numerical effects of SV+NSV higher order predictions by providing the K-factor values for central scale µ R = µ F = M Z . We find that there is an enhancement of 4.9%, 3.98% and 1.24% at LO+LL, NLO+NLL and NNLO+NNLL respectively by the inclusion of SV+NSV resumed results. Also, there is an improvement in the perturbative convergence over the fixed order results till NNLO+NNLL accuracy. The sensitivity of our predictions to the unphysical scales µ R and µ F is studied using the canonical 7-point scale variation approach. We have given the plot of 7-point scale variation for two values of invariant mass, q = M Z and q = 2 TeV. We find that at q = M Z , the uncertainty of resummed predictions is more than the corresponding fixed order results till NNLO. However, at q = 2 TeV, the scale sensitivity at NLO+NLL is decreased over the entire rapidity region whereas at NNLO+NNLL, it gets reduced around the central rapidity region. Thus, by doing a comparative study of the scale uncertainties at two different q values, we infer that the resummation effects become prominent as we go to higher values of q. Nevertheless, there is a systematic reduction in the uncertainty of the resummed results while moving to higher logarithmic accuracy for both q = M z and q = 2 TeV.
Further analysis of the scale dependency revealed that the 7-point scale uncertainties of resummed predictions are largely governed by the factorisation scale µ F especially at NNLO + NNLL. Moreover, the comparative study of SV and SV+NSV resummed results shows that the NSV part of the resummation increases the uncertainty due to µ F scale variations. We know that different partonic channels mix under factorisation scale variations when they are convoluted with appropriate PDFs. Therefore, the absence of NSV contributions coming from the off-diagonal qg channel increases the sensitivity to µ F scale at the hadronic level. However, this missing compensation from qg channel is more evident at NNLO level due to considerable contribution from qg channel at this order. This suggests that the NSV resummation corresponding to qg channel is necessary to improve the predictions as we go to higher orders in perturbaton theory. In addition, as far as the µ F scale variation is concerned, resummed PDFs are also useful to be included for better results.
The independent study of renormalisation scale variation shows that the improvement in the scale uncertainty at NLO + NLL is not quantitatively significant, however at NNLO + NNLL, there is a substantial decrease in µ R scale sensitivity as compared to the corresponding fixed order results. Note that the µ R scale uncertainty at NNLO is reduced from (-0.56%, +0.5%) to (-0.16%, 0%) for the central rapidity region by the inclusion of NNLL. From the comparison of SV and SV+NSV resummed results, we find that it is the inclusion of NSV resummed corrections at NNLL accuracy to its SV counterpart which brings down the µ R scale dependency to a great extent. This is expected because different channels being renormalisation group invariant, do not mix under µ R scale variation.

V. ACKNOWLEDGEMENTS
We thank J. Michel and F. Tackmann for third order DY results of rapidity for comparing purposes and C. Duhr and B. Mistlberger for providing third order results for the inclusive reactions. Further, we are grateful to Marco Bonvini for notifying some typos in the paper. A. A. H is supported by the French ANR under the grant ANR-20-CE31-0015 ("PrecisOnium"). In addition we would also like to thank the computer administrative unit of IMSc for their help and support.
Appendix A: NSV Resummation exponents g q d,i (ω) The NSV resummation exponents g q d,i (ω) given in (10) are provided below.