NNLO QCD Corrections to the Drell-Yan Cross Section in Models of TeV-Scale Gravity

The first results on the complete next-to-next-to-leading order (NNLO) Quantum Chromodynamic (QCD) corrections to the production of di-leptons at hadron colliders in large extra dimension models with spin-2 particles are reported in this article. In particular, we have computed these corrections to the invariant mass distribution of the di-leptons taking into account all the partonic subprocesses that contribute at NNLO. In these models, spin-2 particles couple through the energy-momentum tensor of the Standard Model (SM) with the universal coupling strength. The tensorial nature of the interaction and the presence of both quark annihilation and gluon fusion channels at the Born level make it challenging computationally and interesting phenomenologically. We have demonstrated numerically the importance of our results at the Large Hadron Collider (LHC) energies. The two loop corrections contribute an additional 10\% to the total cross section. We find that the QCD corrections are not only large but also important to make the predictions stable under renormalisation and factorisation scale variations providing an opportunity to stringently constrain the parameters of the models with a spin-2 particle.


Contents 1 Introduction
At hadron colliders, the production of a pair of leptons from the decay of electroweak gauge boson is not only a clean process but also it is immensely important for physics studies at the LHC [1,2]. The experimental signature involves two high p T leptons as a result of a neutral gauge boson decay or a single high p T lepton and missing transverse energy in the case of charged counterpart. The parton model ideas intended for the deepinelastic scattering of lepton-proton were formally extended to the proton-proton collisions to produce a pair of leptons (Drell-Yan (DY) process) [3].
The massive electroweak gauge bosons (W ± and Z) were subsequently discovered using this production process. The high production rate and clean experimental final state make the DY process a very important experimental tool and can be used to determine electroweak model parameters. For example, measurements of the W boson production at the Tevatron [4] lead to an accurate determination of the W mass and width. DY processes play an important role in constraining the parton distribution functions (PDF) [5][6][7] of the proton and also serve as luminosity monitor of hadron collider.
While Run-I at the LHC culminated in the discovery of the Higgs boson [8,9], Run-II is currently in operation and the SM is being scrutinised at unprecedented levels of precisions. To fully benefit from the experimental program at the LHC, precise theoretical predictions for both signals of new physics and SM background are very essential. The dominant nextto-leading order (NLO) corrections to the leading order (LO) DY result come from QCD and are large at LHC energies. In addition, an estimate of the theoretical uncertainties due to truncation of the perturbative expansion in the strong coupling constant, a s , reduces on the inclusion of the higher order terms in a s . For the DY process, the NNLO corrections in QCD are available for inclusive cross section [10], rapidity distributions [11,12], fully exclusive distributions including γ-Z interference, the leptonic decay of gauge bosons and finite width effects are also included [13][14][15]. The current accuracy of the DY process is next-to-next-to-next-to-leading order (N 3 LO) corrections to the production cross section near the partonic threshold [16][17][18].
Searches for physics beyond the SM involve looking for deviations from the SM predictions. The excess in the di-photon channel reported by the LHC collaborations [19][20][21][22] triggered enormous interest among theorists to interpret it in terms of a new resonance of mass 750 GeV. While several models with a new particle of mass 750 GeV explaining this excess have already been proposed, the conclusive and the most plausible interpretation is possible only with more data. Though the interpretation with a heavy spin-0 particle could explain the excess in most of the scenarios, the data do not rule out the possibility of a spin-2 particle decaying into a pair of photons. Massive spin-2 particles have been phenomenologically well studied in the context of models with extra spatial dimensions which could be flat as in the large extra dimension model, namely ADD [23][24][25], or warped as in the RS model [26] or any other new physics scenario with spin-2. They couple to all the SM particles universally through energy-momentum tensor of the SM. There are also some studies with non-universal coupling of a spin-2 particle with the particles of the SM [27]. A generic spin-2 particle can also contribute to other production channels, namely di-lepton or di-vector boson productions at the LHC. In this article, we will restrict ourselves to study the invariant mass of di-lepton pair in the ADD model with spin-2 particle. The extension to other production channels is straightforward.
To match the theoretical accuracy of the SM DY process, the di-lepton final states including a spin-2 intermediate state should also be calculated to the same order of accuracy in QCD. Presently for the ADD and RS model, NLO QCD corrections are available for most of the di-final state process with a trivial colour flow viz.: di-lepton [28][29][30], di-photon [31,32], ZZ [33,34] and W + W − [35,36]. In addition, these processes have been extended to NLO+Parton Shower accuracy [37][38][39][40]. In going from LO to NLO the theoretical uncertainties gets reduced, but for most of these processes the renormalisation scale (µ R ) dependence begins at the NLO level and to compensate the µ R -dependence, going beyond NLO is inevitable. Unlike the SM DY the gluon-gluon subprocess starts at LO itself for spin-2 process and so the NLO corrections are large where the spin-2 effects are dominant as compared to the SM.
To go beyond NLO to NNLO, it is prudent to take incremental steps. A general feature of the production of a large invariant mass system in hadronic collisions is that the dominant contributions are often given by the threshold approximation. In [41], the relevant form factors such as the gluon-gluon → spin-2 and quark-antiquark → spin-2 at two-loop level in QCD [42] were computed to obtain threshold corrections at NNLO in QCD to the invariant mass distribution of di-leptons at hadron colliders in ADD model and to a resonant production of a graviton in RS model. In [43], three loop QCD corrections were computed for these form factors in order to study the universal infra-red structure of the QCD amplitudes involving spin-2 particle in the external states. In [44], the two-loop QCD corrections to the amplitudes of massive spin-2 resonance → 3 gluons relevant for the production of a spin-2 particle plus jet were carried out.
Going beyond threshold corrections is inevitable in order to make accurate predictions. The contributions resulting from hard part of the cross section in quark annihilation and gluon fusion channels and those from other partonic channels may contribute significantly. We will demonstrate in this article that this is indeed the case for the invariant mass distribution of di-leptons by explicitly computing the full NNLO QCD corrections. We also find that the contributions from quark-gluon initiated processes both at NLO as well as NNLO levels are not only negative but also large, hence affects the threshold approximation. This is one of the main results of the present paper.
The paper is organised as follows: in Sec. 2, we introduce the effective action that describes the interaction of the spin-2 particles with the SM fields, in particular, the part that is relevant for our computation and then present the theoretical framework to compute the invariant mass of the di-leptons at hadron colliders up to NNLO level in QCD. Sec. 3 is devoted to the methodology employed to compute all the partonic cross sections that contribute. In Sec. 4, we present the numerical impact of our new results. Appendix A and B contain partonic coefficient functions resulting from all the channels up to NNLO level along with some useful identities involving multiple polylogarithms.

