Resummed inclusive cross-section in Randall-Sundrum model at NNLO+NNLL

The complete next-to-next-to leading order (NNLO) QCD correction has been studied to the di-lepton invariant mass distribution within the Randall-Sundrum (RS) framework. In addition, the soft-virtual (SV) cross-section at next-to-next-to-next-to leading order (N3LO) as well as threshold resummation to next-to-next-to leading logarithms (NNLL) level have been presented. The analytical coefficient for SV production has been obtained up to three loops very recently along with the process-dependent coefficients needed to perform resummation up to NNLL. These coefficients are universal for any universal spin-2 model where spin-2 particle couples to the Standard Model (SM) particles with equal strength. We use these coefficients in predicting N3LO SV results as well as matched NNLO+NNLL results for invariant mass distribution for Drell-Yan (DY) production in RS model. We performed a detailed phenomenological analysis and present our results in terms of mass dependent K-factors for the 13 TeV centre-of-mass energy at the Large Hadron Collider (LHC) for the search of such RS Kaluza-Klein (KK) resonances. The NNLO cross-section adds about 21% correction to the next-to-leading order (NLO) results. We found that the SV correction at the N3LO order decreases the cross-section by 0.7% near the first KK resonance (M1 = 1500 GeV) whereas the resummed result shows an increment over NNLO by 7% of LO. We performed a detailed analysis including scale variation and parton distribution function (PDF) variations. These new results provide an opportunity to stringently constrain the parameters of the model in particular in the search of heavy spin-2 resonances at the LHC.


