Precision QCD phenomenology of exotic spin-2 search at the LHC

The complete next-to-next-to leading order (NNLO) QCD correction matched with next-to-next-to leading logarithm (NNLL) has been studied for Drell-Yan production through spin-2 particle at the Large hadron collider (LHC). We consider generic spin-2 particle which couples differently to the quarks and the gluons (non-universal scenario). The threshold enhanced analytical coefficient has been obtained up to third order exploiting the universality of the soft function as well as the process dependent form factors at the same order. We performed a detailed phenomenological analysis and give prediction for the 13 TeV LHC for the search of such BSM signature. We found that the matched correction at the second order gives sizeable corrections over wide range of invariant mass of the lepton pair. The scale variation also stabilizes at this order and reduces to 4%. As a by-product we also provide ingredients for third order soft-virtual (SV) prediction as well as resummation and study the impact on LHC searches.


Contents 1 Introduction
Following the discovery of the Higgs boson at the Large Hadron Collider (LHC) [1,2], much of the focus is given to the search for signals of possible beyond Standard Model (BSM) Physics scenarios. There have been many such models proposed to address various problems in the Standard Model (SM) viz. supersymmetry, extra dimensions, little Higgs, technicolor etc. [3]. These signals can be seen either in the form of smooth deviations from the SM predictions due to contact interactions, or in the form of resonances due to the existence of new heavy particles that these models predict.
In the context of LHC, there have been various observable like invariant mass and transverse momentum of the final states involving either leptons or photons that have a very large detection efficiency (≥ 90%) at the ATLAS and CMS detectors. Among such observable, the dilepton invariant mass is of particular interest as the lepton (electron and muon) signals are very clean and the invariant mass (Q) of the lepton pair has been measured up to 2 TeV to a very good accuracy [4][5][6]. It is to be noted that with the availability of such increasing precision in the experimental data, an adequate theoretical prediction to that accuracy is necessary. After the computation of the Higgs production to N 3 LO accuracy [7], such a high precision has been achieved for the dilepton production at hadron colliders, thanks to the recent computation in the photon mediated channel [8] and the massive gauge boson W ± production [9]. Such a precision both in the theoretical as well as in the experimental frontiers facilitates a robust search for BSM signals.
In this work, we focus on one particular BSM scenario where generic massive spin-2 fields interact with the SM ones. Such a model is well motivated in the context of search for spin-2 particles (graviton) that couple to SM bosons and fermions differently. The phenomenology of such a graviton is similar to some extent to that of the RS model [10] where the massive graviton couples with equal strength to fermions and bosons in the SM. However unlike the RS model, the parameter space of the non-universal case is much flexible and less constrained at the LHC searches so far. In the context of Higgs characterisation, this model has been studied extensively in the di-boson channels [11]. Later, a complete automation has been done at NLO [12] in the FeynRules [13]-MadGraph5 aMC@NLO [14] framework. A detailed phenomenology has been performed there for arbitrary values of mass of graviton and its couplings to the SM including parton shower effects as well. It has been observed that the K-factors for different channels give sizeable contribution depending on model parameters as well as phase space region. This necessitates further higher order QCD corrections for this model. The first complete next-to-next-to leading order (NNLO) computation has been performed [15] in the dilepton channel. A detailed analysis has been presented there along with the comparison with universal scenarios. It has been observed that the K-factors at NLO and NNLO are significantly different from those of the ADD [16] or RS [10] cases. Note that both in ADD and RS models, the Kfactors for the dilepton production are found to be much larger than those in the SM. This behaviour is well understood because of the presence of the additional gluon fusion process (Higgs-like) at the born level in addition to the quark annihilation process (Drell-Yan-like) for graviton production, which is absent in the SM where the vector boson production takes place via only quark annihilation process at the LO. In the present model, the graviton production takes place at the LO via Higgs-like and DY-like processes similar to the ADD or RS model, but with different couplings to quarks (k q ) and gluons (k g ). This particular feature of this model controls the size of the contributions from different channels at higher orders in different parameter space.
The observed discrepancy between the K-factors for the graviton production in the universal and the non-universal cases can be explained to some extent with the help of additional parton radiations at higher orders. For example, a DY-like process with small k q coupling can receive large corrections when the graviton couples with large coupling k g to gluon emitted from the quark lines at higher orders. Similar is the case for Higgslike graviton process. The contribution of other parton level subprocesses initiated by qg starting from NLO and qq from NNLO onwards will have further noticeable contributions. For example at NLO, the qq subprocess gets corrections a few times larger than the LO subprocess for the choice of parameters where k q k g simply due to the fact that high k g contribution appear at NLO from qq subprocess itself. The contribution from the qg channel could be as large as qq for different parameter choices. The presence of such model dependent large contributions at higher orders for an inclusive process can not be estimated with a simple scaling of the lower order cross sections with the conventional (NLO/NNLO) K-factors computed either in the SM or in the RS model. This motivates not only a detailed phenomenological study at higher orders in QCD but also questions the convergence of the perturbation theory. This is particularly important in the higher mass region where one would normally expect the BSM effect. This region suffers from large threshold corrections and needs a procedure to resum them to all orders. The resummation procedure has been pioneered [17][18][19][20][21][22][23][24] over the last few decades to systematically include those large contributions in the threshold region from all orders providing better predictability of the perturbative series.
At the LHC energies where the parton fluxes are large, the threshold logarithms can give sizeable contribution and through resummation these lead to better perturbative convergence as well as theoretical uncertainties for various processes [19,21,23,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. Recently the effect of these large threshold logarithms has been studied in ADD and RS models [43,44] and better perturbative results have been predicted at N 3 LL(NNLL) accuracies. In view of this, we systematically include the threshold resummation effects in DY production for non-universal spin-2 model and present a detailed phenomenological results up to NNLL accuracy. Moreover we have also computed the necessary ingredients to perform resummation up to N 3 LL.
The paper is organized as follows. In the sec.
(2), we develop the theoretical formalism and collect important formulas needed for DY production at the LHC. In sec. (3), we present detailed numerical results for the resummed cross-section matched to the fixed order at the NNLO level. We also present a quantitative estimate of the threshold corrections at the fixed order N 3 LO level and finally give a summary of findings in sec. (4).

Theoretical Framework
The interaction of spin-2 field (h µν ) with the SM ones is described through the following effective action, Sum over repeated indices is implied in all equations. The interaction action is splitted in a way so that I = G contains purely gauge terms whereas I = Q contains the fermionic sector and its gauge interactions terms. This decomposition however is not unique and one can shuffle gauge invariant terms between these two.k I are the unrenormalized coupling constants with which the spin-2 field couples to the operatorsÔ I,µν . These gauge invariant operators have the following expressions in terms of unrenormalized quarks and gluon fields 1 , .ω andξ are the ghost field and gauge fixing parameter respectively. Notice that the sum of these operators are protected against radiative corrections to all orders due to the fact that it corresponds to the conserved energy momentum tensor of QCD. However individually they are not conserved and requires additional UV renormalization. These operators are closed under UV renormalization and hence renormalization can be performed [45,46] through a mixing matrix Z.
The unrenormalized operatorsÔ µν I (x) can be renormalized in a closed form [45,46] using renormalization matrix Z IJ as We use dimensional regularization in d = 4 − 2 dimensions to regulate both UV and IR divergences. In d-dimensions, the bare strong coupling can be related to the renormalized one asâ where S d is the spherical factor. Z(a s (µ r )) is the strong coupling renormalization constant which takes the following form in d = 4 − 2 dimensions.
The UV renormalization constants Z IJ satisfy the following RGE where γ IJ are the UV anomalous dimensions. The interation action can be expressed in terms of renormalized quatities as well i.e. in terms of renormalized couplings (κ I ) and renormalized operators (O I ) 2 as The unrenormalized coupling constants are then related to the renormalized ones with the transpose of the same renormalization constants 3 Z IJ aŝ The solution of UV renormalization constants RGE (eq. (2.7)) matrix leads to the following expansion in terms of renormalized strong coupling up to the third order, Note that in a pure gauge theory (n f = 0), the operator O µν G is conserved. 3 We confirmed a minor typo in eq. (2.11) of [47] which however finally does not affect the results there.
All the relevant anomalous dimensions (γ IJ ) are extracted to three loops [15] from the bare quark and gluon form factors and by simply claiming the universality of the infrared divergences at the same order. For the sake of completeness we have collected those in the appendix appendix A. Note that one can equivalently perform coupling constant renormalization using eq. (2.9) instead of operator renormalization in eq. (2.4) to remove the UV divergences from the form factor. We now turn our discussion into the DY production cross-section at the LHC which takes the following form in terms of partonic coefficient function ∆ I and luminosity L, The luminosity function consists of non-perturbative parton distribution functions. The prefactor F (0) I is given for SM and for spin-2 (denoted as GR) channels as follows Here α represents the fine structure constant, c w , s w are the sine and cosine of the Weinberg angle and M Z and Γ Z are the mass and width of the Z boson respectively. The vector and axial coupling of the weak boson is given as Λ GR in eq. (2.12) is the cut-off scale of the spin-2 theory and k I are introduced for convenience and are defined as k I = √ 2κ I /Λ. The propagator for GR theory with mass M G is given as The analytic form of the partonic coefficient function is known for some time to the second order for all subprocesses. The partonic coefficients can be decomposed into the following form for the gluon and quark-antiquark initiated processes: ∆ I ab = ∆ I ab,sv + ∆ I ab,reg , ab ∈ {gg, qq} . (2.16) The first term in the above equation is termed as the soft-virtual term which consists of form factor contribution as well as soft gluon radiation. The perturbative expansion of the SV coefficients are given as, N I are some overall prefactors taken out from the SV coefficient. The N I prefactors are given as .

(2.18)
Note that with this normalisation in eq. (2.18), the SV coefficient at LO simply becomes δ(1 − z). The main ingredients to compute the SV coefficients (after the PDF renormalization) are the form factor and universal soft distribution functions. As mentioned earlier, the form factors are available in the literature up to three loops [47]. We have used the eq. (2.9) to remove the UV divergences from the bare form factors. The UV finite form factor then contains the IR divergences originating from the soft-collinear region. After performing collinear renormalization of parton distribution functions, the remaining IR divergences are cancelled upon inclusion of universal soft distributions. The universal soft distributions are already known in the literature up to three loops [23,[48][49][50][51]. This finally leads to finite SV coefficients up to third order which contains δ(1 − z) and plus-distributions in partonic threshold variable z = Q 2 /ŝ. The new third order SV coefficients are collected in the appendix C. In order to better describe the fixed order cross-section particularly in the threshold region, one needs to resum threshold enhanced logarithms to all orders. In the threshold region the partonic z → 1 induces large singular terms from delta function and plus distributions. The resummation is performed in Mellin space where the threshold limit z → 1 translates into N → ∞ and the large logarithms appear in the form of ln i N (with ). The partonic coefficient takes the following exponential form (up to normalisation by born factor.) The born normalisation factor is given as The prefactor g I 0,p can be expanded in terms of strong coupling as The universal resummed exponent in eq. (2.19) can be found from the Mellin-transform of well known universal cusp anomalous dimensions A I p [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66] and constant D I p [19,21,43], This integration can be performed after expanding the anomalous dimensions in strong coupling. The resulting function can be organised in a perturbative series as where ω = 2β 0 a s ln N . Note that in the context of resummation a s ln N ∼ O(1). The successive terms in the above series along with the same for g I 0,p in eq. (2.21) define the resummation accuracies LL, NLL etc. The required exponents can be found in [19,21,40,41]. The coefficients g I 0,p which contains the process dependent information can be extracted from the newly computed SV coefficients by just transforming all the δ(1 − z) and plus distributions into the Mellin space and dropping any N −dependent terms. Note that these coefficients differ significantly compared to the universal scenario due to the fact that they now explicitly depend on the couplings. However in the limit when these couplings (k q , k g ) are same, they should reproduce the exact coefficients for the universal case. This serves a crucial check for our computation. We collect these newly computed coefficients up to the third order in the appendix B.
The resummed expression in eq. (2.19) obviously does not contain any contributions which are not enhanced in the threshold limit particularly the regular contributions. These regular pieces from the quark or gluon initiated process or from other subprocesses are also important to describe the full phase space. A consistent way to include those contribution is generally performed through matching procedure as The first term on the right hand side is the exact fixed order cross-section at N n LO. The terms in the parentheses resums all the large logarithms avoiding any double counting from the fixed order counter part. Note that resummation is performed in Mellin space and finally one has to perform a Mellin-inversion. In particular one has to avoid the Landau pole while performing the inversion. This can be done using Minimal prescription [67] or Borel prescription [68,69]. We followed the former 4 approach and fixed the contour accordingly.
Finally we conclude this section with the note that although generally the large Mellin-N logarithms are exponentiated, it is also possible to exponentiate some part or all g I 0,p [31,40,41,[71][72][73]. However in this article we sticked to the standard approach and studied its effect at the LHC.

Numerical Results
We now turn to the discussion on resummed numerical result up to NNLO+NNLL accuracy in QCD for dilepton production via spin-2 particle at the LHC with non-universal coupling of spin-2 particle with gauge bosons and fermions. We remind the coupling strength of spin-2 particles to bosons through κ G = √ 2k g /Λ and to fermions through κ Q = √ 2k q /Λ. Here, Λ, the scale of the theory, is of the order a few TeV and we choose 0 < (k q , k g ) < 1. For this analysis we use PDF4LHC15 [74] parton distribution functions (PDF) throughout from LHAPDF [75] subroutine unless otherwise stated. For the fixed order as well as resummed computation, we convolute NLO(NNLO) level partonic coefficient functions with NLO(NNLO) PDF. The corresponding strong coupling constant α s (µ 2 r ) is also provided by LHAPDF subroutine and for convenience we have defined a s (µ 2 r ) = α s (µ 2 r )/(4π). The fine structure constant is taken to be α em = 1/128 and the weak mixing angle is sin 2 θ w = 0.22343. Here we present the results for n f = 5 flavors in the massless limit of quarks. Our default choice for the center of mass energy of the LHC is E CM = 13 TeV. Except for the study of scale variations, we set renormalization (µ r ) and factorization (µ f ) scales equal to the dilepton invariant mass, i.e. µ r = µ f = Q. Our default choice of the cut-off scale of the theory is Λ = 3 TeV.
Before we present the resum result, we discuss some distinctive features of this model that are not observed for the case of universal couplings. For this, in fig. (1) left panel we present the dilepton invariant mass distribution at the resonance Q = M G for different choices of the model parameters k q and k g . For the results presented just in fig. (1), we use NNLO PDF at all orders. We observe that even at the resonance region, the cross section depends on the choice of k q and k g , which is not the case for the RS model where the height of the peak is independent of the corresponding coupling c 0 [44]. This particular nature can be understood from the parton level Born cross sections, eq. (3.2). In fig. (1) right panel we present the corresponding K-factors (K 1 , K 2 ) defined with respect to LO as, Here we observe that for (k q , k g ) = (1.0, 0.1) at low Q = M G values, the NLO K-factors K 1 can be as high as 3 while the corresponding K-factors for the universal case are about 1.5 [44]. Because k q is much larger than k g , at LO one can expect dominant contributions from Drell-Yan like process while the Higgs-like process is suppressed because of smaller coupling k g . For the same reason, at NLO the gluon emissions from the underlying born level parton processes is also expected to give smaller contribution. However, at NLO the presence of additional subprocesses like qg with large parton fluxes, give a significant model dependent contributions at NLO, it is necessary to include the QCD corrections at second order and beyond in the perturbation series. We present the NNLO K-factors (K 2 ) in the right panel of fig. (1) and it can be seen that these second order contributions are smaller than the corresponding NLO ones for different extreme choices of the non-universal couplings and confirms that at NLO all possible dominant contributions are considered giving typical sizeable corrections from NNLO onwards. In this regard, it is convenient to study the higher order QCD corrections in terms of the K-factors R nm defined with respect to NLO. In the rest of the analysis, we present the results in terms of these K-factors R nm . Next, we present the dependence of model parameters k q and k g on the FO cross sections at the resonance for Q = M G = 1 TeV through the contour plot in fig. (2). The cross section is large for larger k q values and becomes maximum around the universal line (k q ∼ k g ). The behaviour can be better understood from the dependence of born level partonic subprocesses on k q and k g as, (3.2) Here the coefficients A and B contain the contributions for the decay of the spin-2 particle to fermions and bosons respectively and C(Q 2 ) = (Q 2 −M 2 G ) 2 . At resonance where C(Q 2 ) = 0, the cross section will increase with increasing k q for a fixed k g and k g effect is mild here. For any given k q , the cross sections are maximum when k g k q .
We now discuss the resummed effects for the non-universal couplings. We have performed resummation up to NNLL accuracy and matched them with the fixed order NNLO level as described in sec. (2). In fig. (3), we have shown the dependence of the cross section on k q and k g including the resummation effect at NNLL. In the left panel, we present the  contour region at the resonance Q = M G = 1 TeV and in the right panel we present the same but at the off-resonance region Q = 1.5 TeV. Note that unlike the RS case, here we find non-negligible GR contribution away from the resonance as well. The dependence on the model parameters is found to be the similar as that observed in the FO case fig. (2). However, the maximum cross section region for the resummed results is wider. We have also studied the cross-section regions for varying the spin-2 mass M G and either of the couplings while keeping the other coupling constant. In fig. (4), we present such behaviour at the resonance region Q = M G to NNLO+NNLL accuracy for different values of M G and k q for a given k g = 1.0 in the left panel. We observe that with increasing M G the cross  section is falling for a fixed k q . This could be described due to the decreasing parton fluxes. The cross section is increasing with k q for a fixed value of M G as discussed in eq. (3.2). In the right panel, we present similar results but by varying M G and k g for a given value of k q = 1.0 and we observe the mild effect of k g on cross section when both M G and k q are fixed. To quantify the resummation effect we define K-factor with respect to NLO as discussed before. We define, We present in fig. (5), the ratio of resum results matched to the FO ones as given in eq. (3.3). We notice that there is around 4% enhancement due to NLL resummation over the NLO FO result at the resonance region while the enhancement due to NNLL resummation is around 2.5% over the NNLO results. In general we notice that in the region Q > M G , the resummation effects keep increasing with Q and are dominant at NLO+NLL level, whereas they are almost constant and are about 1% at NNLO+NNLL compared to NNLO. In fig. (6), we present the invariant mass distribution to NNLO+NNLL accuracy in the left panel and the corresponding K-factors in the right panel for the default choice of model parameters. We observe a significant enhancement in the cross section from lower order to higher order. In order to study different contributions coming from SM, pure GR and the signal, we present fig. (7). We observe that the contribution of gravity in the invariant mass distribution is negligible for Q < 900 GeV. However, after the resonance the gravity contribution is also prominent which is in contrast to RS scenario as mentioned earlier.
In the right panel we present the corresponding K-factors for the SM, pure gravity and the signal. At low Q value the signal K-factor is equal to that of SM. At the resonance, most of the signal contributions are coming from the gravity, therefore, the signal K-factor is almost equal to that of gravity. In fig. (8), we show the NNLO+NNLL K-factors for different choices of couplings. It can be seen that the higher cross-section is achieved when k g k q , We also estimate various theoretical uncertainties in our analysis. We first consider the uncertainties due to the unphysical scales µ r and µ f . To quantify these uncertainties, we use  Figure 8. K-factor of signal with respect to NLO at NNLO+NNLL level for different choice of k q and k g the conventional canonical 7-points scale variations by varying µ r and µ f simultaneously from Q/2 to 2Q subject to the constraint that the ratio of unphysical scales should not be than 2 and taking the maximum absolute deviations. We put the following constraints on the variation in order to remove extreme combinations, In fig. (9), we present these scale uncertainties both in the fixed order results (left panel) as well as in the resummed results (right panel). At NNLO this uncertainty is estimated to be around ±0.6% for Q < M G , around ±6.3% at the resonance region Q = M G and is about ±2.0% for Q > M G . However, the corresponding scale uncertainties for the signal at NNLO+NNLL accuracy found to get significantly reduced, respectively, to around ±0.4%, ±3.6% and ±1.0%. We have also estimated the uncertainties only due to the remormalization scale µ r by varying it from Q/2 to 2Q and keeping µ f = Q fixed. We observe a significant reduction in the renormalization scale uncertainty from ±1.4% at NNLO to about ±0.5% at NNLO+NNLL accuracy. The corresponding factorization scale uncertainties obtained by varying µ f from Q/2 to 2Q and keeping µ r = Q fixed are found to get reduced from 5% at NNLO to 2.6% at NNLO+NNLL accuracy at the resonance. Apart from the unphysical scale uncertainties, we also estimate the uncertainties in the fixed order NNLO result as well as the resummed NNLO+NNLL predictions coming from non-perturbative PDFs. We have estimated the intrinsic PDF uncertainty for the default choice of PDF4LHC15 PDFs using the recommendation in [76], and present these results in tab. (1). We observe that the PDF uncertainty is increasing with invariant mass (Q). However, there is no significant improvement in the PDF uncertainty after inclusion of threshold logarithms. Furthermore, we also compute the cross sections at NNLO+NNLL accuracy for the central set i = 0 of different PDF groups namely, MMHT2014nnlo68cl [77], CT14nnlo [78], ABMP16 5 nnlo [79] and NNPDF31 nnlo as 0118 [80] at resonance Q = M G = 1 TeV and at Q = 1500 GeV. The corresponding results obtained are found to be different (as much as 3%) in some cases than the our default choice PDF set. Although the main focus of our phenomenological study is on the threshold resummation, we have also studied the soft plus virtual (SV) corrections at third order in QCD, N 3 LO SV . We have computed these third order SV coefficients using the three-loop forms factors [47] and exploiting the universality of soft radiations. The new results are presented in appendix C. We finally then give a numerical estimate of these third order corrections by using the running strong coupling constant at 4-loop [81] level. We use the same PDF4LHC15 NNLO set at this order and find that the third order SV result contributes an additional 1% to the NNLO result at the resonance region. The renormalization scale uncertainty reduces to ±0.2% at resonance for a canonical variation 5 within µ r ∈ {1/2, 2}Q keeping µ f = Q fixed. However for a proper estimation of third order QCD corrections (particularly in the region away from threshold region viz. in the lower invariant mass region) one needs to include the regular terms pieces as well as PDF at the third order.

Summary
We have studied a detailed phenomenology for a generic spin-2 particle production at the 13 TeV LHC. We assumed non-universal coupling of the spin-2 particle with the SM fields. This needs additional UV renormalization for individual operators. After performing UV renormalization of the form factors, we have computed the new SV coefficients at the third order. From these coefficients, we were able to extract the process dependent coefficients needed for threshold resummation. Using the universal threshold exponent, which are already available in the literature as well as the newly computed process dependent coefficients in this article, we find all the ingredients to perform resummation up to NNLL accuracy. We observe a better perturbative convergence after inclusion of these threshold logarithms. Compared to the fixed order, the cross-section increases by 2.5% at the NNLL level at the resonance. We show that inclusion of these threshold logarithms are indeed important in taming the theoretical uncertainties to as small as 3.6% near the resonance. We also discuss the impact of the third order SV coefficients and the corresponding renormalization scale uncertainty. We stress that the K factors in the non-universal case strongly depend on the higher order corrections and a naive scaling from LO will heavily undermine 5 Note that it is also possible to estimate theoretical scale uncertainty through a probabilistic description as in [82,83]. However we refrain from this study in this article.
the correct result at NLO(NNLO) level. Our mass-dependent K factors are thus expected to be useful in the search of such spin-2 resonances at the LHC.

A Anomalous Dimensions
The UV anomalous dimensions appearing in Eq. (2.10) can be written in terms of renormalized strong coupling as Up to three-loop the anomalous dimensions are extracted by imposing the universal IR structure of on-shell Form factor in [47] and also collected below for completeness, Note that the conservation of the sum of the operators fixes the remaining anomalous dimensions i.e.

B Process dependent resum coefficients
The process-dependent resum coefficients are collected here for both quark and gluon channels with exact dependence of renormalization and factorization scales. Notice that the coefficients now explicitly depend on the gluon and quark coupling to spin-2 particle through the ratio r q = k g /k q for quark initiated subprocess and through r g = k q /k g for gluon initiated subprocess. At each order the coefficients take the following form, Notice that the universal case is recovered after realising r p = 1. Defining L qr = ln(Q 2 /µ 2 r ) and L f r = ln(µ 2 f /µ 2 r ), Note that all the coefficients for the SM DY process up to third order have been obtained by some of us previously in [41].

C Soft-Virtual coefficients at the third order
Here we collect the new third order soft-virtual coefficients for both gluon and quark channels.