The Effective Action
In an effective theory, the spin-2 field, h µν , couples to the SM ones through the conserved SM energy momentum tensor, T SM µν . The effective action [23][24][25][26] describing this interaction reads: where, S SM and S h represent the actions of the SM and spin-2 fields, respectively. κ is a dimensionful coupling constant and T QCD µν is the conserved energy momentum tensor of QCD which is given by g s is the strong coupling constant and ξ is the gauge fixing parameter which is set to 1 for working in Feynman gauge. T a and f abc are the Gell-Mann matrices and structure constants of SU(N) gauge theory, respectively. Presence of the ghost fields, ω a , in the interaction part of the action is a reflection of the fact that spin-2 fields couple to everything democratically. In the action, Eq. (2.1), we have presented only the QCD interaction term as we are interested only in this regime.

Invariant Lepton Pair Mass Distribution dσ/dQ 2
We consider the production of a leptonic pair, l + and l − , through the scattering of two hadrons, represented by H 1 and H 2 : X denotes the final inclusive state. The terms inside the parentheses represent the 4momenta of the corresponding particles. In the QCD improved parton model, the hadronic cross section is related to the partonic one through In the above expression, S is the square of the hadronic center of mass energy which is related to the partonic one,ŝ, throughŝ = x 1 x 2 S. The invariant mass square of the final state leptonic pair, m 2 l + l − is represented through Q 2 . f a and f b are the partonic distribution functions of the initial state partons a and b, respectively. The other parameters are defined as The underlying partonic process corresponding to the hadronic one (2.3) is where, j can be photon (γ * ), Z-boson (Z) or spin-2 particle. X i stands for the real QCD hard radiations from the initial state partons a and b. In perturbative quantum field theory (pQFT), the cross section for the Drell-Yan process can be factored out into partonic (ab → j) and leptonic (j → l + l − ) parts: where, the m + 1 body phase space in n-dimensions is defined as and the quantity L jj ′ →l + l − is given by M ab→jj ′ and M jj ′ →l + l − are the partonic and leptonic part of the matrix elements, respectively. j = j ′ reflects the interference terms between the channel j and j ′ . In the above Eq. (2.6), the sum over Lorentz indices between matrix element squared and the propagators is implicit through a symbol 'dot product'. The propagators are η µν = diag[1, −1, −1, −1, · · · ] and D(Q 2 ), the summation over the virtual Kaluza-Klein (KK) modes in the time like propagators [45] in (4 + d)-dimensions, is The integral I is regulated presumably by a cutoff of the order of M S in ultraviolet (UV) region [45]. This cutoff sets the limit on the applicability of the effective theory. For the DY process, this implies Q < M S . The quantity q 2 = Q 2 = m 2 l + l − . Hence, the computation of the partonic level cross section boils down to the evaluation of partonic and leptonic parts. The leptonic part comes out to be and g µν (q) ≡ η µν − q µ q ν q.q . (2.13) In the above equation, α is the fine structure constant, c W ≡ cos θ W , s W ≡ sin θ W and θ w is the weak mixing angle. g V f and g A f can be expressed in terms of charge Q f of the fermions (f ) i.e. quarks, leptons and weak isospin T 3 f : Hence, the hadronic cross section (2.4) can be rewritten as where, the hadronic structure function W is As a consequence, the computation of the Q 2 distribution of the di-lepton pairs requires the evaluation of the integrals in a suitable frame over dP S m+1 and dz after substituting the matrix element squared |M ab→jj ′ | 2 T jj ′ (q) in Eq. (2.16). We define the bare partonic coefficient function∆ jj ′ ab z, Q 2 , 1/ǫ througĥ The partonic cross section receives contributions from two different classes of processes: first one happens through a virtual photon or a Z-boson whereas the second one contains a spin-2 particle in the intermediate state. Interestingly, on performing the phase space integration, the interference term between the two classes of diagrams up to NNLO identically vanishes, this was earlier noted to NLO [28]. Hence, our result does not receive any contribution from the interference terms.
In case of spin-2 appearing as an intermediate state, at leading order (LO) we can have gluon initiated process as well, in addition to the quark initiated one. Upon inclusion of the spin-2 contribution, at the LO we have (See Fig. 1) (2.20) Figure 1: Leading order processes for the DY Beyond LO, the contributions arise from virtual as well as real emission diagrams. At next-to-leading order (NLO) in QCD, we have Similarly at the next-to-next-to-leading order (NNLO) in QCD, the contributions can arise from three different categories: double-real emissions, real-virtual and double virtual diagrams. The processes which belong to the double-real emissions are The processes which contribute in real-virtual are g +q → γ * /Z/h +q + one loop (2.23) and the pure double virtual diagrams are The evaluation of the virtual as well as real emission diagrams exhibits divergences of three kinds: ultraviolet (UV), soft and collinear (IR). In particular, diagrams involving virtual particles give rise to all the above-mentioned divergences whereas the real emissions cause only soft and collinear ones. We regulate the UV as well as IR divergences using dimensional regularisation where the space-time dimensions n is chosen to be equal to 4 + ǫ. All the divergences manifest themselves as the poles in dimensional regularisation parameter ǫ: 1/ǫ α with α ∈ [1,4]. In MS, the UV poles are removed through strong coupling constant renormalisation usinĝ where, and µ is the scale introduced to keep the unrenormalised strong coupling constantâ s dimensionless in n-dimensions. The corresponding renormalisation scale is denoted by µ R . β i 's are the coefficients of QCD β-function [46][47][48][49][50]. Since, the spin-2 particles couple to the SM ones through the conserved energy momentum tensor, the universal gravitational coupling constant κ is protected from any ultraviolet (UV) renormalisation. Hence, there is no additional UV renormalisation required other than the strong coupling constant renormalisation. The soft divergences arising from virtual diagrams cancel exactly against the same coming from real emission ones, thanks to the Kinoshita-Lee-Nauenberg (KLN) theorem [51,52]. The collinear divergences are removed through mass factorisation, performed at the factorisation scale µ F : In the above expression,∆ ≡σ/z is the bare partonic coefficient function and the corresponding one after performing the mass factorisation is denoted by ∆. Further we have dropped the double index jj ′ from the partonic coefficient function (see Eq. (2.18)) because of the vanishing interference terms between the two classes of diagrams and replace it by the single index i instead. The mass factorisation kernel in the modified minimal subtraction (MS) scheme is given by ab are the Altarelli-Parisi splitting functions [53][54][55][56][57]. The symbol ⊗ stands for the convolution: (2.29) Expanding the unrenormalised coefficient function in Eq. (2.18) and the mass factorised one in Eq. (2.27) in powers of strong coupling constant aŝ and using the Eq. (2.28), we can get all the contributions to NNLO arising from all the subprocesses: To arrive at the above set of results (2.31), we have made use of and the presence of the n f in RHS comes from the sum over flavors. The superscript S in the last two equations of (2.31) denotes the flavour singlet part of the AP kernel. From the results of the bare coefficient functions and the known splitting functions, we can obtain the finite ∆ i ab using the above Eq. (2.31). This in turn gives us the Q 2 distribution of the leptonic pair in the DY process: In the above expression and the renormalised partonic distributions are In this article, we extend this distribution of the DY pair to NNLO QCD from the existing NLO result [28] in models of TeV scale gravity. The contributions arising from solely SM are already available in the literature [10,[58][59][60]. The missing parts, namely, the contributions coming from the presence of the spin-2 particles, ∆ h,(2) ab are computed in this article. In the next Sec. 3, we discuss the methodology of this computation in great details.