Introduction
The Standard Model (SM) of particle physics is now well established after the discovery of scalar Higgs boson [1,2] at the Large Hadron Collider (LHC). The properties of the Higgs boson is being tested at a very high accuracy in the hope of new physics beyond the SM (BSM). A large class of the BSM scenarios are motivated by the large hierarchy between the electroweak symmetry breaking scale and the Planck scale, known as the gauge hierarchy problem. A wide class of theories have been proposed to address this problem through the introduction of large extra dimensions in the TeV scale brane world scenarios. In particular the models with warped extra dimension as proposed by Randall and Sundrum (RS) [3] are attractive candidates to solve this gauge hierarchy problem. In its simplest version, it predicts spin-2 Kaluza-Klein (KK) excitations in the TeV mass range which could be accessible at current hadron collider LHC or in any future hadron colliders or electron-positron colliders.
Search of physics beyond the SM has been an important objective of the LHC physics program. Precision physics plays an important role in this regard to accurately predict the cross-sections and distributions within perturbative framework. The process like Higgs and pseudo-scalar Higgs boson [4][5][6][7][8][9], DY [4,10] productions are already available at NNLO accuracy. The large perturbative corrections for Higgs at NNLO even pushes the accuracy to be calculated to even N 3 LO order [11][12][13]. Recently the DY production has also been calculated to third order in strong coupling [14]. The exclusive observables like rapidity are also being calculated to the same accuracy (see for example [15][16][17][18][19][20][21][22][23][24]).
In order to achieve perturbative stability, it is instructive to go beyond NNLO by computing the SV cross-section at N 3 LO order. The SV corrections constitute a significant JHEP07(2020)040 part of the full cross-section and have been successfully computed for many processes in the SM, for example, Higgs production [25][26][27][28][29][30][31][32], DY production [33,34], associated production [35] to N 3 LO as well as in BSM domain like pseudo-scalar production [36] in 2HDM. Similar accuracy has also been achieved for rapidity distributions [33,[37][38][39]. Partial four-loop results are recently available for Higgs and DY production [40].
In the threshold region where partonic z → 1, the truncated fixed order cross-section however becomes unreliable due to the presence of large logarithms. These large logarithms arise due to constrained phase space available for the soft gluons. In order to get a reliable prediction also in these corners of the phase-space, it is thus essential to resum these large logarithms to all orders. Threshold resummation has been performed successfully to inclusive Higgs production [26,[41][42][43][44][45][46], DY production [26,42], DIS [47] as well as for pseudo-scalar production [48][49][50] up to N 3 LL accuracy. The first results towards N 4 LL corrections are also available recently for DIS in [51]. Moreover for differential observables like rapidity, it is known to NNLL accuracy for many important processes (see for example [52][53][54][55][56][57]).
In the context of large extra dimension, the NLO corrections were known for many important processes at the LHC [58][59][60][61][62][63][64][65][66][67][68] within Arkani-Hamed-Dimopoulos-Dvali (ADD) [69] and RS [3] model. It is observed in the NLO QCD computation [58] that the K-factors in the di-lepton production case are potentially large and range up to 60%. The matched NLO results with parton shower is also known for di-final processes in ADD [70,71] and in RS [72] model. The associated production [73] as well as triple gauge boson production processes [74] are also known. In RS model, the triple neutral gauge boson productions are available [75] in ME+PS accuracy in the MADGRAPH framework. Moreover generic universal and non-universal spin-2 production processes are automatized [76] in FEYNRULES [77]-MADGRAPH5 AMC@NLO [78] framework providing NLO accuracy for inclusive and exclusive cross-sections for all relevant channels at the LHC.
The first attempt to go beyond the NLO accuracy has been seen in [79] calculating the SV corrections at NNLO. This has been possible due to the calculation of spin-2 form factor [80] at the same order. Shortly after, the complete NNLO corrections were computed in [81] using the method of reverse unitarity [6] and phenomenological study has been performed in the context of ADD model. There it has been found that the NNLO correction changes the cross-section by 21% over NLO results and constrains the scale uncertainty to 1.6%. Similar accuracy is also available for non-universal spin-2 production [82] where spin-2 couples with different coupling to the SM fields. The first attempt in calculating the SV corrections beyond NNLO can be seen [83] in the context of ADD model in DY invariant mass distribution after the completion of three-loop quark and gluon form factor [84]. The perturbative coefficients are same for any spin-2 production with universal coupling to the SM. There it has been noticed that the N 3 LO SV cross-section changes the NNLO by −0.7% at Q = 1500 GeV (Q being the invariant mass of the lepton pair). Moreover the authors also performed threshold resummation up to N 3 LL accuracy and the corrections are found to be around 1% over NNLO with scale uncertainty reduces to 1.5%.
In this article, we focus on massive KK production in the RS framework. We note that the search for the spin-2 resonances at hadron collider experiments has been of significant JHEP07(2020)040 interest over the last many years. The production of spin-2 particles in hadron collisons proceeds via both quark antiquark annhilation (SM DY-like) as well as the gluon fusion (Higgs-like) channels. In the case of Higgs production via gluon fusion channel, it is well known in the literature that the QCD corrections at NNLO and even beyond contribute substantially to the cross section. In the gluon fusion channel, particularly at the LHC energies, these higher order corrections are absolutely necessary for the convergence of the perturbation theory as well as for the minimization of the factorization and renormalisation scale uncertainties. It is in the same spirit such a precise QCD calculations are required for the spin-2 production cross sections as their production rates receive dominant contributions from the gluon fusion channel. These higher order corrections presented in the later sections are expected to improve by more than a few percent the bounds on the RS model parameters that are presently obtained from the LHC data using the known NLO K-factors.
Since the spin-2 RS KK excitations also couples universally to the SM stress-energy tensor, the analytical perturbative coefficients are same as the generic universal spin-2 case like ADD. Phenomenologically however the RS KK states provide very distinctive signature from that of ADD model at the LHC. Where the di-lepton invariant mass distribution in the ADD model provides a continuum distribution, in the RS model one finds well-separated massive KK resonances. Using the coefficients already obtained in ADD scenario, we first study the invariant mass distribution for DY production at NNLO accuracy in the RS model. Next we study the impact of N 3 LO SV correction as well as the NNLL resummed effect over the NNLO correction within this model.
The article is organized as follows: in section 2.1 we briefly describe the RS model and present the interaction lagrangian. In section 2.2 we set up the theoretical framework for invariant mass distribution for di-lepton production in RS model followed by the discussion on SV cross-section in section 2.3 and threshold resummation in section 2.4. In section 3, we present the distribution at NNLO and the results for N 3 LO SV as well as the resummed results matched at NNLO+NNLL accuracy. Finally we conclude in section 4.

The model
The RS background is a warped metric and can be parametrized [85] as where η µν is the flat Minkowski metric and φ is the extra dimension with periodicity 0 ≤ φ ≤ π and is compactified on a S 1 /Z 2 orbifold with radius r c . The curvature of the AdS 5 space-time is denoted as κ. In the RS model, there are two 3-branes located at two orbifold fixed points on the coordinate of the fifth dimension φ = 0 and φ = π known as the 'Planck brane' and 'TeV brane' respectively. All the SM particles are confined in the TeV brane and only the gravity is allowed to propagate into the fifth dimension. With this set-up, the hierarchy between the electroweak scale and the Planck scale is understood reasonably if the compactification radius (r c ) and the AdS curvature (κ) satisfy a condition, JHEP07(2020)040 κr c ∼ O(10) [86,87]. The higher KK modes therefore will produce observable effects on the LHC processes at the TeV range. These massive KK states (Y (n) µν ) interact with the SM fields through stress-energy tensor (T µν ) and the interaction Lagrangian is given [88,89] as, The interaction with zeroth KK mode (Y , m 0 = κe −κrcπ . The masses of the KK modes are given by M n = x n κ e −πκrc , with x n being the zeros of the Bessel function J 1 (x). The effective graviton propagator can be found after summing over all the massive KK modes except the zeroth one and it takes the form [62,90], where s ij = (p i + p j ) 2 , x = √ s ij /m 0 and Γ n denotes the width of the resonance with mass M n (see [88,91]). In the RS model, the individual KK resonances are well-separated and can be probed for example in the invariant mass distribution of lepton pairs in DY production.

Drell-Yan invariant mass distribution
The invariant mass distribution for DY production at the hadron collider is given by, Here S andŝ denote the centre-of-mass energy in the hadronic and partonic frame respectively. The mass factorized partonic coefficient function ∆ I ab is convoluted with the parton luminosity distribution L ab consisting of parton distribution functions f P 1 a (x 1 , µ 2 f ) and f P 2 b (x 2 , µ 2 f ) respectively for two incoming protons. The summation over I takes care of the SM and the RS contributions. The hadronic and partonic threshold variables τ and z are defined as They are thus related by τ = x 1 x 2 z. To the all order in strong coupling, the partonic cross-section in the above eq. The pre-factor F (0) takes the following form for neutral vector bosons in SM and spin-2 (RS) boson respectively, where α is the fine structure constant, Q is the invariant mass of the lepton pair, M Z and Γ Z are the mass and the decay width of the Z-boson, c w , s w are sine and cosine of Weinberg angle respectively. The vector and axial-vector part of the weak boson coupling is given as,

8)
Q a being electric charge and T 3 a is the weak isospin of the electron or quarks. Note that the SM contribution consists of contribution from γ and Z as well as their interference. For the invariant mass distribution, however the spin-2 production is decoupled from the SM one [58] and thus there is no interference of them.
The complete SM contribution to DY invariant mass distribution is known up to second order in the strong coupling [10,[92][93][94]. Up to two loops the contribution from RS spin-2 can be written as, n=2 a n s ∆ RS,(n) Notice that computation of the partonic coefficients at the second order requires evaluation of matrix element as well as the proper phase space for the di-lepton pair. Using the method of reverse unitarity [6], where the phase space integrals were converted to loop integrals, the later has been performed very recently in the case of generic spin-2 production [81]. The advantage is that, one then can re-use all the techniques developed for the multi-loop computation. The analytical result obtained in this way is useful for any spin-2 production with universal coupling to the SM. We use these coefficients to predict complete NNLO cross-section for the RS model.

Soft-virtual cross-section
It is important to consider corrections beyond NNLO in order to obtain perturbative stability. The first step as discussed earlier is to calculate the third order soft-virtual crosssection. This can be achieved by calculating the spin-2 form factor at three loops [84] as well as the soft function at the same order. The soft function being maximally non-abelian up to three loops, can be extracted from the known Higgs [25,31] and DY [34,42] results. Using these informations, the third order coefficients for the SV corrections have been obtained [83] for generic spin-2 coupling and have been applied to ADD model to predict the DY distribution to N 3 LO. The analytical coefficients can be also used to predict the SV cross-section in the RS graviton production at the same perturbative order. The SV cross-section at each order consists of plus-distributions and delta function which are most singular terms in the partonic coefficient function. One can write the SV coefficient in terms of perturbative expansion in strong coupling, (2.11) The SV coefficients arise only in flavor diagonal channels i.e. either in qq or gg born process. At any order n of the strong coupling, the SV coefficients have the following structure in terms of plus distributions and delta function,

JHEP07(2020)040
In the SM, only qq contributes wheres for the RS scenario both qq and gg channels contribute to the SV coefficient. Here we present only the leading order term (i.e. n = 0) in this series to follow the overall normalization to the coefficient, Up to two loops, these are already known for quite some time and can be found in [79].
Recently we calculated the three loop pieces [83] for generic spin-2 couplings. In this article we use these three-loop coefficients for the third order phenomenological prediction in the RS model.

Threshold resummation
The NNLO cross-section can be improved with the contribution from threshold logarithms at all orders. In particular when partonic z → 1 the contribution from these singular terms becomes large and unreliable and thus needs to be resummed to all orders. The resummation is usually performed in conjugate space where all the convolutions become simple product and therefore easy to calculate. We follow the standard approach where we evaluate the resummed coefficients in the Mellin-N space. The threshold limit translates there into N → ∞ with N = N exp(γ E ), N being Mellin conjugate to z and γ E is the Euler-Mascheroni constant. Up to the born factor, the partonic resummed cross-section can be organized as [47,52,95], (2.14) The normalization (dσ LO /dQ) is given as, The exponent can be found through the following representation in Mellin space in terms of universal constants the cusp anomalous dimensions A I [96-98, 98-100, 100-108] and constants D I [41,47,83], Performing the integration, one can organize it as a resummed series realizing ln N ∼ 1/a s ,

JHEP07(2020)040
where ω = 2β 0 a s ln N . The expressions for g I i required up to N 3 LL resummation can be found in [46,51,109]. The process dependent coefficient g I 0 can be extracted from the SV results in Mellin-N space and can be written as a perturbative series in strong coupling as, (2.18) We have extracted these coefficients up to the third order from the known SV coefficients at the same order for a generic spin-2 coupling, which can be found in [83]. The first term in eq. (2.17) along with first term in eq. (2.18) define the leading logarithmic (LL) accuracy. Similarly the first two terms in eq. (2.17) along with the terms up to O(a s ) in eq. (2.18) define next-to-leading logarithm (NLL) order and so on. Note that the expansion in eq. (2.17) is different from [41,47], where one organizes the series in terms of N instead of N . Physically both are equivalent in the large-N limit, however numerically it has been seen that the N exponentiation shows better perturbative convergence for DIS [51] as well as in DY [109]. This is because in the second case, the γ E terms are also exponentiated along with the Mellin-N . It is already shown in the case of DIS [51], that the N exponentiation differs as large as 15% at LL compared to N exponentiation, however at the higher logarithmic accuracy these differences get minimized. Recently in the case of DY production [109], it has been observed that the resummed K-factors (at Q = 500 GeV) in the N -space are K LL = 1.133, K NLL = 1.378 and K NNLL = 1.400 while in the N -space they are 1.029, 1.363 and 1.399 respectively. This clearly shows a faster convergence for the N scenario.
The partonic resummed cross-section has to be Mellin-inverted with suitable N -space PDF and finally has to be matched with the known fixed order results. The general expression for the matched cross-section can be written as below: The first term in the above eq. (2.19) represents the fixed order known to N n LO. The last term in the bracket represents the resummed result truncated to the fixed order accuracy i.e. to N n LO to remove all double counting of the singular terms that are already present in fixed order. f i,N (µ 2 f ) are the Mellin space PDFs which can be evolved using publicly available code QCD-Pegasus [110], however we followed the prescription provided in [95] relating N -space PDFs to the derivative of x-space PDF for simplicity and flexibility to using lhapdf [111] libraries. To avoid the Landau pole problem in the Mellin-inversion integration, we have followed the minimal prescription [95] and chosen the contour accordingly. All the necessary analytical ingredients are now available to perform the numerical study which we report in the next section. In our work, all the algebraic computations have been done with the latest version of the symbolic manipulation system Form [112,113].

JHEP07(2020)040 3 Numerical results
In this section we present our numerical results for the di-lepton production cross section in the RS model at LHC. The LO, NLO and NNLO parton level cross sections are convoluted with the respective order by order parton distribution functions (PDF) taken from lhapdf [111]. The corresponding strong coupling constant a s (µ 2 r ) = α s (µ 2 r )/(4π) is also provided by lhapdf. The fine structure constant and the weak mixing angles are chosen to be α em = 1/128 and sin 2 θ w = 0.227 respectively. Here the results are presented for n f = 5 flavors in the massless limit of quarks. The default choice for the centre-of-mass energy of protons is 13 TeV and the choice for the PDF set is MMHT2014 [114]. Except for the scale variations, we have used the factorization (µ f ) and renormalisation (µ r ) scales to be the invariant mass of the di-lepton, i.e. µ f = µ r = Q. We also note that there have been several experimental searches at the LHC for warped extra dimensions in the past, yielding stringent bounds on the RS model parameters, the mass of the first resonance mode (M 1 ) and the coupling strength (c 0 ) [115,116]. Such analyses have already used the K-factors that have been computed in the extra dimension models. In the di-lepton channel using the combined 8 and 13 TeV data at the LHC, the observed 95% CL lower limit on the RS resonance is 1.38 (2.98) TeV forc 0 = 0.01 (0.1) [115]. On the other side, in the di-photon channel using 13 TeV LHC data, the lower limit at 95% CL on the RS first resonance mass is found to be 4.1 TeV forc 0 =0.1 [116].
Here in this work, for our phenomenological study to assess the impact of QCD corrections, we choose M 1 = 1.5 TeV andc 0 = 0.05. The computational details of the QCD corrections presented here are model independent, and a numerical estimate of the theory predictions for any other choice of the model parameters is straight-forward. For completeness, we also study the dependence of the invariant mass distributions on the model parameters considering the recent bounds on M 1 for differentc 0 values.

Fixed order corrections
First we present in figure 1 the contribution from different subprocesses for the pure RS graviton (GR) at NNLO level right at the resonance region by varying the first resonant mass M 1 and keepingc 0 = 0.05. At this order in QCD there are six different subprocesses that contribute for GR case, viz. qq, gg, qg, qq, q 1 q 2 and q 1q2 . Here, the dominant contribution comes from the gg-subprocess and it remains dominant for resonance values as large as 4.5 TeV. The next dominant contribution comes from qg-subprocess but it is negative for this entire mass ranges. This is followed by quark initiated processes with qq being the largest in this category. For a typical choice of first resonance M 1 = 2500, we find that the total cross-section is 0.63 × 10 −5 pb in which the dominant gg subprocess overshoots the value by 151%. The qq, qq, q 1 q 2 and q 1q2 channels contribute in addition 24.7%, 2.7%, 2.2%, 0.7% respectively of the total cross-section. As stated earlier only the qg channel contributes negatively of about −82% of the total cross-section.
Next, we present in figure 2 the di-lepton invariant mass distribution (dσ/dQ) as a function of the invariant mass of the di-lepton Q for GR and for the signal (SM+GR). The width of the resonance depends onc 0 and near the resonance region the signal receives most of the contribution from the pure RS graviton. Far away from this resonance region, the RS contribution is found to be comparable to that of the SM background for Q > 3500 GeV.
In figure 3 we present the K-factor, defined in eq. (3.1), for both GR and the signal cases. In an earlier work we had presented third order SV and resummed results for dilepton production at LHC in the ADD model [83]. We note that it is the same virtual graviton exchange process that contributes both in ADD and RS model.  the difference between these two models arise because of the difference in the summation over the tower of KK gravitons and also in the overall wrapped factor. Consequently the relative weight of the contribution from the gravitons in these two models will be different for different invariant mass region. This results in different mass-dependent K-factors in the ADD and RS model. The NLO corrections for pure RS case at Q = 1000 GeV are found to contribute by about 57% of LO, while NNLO corrections add an additional 18% of LO to the total invariant mass distribution. In table 1 we present the signal K-factors up to NNLO QCD for different Q values. For signal case, the NLO corrections at Q = 1000 GeV contribute by about 34% of LO and NNLO corrections add an additional 6% of LO to the total invariant mass distribution. However, right at the resonance region, these NNLO corrections are found to enhance the production cross section by an additional 20% of LO results. This shows that NNLO corrections are indeed essential for this process in order to make any reliable predictions. In figure 4 we present the di-lepton invariant mass distribution for SM, GR and signal cases. The behavior of the signal K-factor is governed by the respective coupling constants in SM and RS as well as the parton fluxes. As discussed earlier, in the RS case gravity contribution is significant near resonance region and therefore the whole signal K-factor is controlled by RS. In the off resonance region at high Q, both RS and SM contributions are comparable and hence the signal K-factor receives contributions from both RS and SM.   Hence the behavior of the mass dependent K-factor for the signal in the RS model is very distinct from that in the ADD model [83]. We also study the dependence of our results on the RS model parameters, M 1 andc 0 . In figure 5 we present the di-lepton invariant mass distribution by varyingc 0 from 0.03 to 0.1 keeping M 1 fixed at 1.5 TeV in the left panel. We also present the results in the right panel by varying M 1 from 1.5 TeV to 3.0 TeV for a fixedc 0 (0.05). The width of the resonance depends onc 0 , however, right at the resonance this dependence of the production cross section on this couplingc 0 cancels and the height of peak for any given M 1 will be independent ofc 0 . Consequently, the respective NNLO K-factors near the resonance region depend onc 0 but right at the resonance they do not. These signal K-factors are presented in table 2 right at the resonance region for different M 1 values. The NNLO corrections increases the K-factors substantially compared to NLO implying the importance of the higher order correction for this process.  We have considered different sources of theoretical uncertainties in our analysis. First, we considered the uncertainties due to the presence of two unphysical scales µ r and µ f in the theory and then those coming from the non-perturbative parton distribution function in the calculation. For the scale uncertainties we vary µ r and µ f simultaneously from Q/2 to 2Q by putting the constraint that the ratio of unphysical scales is less than 2, as

JHEP07(2020)040
The last condition in eq. 2.91×10 −6 (±6.5%) 3.78×10 −6 (±20.0%) 4.26×10 −6 (±8.0%) 3.20×10 −6 (±11.9%) 3.96×10 −6 (±11.2%) Table 3. Intrinsic PDF uncertainties in the signal at NNLO QCD for different PDF choices are given right at the resonance for different M 1 values. All the results are presented for 13 TeV LHC. The cross sections are given for the central set (n = 0) for each PDF group along with the corresponding intrinsic uncertainties in terms of the percentage. Q = 1500 GeV i.e. right at the first resonance, the scale uncertainties at LO are ±13.3%, at NLO they are ±5.7%, at NNLO it further reduces to ±2.3%. Away from resonance, the uncertainty also decreases order by order. For example at Q = 3000 GeV, the scale uncertainties get reduced from ±14.8% at LO to about ±2.5% at NNLO. In the offresonance region the uncertainty in general increases with increasing Q which can be tamed with inclusion of further higher order terms in the perturbation theory. We also estimate the uncertainties coming from the non-perturbative PDFs. For this we calculate the uncertainty due to the intrinsic errors in the PDFs that result from various experimental errors from the global fits. In this cases we use the PDF sets ABMP16 [117], CT14 [118], MMHT2014 [114], NNPDF31 [119], and PDF4LHC15 [120] provided from the lhapdf. The central predictions for these different PDF groups also differ due to different underlying assumptions in global fits for different groups. We calculate the intrinsic PDF uncertainties using 51 sets for MMHT2014, 57 sets for CT14, 101 sets for NNPDF31, 30 sets for ABMP16 and 31 sets for PDF4LHC15. To this end we use all PDF sets extracted at NNLO level. In table 3 we present these uncertainties for the di-lepton invariant mass distribution to NNLO. We find that around the resonance M 1 = 1500 the PDF uncertainty is well within 5% except for the CT14 which shows relatively increased uncertainty. In the high invariant mass region the uncertainty however increases due to unavailability of sufficient data in those region.
From the above observation we notice that the NNLO corrections for RS production are large enough to truncate the perturbation theory at this order and necessitates the computation of higher order corrections for the convergence of the perturbation series. As a first step beyond NNLO we studied the three-loop SV correction for di-lepton production channel using the universal property of the SV coefficients for generic spin-2 couplings. In figure 7 we present these three loop SV corrections in terms of the corresponding K-factors up to N 3 LO sv as a function of the di-lepton invariant mass for pure RS case (left panel) as well as for the signal (right panel). We use MMHT2014nnlo set for this analysis. These threeloop SV corrections are found to contribute an additional −0.7% of LO to the NNLO result at first resonance M 1 = 1500 GeV for pure RS case, demonstrating a very good convergence of the perturbation theory at this order. In the high invariant mass region away from the RS resonance however we see correction due to third order SV terms is about 1% of LO cross-section. We also note that the three-loop SV corrections are negative in the low Qregion while in the high Q-region they are positive because of threshold enhancement. The 7-point scale uncertainty is seen to increase at the lower invariant mass region whereas as we reach higher invariant mass region this becomes better. To further constraint the scale uncertainty the N 3 LO PDFs are essential at this order. Also the missing sub-leading pieces are important in particular in the low-Q region (see e.g. [40,121]). The µ r uncertainty however is seen to improve in the whole invariant mass region. Keeping µ f = Q = M 1 we observe an uncertainty of ±0.9% around the resonance M 1 = 1500 for the variation in the range (1/2, 2)M 1 .

Resum result
We now move to study the effect of threshold logarithms by resumming them to NNLL accuracy and match to the computed NNLO cross-section in the section 3.1. For this, the same choice of SM and RS model parameters has been used as in the fixed order computation. For the inverse Mellin transformation eq. (2.19), we use c = 1.9. In figure 8 we present the di-lepton invariant mass distribution for GR and for the signal at different logarithmic accuracy.
To quantify these resummation effects, we define the resum K-factors in eq. (3.3) and present the same in table 4 for different Q values. The enhancement due to threshold logarithms for the signal is significant for all Q values, however it is more significant at the resonance region. This is because of the underlying born processes for the graviton production in the RS model. At the born level, the RS graviton can be produced via quarkantiquark annihilation process (DY-like) as well as gluon fusion channel (Higgs-like). It is well known that the QCD corrections, particularly, the threshold enhancement in these two channels are different and are more pronounced for gluon fusion channel. Here, the signal receives contribution from RS (DY-like as well as Higgs-like) and the SM background (DYlike). However, at the resonance region GR dominates over the SM background by several   orders of magnitude and hence the threshold enhancement due to the gluon fusion channel becomes prominent. Far off the resonance region, the signal is essentially dominated by the SM background and assumes DY-like threshold enhancement. For completeness, we present these resummed K-factors for GR case in figure 9.
In order to further study the enhancement due to threshold resummation for the signal, we consider the ratios of the resummed results to the fixed order results defined in eq. (3.4).
We observe that at resonance (Q = M 1 = 1500 GeV), NNLO+NNLL contributes additional 4% enhancement over NNLO. These ratios are presented in figure 9. Moreover, in table 5, we present these resum K-factors right at the resonance for different values of resonance mass M 1 . Next, we estimate the theoretical uncertainties in resummed predictions due to the unphysical scales µ r and µ f as well as due to the non-perturbative PDFs. The conventional 7-point scale uncertainties for the signal are presented in figure 10 for different logarithmic accuracy. At the resonance region, these scale uncertainties are estimated to be about ±17.5%, ±8.1% and ±3.4% at LO+LL, NLO+NLL and NNLO+NNLL respec-JHEP07(2020)040  tively. Moreover, these uncertainties are bit larger than the corresponding ones for the fixed order results presented in figure 6.
To estimate these uncertainties at the NNLO level and beyond, we contrast these scale uncertainties in the resummed results against those in the fixed order results in figure 11 (left panel). For the resummed case, the scale uncertainty at NNLO+NNLL for Q = M 1 is about 3.4% and is larger than the one 2.3% at NNLO level. This increase in scale uncertainty can be understood from the fact that in the resummation formalism only threshold logarithms that are significant in the limit z → 1 have been resummed to all orders in QCD but not the logarithms of unphysical scales. Moreover, it is observed that for Higgs-like processes the resummation does not improve the scale uncertainties over the fixed order ones [43] for any choice of central scales. In the present context, the graviton production at the resonance receives significant contribution from this Higgs-like gluon fusion process and hence the associated large scale uncertainty. The scale uncertainties only due to the renormalization scale µ r are found to get reduced from fixed order NNLO level ±1.2% at Q = 1500 to the resummed NNLO+NNLL level ±0.5% (see right panel of figure 11). However, the factorization scale uncertainty for Q = 1500 GeV is found to increase ±0.88% at NNLO to ±2.58 at NNLO+NNLL level. The overall seven-point scale uncertainties are thus dominated by the factorization scale at this order. Further, we estimate in our predictions the uncertainty due to the non-perturbative PDF inputs. These uncertainties are obtained for each PDF group by systematically calculating the cross section for each of the available sets. These PDF uncertainties are presented for different resonance mass M 1 values in table 6. This uncertainty for the kinematic range considered and the PDF groups studied, is smallest at M 1 = 1500 GeV for NNPDF31 and largest for CT14 at M 1 = 3500 GeV.  Table 6. Intrinsic PDF uncertainties in the signal at NNLO+NNLL QCD for different PDF choices are given right at the resonance for different M 1 values. All the results are presented for 13 TeV LHC. The cross sections are given in terms of pb for the central set (n = 0) for each PDF group along with the corresponding intrinsic uncertainties in terms of the percentage.
Finally, we also note that at this precision level both the electroweak corrections as well as the finite quark mass effects need to be included. They are expected to contribute an additional few percent to the signal cross section.

Conclusions
In the absence of any signature of new physics at the LHC, it is high time to explore possible scenarios where we could make potential discovery of new physics beyond the SM.
In particular the RS model provides to be a very good candidate in the search of massive spin-2 resonances. In the literature, it is found that the NLO QCD corrections for the dilepton production process in this model are quite substantial, implying the need for higher order corrections that augment the search for RS gravitons at collider experiments. In this work, we have studied the NNLO QCD corrections for the di-lepton production process through graviton propagator and have presented the results for the di-lepton invariant mass distribution up to Q values as high as 3.5 TeV. The underlying born contributions for this process receive both DY-like as well as Higgs-like contributions and hence the corresponding QCD corrections for the signal at the resonance region are very significant, while the QCD corrections off the resonance are mostly SM DY-like. This results in K-factors that are strongly dependent on the invariant mass of the di-lepton. We have presented these mass dependent K-factors at NNLO and beyond for 13 TeV LHC. We find that while NLO correction is about 53% of LO, the NNLO correction increases the cross section by additional 21%. The scale uncertainty in the NNLO result at the resonance region also got significantly reduced to as small as 2% for Q = M 1 = 1500 GeV.
Further, we have extended our work to include the important SV corrections at the N 3 LO level. We find that the SV contribution at this order for Q = 1500 GeV is about 0.7% of LO in magnitude but negative in sign, thus demonstrating a very good convergence of the perturbation series. In addition we also studied the threshold resummation by resumming the threshold logarithms to NNLL accuracy and then matching to the fixed order NNLO ones. We find that these resummed results contribute an additional 7% of LO to the NNLO ones. To conclude, we note that our results are most precise theoretical predictions available to date and that these mass dependent K-factors will be useful in the search for RS graviton resonances in the experimental data analysis using di-lepton events at the LHC.