Methodology
The computation of the partonic cross section beyond leading order consists of the evaluation of the loop integrals arising from the virtual diagrams and the phase space integrals.
The developments of the techniques to evaluate the former one takes place quite rapidly compared to the latter one. In the very first computation of the NNLO QCD correction to the DY pair production in [10], the phase space integrals were performed through evaluation of the two parametric and two angular integrations in three different frames. Later, to calculate the inclusive production cross section of the Higgs boson three different techniques were employed. In [61], the partonic cross section was obtained by performing an expansion around the soft limit. In the meantime a completely new and elegant formalism was developed in [62] by Anastasiou and Melnikov to get the same result. The phase space integrals were converted to loop integrals by using the idea of reverse unitarity. So, the evaluation of the phase space integrals boils down to the evaluation of the loop integrals. Hundreds of different loop integrals were reduced to only a few number of master integrals (MIs) by making use of the integration-by-parts (IBP) [63,64] and Lorentz invariance (LI) [65] identities. The resultant MIs were computed using the techniques of differential equations to arrive at the final result. The same result was again reproduced in [66] using the conventional method of evaluating loop and phase space integrals. The method of reverse unitarity was latter employed to obtain the state-of-the-art result, namely N 3 LO QCD corrections to the inclusive Higgs boson production [67][68][69]. In this article, we use the formalism developed in [62] to calculate the partonic cross section of the DY pair production through intermediate spin-2 particle at NNLO QCD. In this section, we demonstrate this methodology in great details.
At this order we need to calculate three different contributions which are mentioned in the last section: • double-real: the self-interference of the tree level amplitudes for the process listed in Eq. (2.22). For example, for the process q +q → h + q +q, presented through the    All the required Feynman diagrams are generated symbolically using computer package QGRAF [70]. The raw output is converted to a suitable format using in-house code written in FORM [71,72] for our further computation. In subsequent subsections, we describe the techniques to evaluate the above three categories.

Double-Real
We take a sample double real emission diagram for illustrating the methodology [62,73] to handle phase space integrals: q +q → h + q +q where, δ + (q 2 − m 2 ) ≡ δ(q 2 − m 2 )θ(q 0 ). According to Cutkosky rules [74], the δ + functions can be replaced by the difference between two propagators with opposite prescriptions for their imaginary parts: with ε → 0. Upon this substitution, the square of the diagram, depicted through Fig. 5 becomes equivalent to the forward scattering amplitude, presented in Fig. 6, where, the blue dotted line denotes the cut propagators which should be replaced by the RHS of Eq. (3.1). We begin our computation by evaluating the normal Born square of the above diagram (5-external onshell legs) where the sum over colors and spins are performed. With the final answer, we multiply the phase space factor which contains the three δ + functions corresponding to the final state particles. Moreover, to convert it into a cut two-loop Feynman diagram through the application of reverse unitarity, we replace the δ + functions by the difference of the corresponding propagators using Eq. (3.1). As a consequence, We make use of the IBP and LI identities to reduce this two loop diagram into a set of MIs. Since, the sign of the imaginary parts of the cut propagators are irrelevant for the above identities, the two terms of those propagators which are differed by the different prescriptions of the imaginary parts give rise to same IBP relations. Each of these two terms have the same form of the IBP relations as the original two-loop integral without the cut. Hence, instead of considering the two terms, we can take only one term. This is equivalent to substituting the δ + functions by its first propagator from the RHS of Eq. (3.1). Once the reduction is done, we must put those MIs to zero which do not contain any of the three cut propagators. In other words, the MIs which contain all the three cut propagators are the only ones to contribute to the original phase space integrals owing to the Eq. (3.1). While performing the reduction using the Mathematica based package LiteRed [75,76], we make sure not to apply any transformation on the momenta of the cut propagators which essentially helps to keep intact the cut propagators in its original form even in the MIs. At the end, the δ + functions need to be reinstated in place of all the cut propagators which leads us to the final set of phase space MIs. These integrals are identified with the ones appearing as phase space MIs for the evaluation of the NNLO QCD correction to the inclusive production cross section of the Higgs boson which are obtained in the article [77]. Same set of MIs were also evaluated in [78].

Real-Virtual
The real-virtual diagrams can also be treated in similar way as double real ones. We illustrate it through the example q +q → h + g + 1-loop. Among all the diagrams let us consider the one depicted through the Fig. 7. In the diagram, k 1 is the loop momentum for the virtual loop.
To compute the interference term all the external legs are treated as onshell particles and the polarisation sum of the external gluons is carried out in axial gauge to ensure the exclusion of the unphysical degrees of freedom. We include the ghost loops to cancel the unphysical degrees of freedom of the internal gluons present in the virtual loops. Similar to the double real emission case, we replace the δ + functions arising from the phase space factor by its corresponding propagators using Eq. (3.1). Consequently, the integral manifests itself as a two loop integral with two cut propagators. This is represented through Fig. 8 where the blue vertical dotted line denotes the presence of the cut propagators. Using IBP and LI identities, we reduce the two loop integral to a set of MIs. The MIs deprived of any of the two cut propagators do not contribute to the phase space integral. Hence, we keep only those MIs which contain both the cut propagators. Throughout the calculation including the process of reduction, we preserve the cut propagators to its original form. After the reduction, both the cut propagators are replaced by its corresponding δ + functions to convert it into phase space MIs which are identified against the ones presented in [77].

Double Virtual
The interference between the double virtual and Born level diagrams can be thought of as integrals with a single cut through the propagator of the spin-2 particle. For the process q +q → h + 2-loops, we consider one of the diagrams depicted through Fig. 9. This can Feynman diagrams. In this present article, the computations of the double real and realvirtual contributions are performed mostly using our in-house codes written in FORM [71,72] and Mathematica. The color simplification is done in general SU(N) gauge theory. The Dirac and Lorentz algebra are carried out in n-dimensions (n = 4+ǫ). After performing the IBP reduction of the phase space integrals to reduce these to a smaller set of MIs following the techniques described above, we borrow the analytical results of these MIs from [77] to get the final answer in powers of ǫ. However, instead if directly using the results of the MIs presented in [77], we make use of some identities to convert the expressions into a form which is manifestly real. The results of the two loop virtual diagrams are available from [42] which were computed by some of us. Using the results of all the subprocesses belonging to the above discussed three categories and performing the appropriate mass factorisation using Eq. (2.27), we get the completely finite partonic cross sections or partonic coefficient functions at NNLO QCD. All the final results of the partonic coefficient functions involving spin-2 particle, Eq. (2.31), are presented in the Appendix A.

Numerical Implications
In this section, we present the numerical impact of two-loop QCD corrections on the dilepton production in ADD model at the LHC. The LO, NLO and NNLO corrected hadronic cross sections are obtained by convoluting the partonic coefficient functions order-by-order with the corresponding parton distribution functions (PDFs) taken from lhapdf [79]. We have used the strong coupling constant a s supplied by the corresponding PDF set. The fine structure constant α em = 1/128 and the weak mixing angle sin 2 θ W = 0.227. The results are presented for n f = 5 flavours and in the massless limit of quarks. Unless mentioned otherwise, our default choice of the PDF set is MSTW2008lo/nlo/nnlo. Except for studying the scale variations, the factorisation and the renormalisation scales are set equal to the invariant mass of the di-lepton, i.e., µ F = µ R = Q. Before proceeding further, we note that there have been a series of experimental searches for large extra dimensions in the past at both the Tevatron and LHC. Consequently, stringent bounds have been obtained on the scale M s of the ADD model as a function of the number of extra dimensions d [80,81]. For the illustration of the impact of QCD corrections, we choose the model parameters to be M s = 4 TeV and d = 3.
Let us begin by discussing the relative contributions of various partonic channels that contribute to the hadronic cross section at NNLO level. The contributions from individual channels are not physical while their sum is. The bare partonic cross sections are ill defined due to the presence of infra-red divergences and are removed by mass factorisation in a scheme dependent way. Hence, the resulting channel-wise contributions depend on the scheme, which in our case is MS. In the Fig. 10, we present the Q distributions for various subprocesses at NNLO in the ADD model along with the contribution from SM at NNLO [10,59,60].
At LO, the quark anti-quark initiated sub-process (qq) contributes both in the SM and in the ADD model. However, the gluon fusion sub-process (gg) starts contributing at the LO in the ADD model unlike in the SM where its contribution begins at NNLO. We note that the contributions arising from the gg sub-process in the ADD model dominates over the rest, because of the large gluon flux at the LHC. Recall that the production cross section for the Higgs boson at the LHC is also dominated by gluon fusion sub-process. The crucial difference between these two production channels is the presence of strong coupling constant a s (µ R ) at the leading order for the Higgs boson production cross section. The other interesting aspect that one can not ignore is the numerical impact of quark-gluon (qg) sub-process beyond LO. The major difference between them is that in the ADD model at NLO, through mass factorisation, it receives collinear subtraction terms due to the presence of qq and gg born sub-processes, whereas in the SM it is due to only qq Born sub-process. Irrespective of this difference, the qg sub-process contribution both in the SM and in the ADD model is found to be negative but significantly large in magnitude. The same trend continues even at NNLO. Particularly, we notice that the NNLO QCD corrections from qg sub-process are considerably larger in magnitude than the sum of all the quark initiated sub-processes (qq, qq, q 1 q 2 , q 1q2 ). The other channels, as can be seen from the Fig. 10, contribute very little to the total inclusive cross section but they are important to stabilise the cross section under renormalisation and factorisation scale variations through renormalisation group equations. A generic pattern in all of these sub-processes is that their contributions increase with Q, simply because of the increase in the number of accessible KK-modes with Q.
We next move on to the Fig. 11 where in the left panel we present dσ/dQ as a function of invariant mass Q at LO, NLO and NNLO for ADD model (i.e. setting the SM contributions to zero). We find that the contribution from the interference terms between the SM and spin-2 is zero. It is also observed that the contributions arising from the O(a 2 s ) increase the NLO cross section moderately. In the right panel, we have plotted the K-factors that are defined as The NLO QCD corrections here increase the LO cross sections by about 68% for Q = 1.5 TeV, while the NNLO corrections that are still reasonably large contribute an additional 12% (K 1 = 1.68 and K 2 = 1.80). With this considerably large contributions, the reliability of perturbative QCD calls for the computations beyond NNLO. The K-factors depend on the invariant mass through the logarithm corrections both in partonic cross sections as well as in the evolution of PDFs. Hence one is discouraged to use the constant Kfactor for constraining the model parameters. Finally, we would like to make a remark that the conservative estimate of the K-factor for the Drell-Yan production in ADD model resembles closely to that of the Higgs boson production. However, because of the large negative contribution from the qg sub-process, the exact values of the K-factors differ in these two cases. In any case, we note that K 2 in ADD model alone is bigger than the corresponding one for the SM simply because of the dominance of gg sub-process over others.
In the left panel of Fig. 12, we present the NNLO cross sections for the SM, spin-2 (GR) and the signal (SM+GR) together with the corresponding NNLO K-factor i.e. K 2 in the right panel. Because of the increasing cross sections with Q in the ADD model, they will dominate the decreasing SM ones at some invariant mass Q 0 for any given choice of the model parameters. For our default choice of model parameters, Q 0 is about 1.4 TeV. This implies that the signal is dominated by SM contributions well below Q 0 and by ADD model contributions well above this Q 0 . In the region closer to this Q 0 , which itself depends on the choice of the model parameters, the signal K-factor receives contributions both from SM and ADD model. From now on, we will focus on the signal contributions together with the corresponding SM background owing to their importance in the experimental searches for extra dimensions. In Fig. 13, we present the results for the Q-distributions in the SM and ADD model order by order in the perturbation theory in the left panel and the corresponding K-factors in the right panel. seen in the left panel of Fig. 14. However, for far beyond this Q 0 , the SM contribution can be neglected altogether and hence the SM+ADD K-factor assumes just the pure ADD K-factor that is insensitive to the choice of the the model parameters. Hence far beyond Q 0 , the SM+ADD K-factors tend to converge to each other as can be seen in the right panel.
We also study the dependence on the number of extra dimensions, see Fig. 15. A similar explanation can be given as for the M s variation except noting that the cross sections in SM+ADD decrease with increase in the number of extra dimensions d. The leading order prediction is only a crude estimate of the true cross section. In our case, the LO prediction depends strongly on the factorisation scale µ F through the parton distribution functions.
It is often mild for the quark initiated processes while it is strong for the gluon initiated process. The dependence on the µ F scale starts getting reduced at higher orders leaving a residual scale dependence that is proportional to a n s , n > 1. At NLO level, for the first time the strong coupling constant a s (µ R ) enters our calculation. Since it depends on the renormalisation scale µ R , the result up to NLO level will now become sensitive to choice of µ R . Hence, at NLO, while the factorisation scale dependence gets reduced, the renormalisation scale dependence crops up. Renormalisation group equation ensures that the inclusion of more and more higher order terms in the perturbation theory will reduce its dependence and it will eventually go away if we know the result to all orders in perturbation theory. A similar statement can be made for the factorisation scale dependence as well thanks to the fact that the factorised hadronic cross section is independent of µ F . In order to demonstrate the reduction in the scale dependence, we have plotted the dσ/dQ in the Fig. 16 at a fixed value of Q = 1.5 TeV, the choice where the new physics dominates, as function of µ R (left panel), µ F (right panel) and then µ = µ F = µ R (see Fig. 17. in the range between Q/10 to 10Q, for wider scale variations. As expected, we find that inclusion of higher terms in the perturbation theory indeed reduce the dependence on these unphysical scales. In Fig. 18, we have presented the predictions for the invariant mass distribution for various center of mass energies, namely 7, 8, 13 and 14 TeV at the LHC. As the energy increases, the parton fluxes particularly the gluon flux will increase and hence the sensitivity to the ADD model will also go up. Consequently, both the NNLO SM+ADD cross sections (left panel) and the corresponding signal K-factors (right panel) will increase with the center of mass energy.
In addition to the choice of scale, the choice of PDFs do affect the predictions significantly. The precise value of the strong coupling constant consistent with a given PDF set also influences the prediction. In order to study these effects, we have plotted the cross sections, in the left panel of Fig. 19, using various PDF sets such as ABM12, CT10, NNPDF3.0 against the one which uses MSTW (this choice is arbitrary). This approximately gives the sensitivity to the choice of PDF sets and a s (µ R ), as well as the estimate of the error on the predictions. It is also worth noting here that although the difference in the cross sections is directly related to the difference in the parton fluxes from different PDF sets, the K-factors may not show a similar pattern as that of the cross sections simply because PDFs of different orders enter in the ration of K-factors, as can be seen in the right panel of Fig. 19.
Finally, we address the impact of soft-plus corrections on our fixed order predictions. Note that for ADD, the numerical impact of soft-plus-virtual (SV) were already reported in [41]. Now that we have a complete result at NNLO level, it is important to study the validity of SV approximation. As mentioned before that the gg initiated sub-process in the pure spin-2 case is similar to the SM Higgs production in gluon fusion channel. For the latter case, the SV corrections (or rather with the modified parton fluxes) are found to be a very good approximate for the fixed order results. This indeed is the case even for our ADD model predictions provided we just take only the gg initiated subprocesses. In addition, if we use the modified SV approximation as described in [82], we find that it is closer to the exact result, resulting from gg subprocesses alone (see Fig. 20). Inclusion of qg initiated sub processes spoil this approximation as their contribution is negative and significantly large. Hence, the SV approximation at a 2 s does not seem to be working very well unlike in the Higgs production in gluon fusion.

Conclusions
In summary, we have performed the very first calculation involving a massive spin-2 particle at NNLO level in QCD for the production of a pair of leptons at hadron colliders. We have included all the relevant sub-processes that can contribute to the invariant mass distribution of the di-leptons. The methodology of reverse unitarity and IBP identities are systematically employed to achieve it. Unlike the DY process within the SM, the spin-2 mediated processes are dominated by the gluon initiated ones due to the large gluon flux at the LHC. In addition, the quark-gluon initiated sub-processes have negative but significantly large contributions at NNLO. The corrections at various orders are quantified through their respective K-factors (1.54 at NLO and 1.62 at NNLO). We find, that the corrections are not only large but also important to stabilise the predictions with respect to the unphysical renormalisation and factorisation scales. The scale uncertainties get reduced to 29% at NLO from 71% at LO, which further gets reduced to about 8% at NNLO. The extensions to the scenarios where a spin-2 particle couples differently to various SM particles are straightforward.

A Results of the Partonic Cross Sections
In this appendix, we present the renormalised and finite partonic coefficient functions involving spin-2 particles, ∆ h,(k) ab z, Q 2 , µ 2 F in Eq. (2.30), up to NNLO QCD (k = 0, 1, 2). The results at NLO are in agreement with the existing ones [28]. The soft-virtual corrections i.e. the contributions arising from the soft gluon emissions at NNLO were computed in [41]. Our results are also consistent with these ones. Below we present all of our findings after normalising the components of the coefficient functions by the leading order results: and all the log i (1−z) 1−z terms should be understood as distributions, D i with The results are obtained as