Approximate NNLO Predictions for the Stop-Pair Production Cross Section at the LHC

If the minimal supersymmetric standard model at scales of around 1 TeV is realized in nature, the total top-squark pair production cross section should be measurable at the CERN Large Hadron Collider. In this work we present precise predictions for this observable, which are based upon approximate NNLO formulas obtained using soft-collinear effective theory methods.


Introduction
One of the main goals of the physics program at the Large Hadron Collider (LHC) is to investigate the existence of supersymmetric partners of the fundamental particles with masses in the TeV range. Within the context of the minimal supersymmetric standard model (MSSM) with R-parity conservation, supersymmetric particles are produced in pairs. At hadron colliders, the supersymmetric particles which are expected to be most abundantly produced are the ones which carry color charge: squarks and gluinos. Since supersymmetry is broken, the mass spectrum of squarks and gluinos plays a crucial part in determining which among the colored supersymmetric partners is experimentally accessible. Within the context of unified supersymmetric theories, the scalar and gaugino masses are evolved from a common high scale down to the energy scale of electroweak symmetry breaking. As a consequence of large Yukawa and soft couplings entering the evolution of the mass parameters, the third generation of squarks can have very different masses compared to the first two generations of squarks. In particular, under the assumption of "natural supersymmetry" [1], top-squarks are expected to be relatively light in order to reproduce the correct energy scale of electroweak symmetry breaking. Therefore, the lightest of the two supersymmetric partners of the top quark is usually expected to be the lightest squark in the mass spectrum. Precise theoretical predictions of the stop pair production cross section are instrumental in setting a lower bound on the lightest stop mass. Moreover, if top squarks will be discovered, accurate predictions of the cross section for stop pair production can be employed to determine the masses and other properties of these particles. For these reasons, the study of the radiative corrections to the production of stop pairs has already a quite long history. The calculation of the next-to-leading order (NLO) corrections to the cross section for stop pair production within the context of supersymmetric quantum chromodynamics (SUSY-QCD) was completed 15 years ago [2]. As expected, it was found that the NLO corrections significantly decrease the renormalizationand factorization-scale dependences of the prediction when compared to the leading order (LO) calculation. Furthermore, NLO SUSY-QCD corrections increase the value of the cross section if the renormalization and factorization scales are chosen close to the value of the stop mass. The NLO SUSY-QCD corrections are implemented in the computer programs Prospino and Prospino2 [3]. The electroweak corrections to stop pair production were studied in [4,5]. While these corrections have a quite sizable effect on the tails of the invariant-mass and transverse-momentum distributions, they only have a moderate impact on the total cross section. The emission of soft gluons accounts for a significant portion of the NLO SUSY-QCD corrections [6], which are large. For this reason, the resummation of the next-to-leading logarithmic (NLL) corrections was carried out in [7]. It was found that these corrections increase the cross section at the LHC by up to 10% of its NLO value, while they further decrease the scale dependence of the prediction. In [8][9][10] next-to-next-to-leading order (NNLO) threshold corrections and Coulomb corrections were derived by means of resummation techniques. In these studies, the resummation is carried out in Mellin moment space.
In the last few years, a formalism based on soft-collinear effective theory (SCET), which allows one to resum soft-gluon corrections directly in momentum space, was developed [11,12] and applied to QCD corrections for several processes of interest in collider phenomenology, such as Drell-Yan scattering [13], Higgs production [14,15], direct photon production [16], and recently slepton pair production [17]. A similar approach was developed independently in [18], where methods of SCET and non-relativistic QCD were used to resum simultaneously the effects of soft and Coulomb gluons. This method was applied to squark and gluino pair production [18,19], where soft and Coulomb-gluon contributions were resummed at NLL order, and to top-quark pair production [20], where resummation up to next-to-next-to-leading logarithmic (NNLL) order was implemented. The soft-gluon corrections to the top-quark pair production cross section at the Tevatron and the LHC were also studied in [21,22], within the framework developed in [11,12]. By employing this formalism, it was possible to resum soft-gluon emission corrections up to NNLL accuracy, as well as to derive approximate NNLO formulas for the total production cross section, the pair invariant-mass distribution, the topquark transverse-momentum spectrum, and the top-quark rapidity distribution.
In fixed-order perturbation theory, differential partonic cross sections involve both singular distributions and regular terms. The singular distributions are functions of a "soft parameter", whose precise definition depends on the kinematics. In [13,21], where pair-invariant mass (PIM) kinematics was employed, the relevant soft parameter is √ s(1 − z), where √ s is the partonic center of mass energy, z = M 2 /s, and M is the pair invariant mass of the energetic particles produced. In this case, the soft (or partonic-threshold) region corresponds to the limit z → 1. In [22], where one-particle inclusive (1PI) kinematics was employed, the relevant soft parameter is s 4 , which is obtained by subtracting the heavy-particle mass squared from the invariant mass of the objects recoiling against the observed heavy particle. In this case, the soft gluon emission region is identified by taking the limit s 4 → 0. Detailed studies of Drell-Yan and top pair production processes showed that, after the convolution of the hardscattering kernels with the parton luminosities, the terms which are singular in the partonicthreshold limit provide the numerically dominant contributions to the hadronic cross sections. In these processes, the singular terms typically account for more than 90% of the total NLO corrections. Furthermore, the dominance of the partonic-threshold regions in the calculation of the cross sections arises dynamically, since it is due to the strong fall-off of the parton luminosities outside the threshold region. This phenomenon is called "dynamical threshold enhancement" [13].
In the partonic-threshold region, the partonic cross sections factorize into the product of hard functions, which contain the virtual corrections to the cross section, and soft functions, which account for the effects of soft gluon emissions. The singular (plus) distributions in the partonic cross section can be obtained by solving the renormalization group equations (RGEs) satisfied by the hard and soft functions. The corresponding anomalous dimensions are known up to NNLO [23,24]. The information extracted from the RGE and the knowledge of the hard and soft functions at NLO allows us to obtain approximate formulas including all of the singular terms up to NNLO. These approximate NNLO formulas can then be matched with exact NLO calculations in order to obtain precise theoretical predictions for the observables of interest.
Within this context, the production of top-squark pairs in the soft-gluon emission limit can be studied in analogy to the production of top-quark pairs. The hard and soft functions are matrices in color space. The soft functions in the stop-quark and top-quark production processes are identical; they were evaluated in [21,22] up to NLO. In contrast, the hard functions for the stop-pair production process are so far unknown and must be computed with NLO accuracy. It is important to observe that the only SUSY parameter which appears in the soft functions is the mass of the produced stop quark, while all of the other SUSY parameters appear exclusively in the hard functions. In this work, we carry out the calculation of the hard functions for top-squark pair production up to NLO. By combining the NLO hard functions with the NLO soft functions and with the anomalous dimensions entering the RGEs satisfied by the various terms in the factorized cross section, it is possible to resum soft-gluon emission corrections up to NNLL order. Here we limit ourselves to re-expand the resummed formulas in order to obtain approximate NNLO formulas for the pair invariant-mass spectrum and the stop transverse-momentum and rapidity distributions. Although our formulas enable us to obtain predictions for these differential distributions, we focus our attention on the total topsquark production cross section, which is the observable of most immediate interest in the top squark searches. In fact, by integrating the approximate NNLO formulas for the differential distributions over the complete phase-space, we obtain predictions for the total top-squark pair production cross section at the LHC, and we comment on the phenomenological impact of the NNLO corrections arising from soft emissions. In a future work, we will evaluate the resummed stop pair production cross section and discuss the phenomenological impact of soft gluon resummation at NNLL accuracy. It must be remarked that the resummation at NNLL accuracy, in Mellin space, is already available for the production of other SUSY particles; the production of squark pairs was carried out in [25], while the production of gluino pairs was studied at this level of accuracy in [26]. Furthermore the NLO hard matching coefficients for squark and gluino hadroproduction, an important ingredient for NNLL studies, were evaluated very recently [27].
The paper is organized as follows: In Section 2 we introduce our notation and conventions, which are very similar to the ones employed in [21,22]. Furthermore, we describe the factorization of the stop pair production cross section in the soft limit, both in PIM kinematics and in 1PI kinematics. In Section 3, we discuss the calculation of the soft and hard functions up to NLO. In Section 4, we discuss the structure of the hard scattering kernels and we present approximate NNLO formulas for the stop pair production process. Predictions for the total top-squark pair production cross section at the LHC can be found in Section 5, together with an analysis of the phenomenological impact of the approximate NNLO corrections on this observable. We conclude Section 5 by comparing our results with the NLO+NLL results of and [7] and [19], which are obtained by means of techniques and kinematic schemes different from the ones that we employ in this work. Finally, we collect our conclusions in Section 6.

Notation
In this paper we study the process where N 1 and N 2 indicate the incoming protons in the case of a proton-proton collider as the LHC, while X is an inclusive hadronic final state. In the rest of the paper, we treat the top squarks as on-shell particles and neglect their decay. The terms neglected in this approximation are of order Γt 1 /mt 1 . At the lowest order in perturbation theory, two partonic  Figure 1: Tree level diagrams contributing to the quark annihilation production channel (first row) and gluon fusion channel (second row). channels contribute to the process in Eq. (1): where the momenta of the incoming partons are related to the momenta of the incoming hadrons by p i = x i P i (i = 1, 2). The diagrams contributing to the two production channels at lowest order in QCD are shown in Figure 1. The relevant invariants for the hadronic scattering are In order to describe the partonic scattering, we employ the invariants In Born approximation s + t 1 + u 1 = 0, and consequently M 2 = s and s 4 = 0. It is well known that the kinematics of the process allows one to define different threshold regions. Here we consider two different cases: PIM kinematics, in which the threshold region is defined by the limit s → M 2 , and 1PI kinematics, in which the threshold region is approached by taking the limit s 4 → 0. Both regions were employed in the study of the top-quark pair production cross section and differential distributions [21,22,28,29]. In particular, by working in PIM kinematics one can calculate the pair invariant-mass distribution, while by working in 1PI kinematics one can evaluate the stop transverse-momentum and rapidity distributions. It should be emphasized that in the PIM and 1PI threshold regions the top squarks are not forced to be nearly at rest, as in case of the threshold region defined by the limit β = 1 − 4mt 1 /s → 0, which is often employed in the calculation of soft-gluon corrections to the total cross section [7,18,19]. In the rest of the paper, we refer to the β → 0 limit as the production threshold region.
Our goal is to employ both PIM and 1PI kinematics to obtain approximate NNLO formulas for the total top-squark pair production cross section. Both approaches include the numerically large contributions arising from the emission of soft gluons. Differences between them, and with the production threshold calculations, arise from the way in which different sets of powersuppressed corrections are treated.

PIM kinematics
We focus first on the case of PIM kinematics. It is convenient to introduce the quantities Consequently, the PIM threshold limit s → M 2 corresponds to the limit z → 1. According to the QCD factorization theorem [30], the differential cross section in M and θ (the scattering angle of the top squark with respect to the beam axis in the partonic rest frame) is given by where µ f is the factorization scale, and the sum runs over the incoming partons. In the following we drop the subscript PIM (and the corresponding subscript 1PI) whenever there is no ambiguity about the kinematic scheme employed, or when a formula applies to both schemes. The parton luminosities ff ij are defined as the convolutions of the non-perturbative parton distribution functions (PDF) for the incoming partons i and j: The functions C ij in Eq. (6) are the hard-scattering kernels, which are related to the partonic cross sections and can be calculated in perturbation theory. In order not to clutter the notation, we do not indicate explicitly the fact that the hard-scattering kernels depend on the top-squark masses mt 1 and mt 2 , the mass mq of the first two generations of squarks and of the sbottoms (which we assume to be all degenerate), the top-quark mass m t , the gluino mass mg, and thẽ t 1 -t 2 mixing angle α. The expansion of the C ij functions in powers of α s has the generic form Only the quark annihilation and gluon fusion channels contribute to C ij at lowest order in perturbation theory; in particular where N = 3, and the Mandelstam invariants t 1 and u 1 can be written in terms of s and θ as In order to calculate higher-order corrections to C qq and C gg one needs to consider virtual and real-emission corrections to the Born approximation. Starting at order α s , new production channels, such as qg →t 1t * 1 q, open up. When working in the threshold limit z → 1, the calculations are simplified by the fact that there is no phase-space available for the emission of additional (hard) partons in the final state. Consequently, both the hard-gluon emission and the additional production channels are suppressed by powers of (1 − z) and can be safely neglected while deriving results within the partonic-threshold limit. By neglecting powersuppressed terms in the integrand, Eq. (6) can be rewritten as In Eq. (11) the quark channel luminosities ff qq and ffq q are understood to be summed over all light quark flavors. The two terms in the second line of Eq. (11) differ in the fact that in the first term the quark (antiquark) comes from the hadron N 1 (N 2 ) in Eq. (1), while in the second term the quark (antiquark) comes from the hadron N 2 (N 1 ), respectively. The total cross section can be obtained by integrating over cos θ in the range [−1, 1] and over M in the range [2mt 1 , √ S]. In the soft-gluon emission limit z → 1, the hard-scattering kernels C ij factor into a product of hard and soft functions: Here and in what follows we employ boldface fonts to indicate matrices in color space, such as the hard functions H ij and the soft functions S ij . For simplicity we drop the the top mass and the SUSY parameters from the list of arguments of the hard functions H ij , as well as the stop mass from the arguments of the soft functions S ij . Throughout this paper, we work in the s-channel singlet-octet basis where a i represent the color index of the particle with momentum p i . We view these structures as basis vectors |c I in the space of color-singlet amplitudes. Inner products in this space are defined through a summation over color indices as This inner product is proportional but not equal to δ IJ , so the basis vectors are orthogonal but not orthonormal. A factorization formula analogous to Eq. (12) for the top-quark pair production was derived employing SCET and heavy-quark effective theory in [21]. A completely analogous procedure can be followed in order to derive Eq. (12), which is valid in the case of top-squark pair production.
The hard functions are obtained from the virtual corrections and are ordinary functions of their arguments. The soft functions arise from the real emission of soft gluons and contain distributions which are singular in the z → 1 limit. Contributions of order α n s to the soft functions include terms proportional to the plus distributions as well as terms proportional to δ(1 − z). In particular, the NLO and NNLO hard-scattering kernels in Eq. (8) have the structure The functions D ij can in principle be obtained from results present in the literature. One of the main results of this paper is the calculation of the coefficients D (2) m,ij (with m = 0, . . . , 3) both in the quark annihilation and gluon fusion channels. We can also evaluate all of the scale-dependent terms in Q (2) 0,ij in both channels, but due to the ambiguity on the choice of the normalization scale in the argument of these logarithms we drop part of these terms in the numerical implementation of our formulas. We will return to this issue below.

1PI kinematics
The 1PI kinematics approach allows one to describe observables in which a single particle, rather than a pair, is considered. One can then write the top-squark rapidity (y) and transverse-momentum (p T ) distribution as once again we will drop the subscript 1PI in most of our equations below. The expansion of the 1PI hard-scattering kernels C 1PI in powers of α s has the same structure shown in Eq. (8) for the PIM case. Obviously, also in this case only the qq and gg channels give non-vanishing contributions at lowest order in α s . The hadronic Mandelstam variables T 1 and U 1 are related to the stop rapidity and transverse momentum through the relations where m ⊥ = p 2 T + m 2 t 1 . Therefore, the variables s, s 4 , t 1 , u 1 , which are arguments of the 1PI hard functions can be expressed in terms of p T , y, x 1 , x 2 . The lower integration limits in Eq. (17) are In order to obtain the total cross section, it is necessary to integrate the double-differential distribution with respect to the top-squark rapidity and transverse momentum over the range In the case of 1PI kinematics, the hard-scattering kernels in the soft-emission limit s 4 → 0 factor into a product of hard and soft functions, in analogy with Eq. (12): As emphasized in [22], the Mandelstam invariants s , t 1 , u 1 can differ from s, t 1 , u 1 by power corrections proportional to s 4 . For example, explicit results for the hard and soft functions can be rewritten by employing either the relation s + t 1 + u 1 = 0 or s + t 1 + u 1 = s 4 . The difference between the two choices is due to terms suppressed by positive powers of s 4 . A detailed description of the way in which we deal with this ambiguity can be found in Section 4 of [22]. As in the case of PIM kinematics, the hard and soft function are matrices in color space originating from virtual and soft-emission corrections, respectively. The 1PI hard functions are identical to the ones encountered in the case of PIM kinematics, provided that the variables s, t 1 , and u 1 are written in terms of M and cos θ. The soft functions are different in the PIM and 1PI schemes, but in both cases they are identical to the ones employed in the calculation of the top-quark pair production cross sections in [21,22].
The 1PI soft functions at order α n s depend on the associated plus distributions It follows that where g is a smooth test function. At NLO and NNLO, the hard-scattering kernel in 1PI kinematics have the structure As in the PIM case, the NLO coefficients D 0,ij , and R (1) ij can be in principle obtained from the literature. In this work we are able to derive exact expressions for the NNLO coefficients D

The hard and soft functions at NLO
In this section we describe the calculation of the soft and hard matrices up to NLO in perturbation theory.

Hard functions
As it was shown in detail in [21], the hard functions are defined in terms of the products of Wilson coefficients with their complex conjugates. In order to obtain the Wilson coefficients, one matches renormalized Green's functions in SUSY-QCD with those in SCET. The matching procedure can be carried out by choosing arbitrary external states and infrared (IR) regulators. The simplest and most common procedure consists in employing on-shell states for the partonic processes (qq, gg) →t 1t * 1 and in using dimensional regularization to regulate both IR and ultraviolet (UV) divergences. With this choice, the SCET loop diagrams vanish because they are scaleless. Consequently, the effective theory matrix elements are equal to their tree-level expressions multiplied by a UV renormalization matrix Z. The SUSY-QCD matrix elements are instead the virtual corrections to the scattering amplitudes of the partonic processes (qq, gg) →t 1t * 1 . Rather than implementing the matching condition for the Wilson coefficients, we directly move to the calculation of the the matrix elements H IJ ij of the hard functions projected onto a certain basis {|c I } of color structures. These matrix elements, and not the Wilson coefficients, are needed in order to obtain the approximate NNLO formulas which represent the main result of the present work. Throughout this paper, we use the s-channel singlet-octet basis defined in Eq. (13). In order to calculate the hard-function matrix elements, we use the fact that the renormalized hard functions can be directly obtained from the corresponding (squared) on-shell QCD scattering amplitudes. The infrared poles can be removed from the QCD amplitudes by employing the prescription of [31,32] M ren where the indices ij are not summed over. Here |M ij ( ) are the dimensionally regularized (and UV renormalized) scattering amplitudes, whose IR poles are removed by the factors The resulting finite amplitudes are expressed in terms of α s with n l = 5 active flavors. The explicit form of the Z ij matrices can be obtained by employing the results derived in [23,24,33].
The perturbative expansion of the renormalized hard functions is defined as where d R = N in the quark annihilation channel and d R = N 2 − 1 in the gluon fusion channel. The desired matrix elements can then be written in terms of the renormalized QCD amplitudes and the color basis vectors |c I as The leading-order result for the qq channel follows from a simple calculation. In matrix notation, it reads while that for the gg channel is In order to calculate the NLO hard matrices H (1) ij , one needs to evaluate the one-loop corrections to the partonic scattering amplitudes, by keeping the various color components separate. Although results for the corresponding one-loop diagrams interfered with the Bornlevel amplitudes are well known [2], their decomposition into the color basis is not available in the literature, and we therefore had to calculate it from scratch. For this purpose, we have used in-house routines written in the computer algebra system FORM [34] in combination with Reduze [35,36]. The results of the calculation are expressed in terms of Passarino-Veltman functions [37]. After UV renormalization, we have derived analytic expressions for the IR poles, and we have evaluated the finite parts of the amplitudes numerically by employing the programs described in [38][39][40]. We have checked our results in several ways. First, by applying the IR renormalization factors Z ij of Eq. (25), we find that the IR poles cancel exactly. Second, we have checked that by multiplying the one-loop hard functions with the corresponding tree-level soft functions and by subsequently taking the trace of the resulting color-space matrices, we reproduce the numerical results for the NLO virtual corrections which can be extracted from the code Prospino [3].
The output of our FORM codes are too long to be presented here in explicit form. The explicit expressions for the NLO hard functions are coded in a Mathematica and a Fortran program which are included in the arXiv version of this work. In particular, these programs allow one to evaluate numerically the LO and NLO hard functions for an arbitrary choice of the input parameters.

Soft functions
The soft functions are vacuum expectation values of soft Wilson-loop operators. These functions are not sensitive to the spin of the particles involved, but they depend on the color structure of the underlying partonic subprocesses. Consequently, the soft functions needed for the calculation of the top-squark pair-production cross section are precisely the same functions employed in the calculation of the cross section for top-quark pair production. The calculation of the PIM soft functions at NLO was described in detail in [21], while the analogous calculation of the 1PI soft functions was carried out in [22]. For the convenience of the reader, the explicit results for the soft functions obtained in those two papers are collected in Appendix A.
4 Structure of the hard scattering kernels 4

.1 Preliminaries
Since the soft functions in (12) and (21) depend on plus distributions, it is more convenient to work with the Laplace-transformed functions in PIM and similarly in the case of 1PI kinematics. It was shown in [11][12][13]29] that the leading singular terms in the hard-scattering kernels can be generated by replacing the Laplace variable L with a derivative ∂ η with respect to an auxiliary variable η, which is later set to 0. To this end, one defines Laplace-transformed hard-scattering kernels in PIM kinematics as For 1PI kinematics one replaces the arguments M and cos θ with s , t 1 , and u 1 . The hardscattering kernel in momentum space can then be recovered through the relations [21,22,29] In order to evaluate the above formulas, one needs to employ analytic continuation to regulate the divergences in z → 1 or s 4 → 0, to take the derivatives with respect to η, and finally to take the limit η → 0. It is possible to show that the result of this procedure can equivalently be obtained by implementing a series of replacement rules on the functionsc ij (L, . . . ) considered as polynomials in a new variable L. In PIM kinematics, one must substitute where L M = ln (M 2 /µ 2 ), and the plus distributions P m (z) have been defined in Eq. (15). In 1PI kinematics, one must substitute instead Since the hard and soft functions are known up to NLO, is easy to determine the NLO coefficient in the expansion ofc in powers of α s : It is important to observe that the trace of the product of the LO hard function and NLO soft function contains the dependence ofc (1) ij on L, and therefore it gives rise to the plus distributions. In order to obtain the complete coefficientc (2) ij one needs to know the hard and soft functions up to NNLO: The NLO hard and soft functions in this equation are know exactly. As explained in detail in [21,22], the scale-dependent part of the NNLO hard function and the L dependent part of the NNLO soft function can be reconstructed by exploiting the information coming from the RGE satisfied by these functions. To this end, one only needs the one-loop hard and soft matrices calculated in Section 3 and in [21,22], respectively, as well as the relevant two-and three-loop anomalous dimensions computed in [23,24,41,42]. After this information on the NNLO hard and soft function has been extracted, one can evaluate the coefficients c

Approximate NNLO formulas
By employing the results described in the previous sections, we are able to obtain approximate NNLO formulas for the hard-scattering kernels in Eqs. (12) and (21). These formulas include the exact expressions of the coefficients multiplying the plus distributions up to NNLO, both in PIM and in 1PI kinematics. The complete results for these coefficients, written in terms of Passarino-Veltman functions, are too lengthy to be reported here. The values of the coefficients multiplying the various plus distributions and delta functions for arbitrary values of the input parameters can be extracted from the Mathematica code mentioned above, which we include in the arXiv submission of this work. 1 As a reference for the reader, we present here the hard scattering kernels evaluated by setting the input parameters at the values listed in Table 1. The SUSY spectrum which we chose corresponds to the benchmark point 40.2.5 in [43]; we employ this benchmark point for the remainder of this paper. By definingĈ (i) = d R C (i) (i = 0, 1, 2), in the quark annihilation channel with PIM kinematics one findŝ PIM,qq (z) = 27.0048 where the ellipses indicate terms which are subleading (and finite) in the z → 1 limit. In the gluon fusion channel, we find instead PIM,gg (z) = 16.7315 PIM,gg (z) = 401.555 In order to completely determine the coefficients multiplying the delta functions in the NNLO hard-scattering kernels, one would need to know the complete NNLO hard and soft matrices. Since those matrices are at present not fully known at NNLO, we are only able to determine the scale-dependent terms in the delta-function coefficients, because those terms are governed by RGEs. However, since the scale-independent parts are unknown, the coefficients of the delta function at NNLO depend on an arbitrary second scale chosen to normalize the scaledependent logarithms. In fact, for any renormalization scale µ and reference scale µ 0 , one can always rewrite where the second term in the r.h.s. can be reabsorbed in the unknown µ-independent piece. In Eq. (37) and Eq. (38), we chose the reference scale equal to the pair invariant mass M . Furthermore, in Eq. (37) and Eq. (38), we decided to separate out the contributions to the delta-function coefficients coming from the NNLO hard functions. In particular, in the round brackets multiplying the delta functions, the second number represents the contributions of the scale-dependent terms in H (2) ij . In our numerical analysis in Section 5, we decided to drop the contributions to the NNLO delta-function coefficients arising from the two-loop hard functions. As was observed in [22], this choice is motivated by the fact that by including these µ-dependent terms one might induce an artificial reduction of the scale dependence, which can lead to an underestimated scale uncertainty.
Similarly, in 1PI kinematics we find for the quark annihilation channel while for the gluon fusion channel we obtain In this context, the ellipses indicate subleading terms in the s 4 → 0 limit. Also, in Eqs. (40) and (41), as in the PIM case, the scale-dependent terms are normalized to the scale s, which is equal to M 2 when s 4 = 0. As in the PIM case, also in Eqs. (40) and (41) the second number in the round brackets multiplying the delta functions in the NNLO kernels represents the contribution of the scale-dependent terms arising from the NNLO hard functions.

Total cross section
In this section we present a brief numerical study of the top-squark total pair production cross section at approximate NNLO accuracy. We stress that our approximate NNLO predictions include the full NLO corrections, as well as the part of the NNLO corrections arising from soft gluon emission, obtained by means of the procedure outlined in the previous sections.
Following what was done in the study of the top-quark pair production total cross section, we employ PIM SCET and 1PI SCET kinematic schemes described in [44]. These schemes also include, on top of contributions which are singular in the soft limit, NNLO terms which are regular in the z → 1 (PIM SCET ) or s 4 → 0 (1PI SCET ) limits, and which naturally arise from the SCET formalism. These terms, which are part of the functions R (2) ij in Eqs. (16) and (24), do not represent the complete contribution to the NNLO cross section which is regular in the soft limit, since this quantity can only be obtained with a full calculation of this observable at NNLO accuracy. However, as shown in [21,22], the regular terms appearing in the PIM SCET and 1PI SCET kinematic approaches arise from the exact definition of the soft-gluon emission energy, and they improve the agreement between exact and approximate formulas at NLO.
Unless otherwise stated, we present results which are obtained by averaging the ones obtained in the two kinematic schemes that we consider. In analogy with what was done in [44] for the top-quark pair production cross section, we adopt a conservative approach and consider the difference between the predictions in the two kinematic schemes as an estimate of the uncertainty associated with the use of approximate NNLO formulas. This is justified by the fact that the two schemes neglect different power-suppressed terms, which are formally subleading but which can nevertheless have a noticeable numerical impact on the total cross section. To account for this uncertainty, the scale variation of the total cross section is obtained by setting the renormalization and factorization scales equal to each other, µ R = µ f = µ, and by varying this single scale between mt 1 /2 < µ < 2mt 1 . We then look at the difference between the largest and smallest values obtained. In summary, the central value and perturbative uncertainties for the combined results at approximate NNLO accuracy are determined by employing the definitions where the subscripts 1PI and PIM indicate that the corresponding quantities are evaluated in 1PI SCET and PIM SCET kinematics, respectively, including the full set of NLO corrections and the contribution of the NNLO terms present in the approximate formulas for that scheme. As will be shown later, the total cross section is strongly dependent on the mass of the produced particle, mt 1 . However, similarly to the slepton-pair production case [17], the dependence of the total cross section on other SUSY parameters is relatively small. In order to show that this is indeed the case, we fix the value of the stop mass equal to mt 1 = 1087.17 GeV, corresponding to the SUSY benchmark point 40.2.5 of [43] (which is defined in Table 1), and we provide predictions for the total cross section for three different sets of the remaining SUSY parameters. The first set of SUSY parameters coincides with the benchmark point   Table 2: Stop-pair production cross sections at the LHC with √ S = 14 TeV for three different sets of the SUSY parameters described in the text. The stop mass is fixed to mt 1 = 1087.17 GeV. The numbers are obtained by using MSTW2008 PDFs [45][46][47].
of SUSY parameters (which we refer to as "half") is mt 2 = 659.93 GeV, mq = 730.15 GeV, mg = 744.99 GeV, α = 34.2 • . The values of the total cross sections, with the relative perturbative uncertainty, for the three different SUSY parameters sets discussed above, can be found in Table 2. The table refers to a collider energy of 14 TeV and the PDFs employed are from the MSTW2008 global fit [45][46][47]. We observe that the numerical values for the cross section are very close to each other, in spite of the fact that the input SUSY parameters are considerably different in the three sets. This is true both at NLO and at approximate NNLO accuracy.
In all of the plots and tables discussed below, the SUSY parameters are set to the values corresponding to the benchmark point 40.2.5 in Table 1. This applies also to the lightest stop mass, with the exception of the figures in which we plot the cross section as a function of the stop mass or when a different choice of the stop mass is explicitly indicated.
All of the numbers and plots are obtained by means of an in-house Fortran code, in which the approximate NNLO formulas are implemented. The NLO calculations, which are one of the elements needed to obtain predictions at approximate NNLO accuracy, have been carried out by modifying the public version of Prospino [3].
As a first step, we compare the full NLO cross section with the approximate NLO cross section given by the leading singular terms. To be specific, in the approximate NLO formulas we include the coefficients D  16) and (24), as well as those terms in R (1) ij which naturally arise in the SCET approach. The purpose of this comparison is to establish to what extent the leading terms in the threshold approximation reproduce the full cross section, or, in other words, if the dynamical threshold enhancement of the soft emission region takes place. This comparison is shown in Figure 2, for the case of a hadronic center of mass energy of 8 TeV. The two rows in the figure refer to two different choices of the PDF set. NLO PDFs are employed in all of the four panels. One observes that the average of the approximate PIM and 1PI NLO formulas reproduces very well the band obtained by varying the factorization scale (which is set equal to the renormalization scale) in the full NLO result. Furthermore, the comparison in Figure 2 supports the fact that the integrals of the differential distributions in PIM and 1PI kinematics over the whole available phase space reproduce to a good accuracy the total cross section, at least for the range of stop masses of interest in this work. It is thus reasonable to expect that also the approximate NNLO formulas reproduce to a good extent the unknown full NNLO corrections. Another comparison between the scale dependence of the approximate NLO cross section and the full NLO cross section is shown in Figure 3, in which we plot separate curves for the cross section obtained in PIM SCET and 1PI SCET kinematics, and for the average of the two. Both for mt 1 = 500 GeV and for mt 1 = 1087.17 GeV, the scale dependence of the approximate cross sections in the range mt 1 /2 ≤ µ ≤ 2mt 1 reproduces to a good extent the scale variation of the full NLO cross section in the same range.
We now turn to the discussion of the approximate NNLO results, recalling first that in the approximate NNLO formulas we do not include the terms proportional to the Dirac delta functions δ(1 − z) (PIM SCET ) or δ(s 4 ) (1PI SCET ) arising from the NNLO hard functions, since the scale-independent parts of the NNLO hard functions are still unknown. As expected, the approximate NNLO predictions for the pair-production cross section show a smaller scale dependence than the NLO results for the same quantity. This is illustrated in Figure 4, where two different LHC center of mass energies (8 TeV and 14 TeV) and two different stop-quark masses (mt 1 = 500 GeV and mt 1 = 1087.17 GeV) are considered. The approximate NNLO line is obtained by plotting the scale dependence of the average between the PIM and 1PI cross sections. In order to show the effect of the approximate NNLO corrections on the scale dependence, both the NLO and approximate NNLO curves are plotted using MSTW2008 NLO PDF sets. While in Figure 4 we plot exclusively the cross section scale dependence, we remind the reader that the perturbative uncertainty employed in all tables as well as in all other figures is obtained by evaluating the second and third line of Eqs. (42). Consequently, it reflects both the scale uncertainty and a kinematic scheme uncertainty, which is associated to the different sets of non-singular terms neglected in 1PI SCET and PIM SCET kinematics.  [43].
A more precise assessment of the impact of the approximate NNLO corrections on the central value of the cross section and on the associated perturbative uncertainty can be obtained by comparing predictions for fixed values of the stop mass. This analysis is presented in Tables 3 and 4, 5 and 6, and 7 and 8, which refer to the LHC with center of mass energies √ S = 7 TeV, 8 TeV, and 14 TeV, respectively. In all tables we consider two different values of the lightest top-squark mass: i) the value associated to the benchmark point 40.2.5, mt 1 = 1087.17 GeV, and ii) a stop mass mt 1 = 500 GeV, close to the current lower bound for this particle as determined by searches at the LHC. In all tables, the first uncertainty is perturbative, determined as explained in Eq. (42), while the second uncertainty is obtained by scanning over the 90% CL sets of the corresponding PDFs and by taking into account the error on α s (m Z ). In particular, the PDF uncertainty is determined by seeing how the average of the PIM and 1PI predictions changes when evaluated with different members of the PDF   In all cases listed in the tables, the inclusion of the approximate NNLO corrections reduces the perturbative uncertainty, when expressed as a percentage of the central value, by more than a factor of 2 with respect to the corresponding NLO prediction. We can summarize the content of the tables as follows: the scale variation in the range mt 1 /2 ≤ µ ≤ 2mt 1 can increase the NLO central value up to + [11,16]% or lower it up to − [13,18]%. At approximate NNLO, the scale variation can increase the central value of the cross section up to +[2, 5]% or decrease it up to − [5,7]%. These considerations are valid both when one employs CT10 PDFs [48] or MSTW2008 PDFs.
In almost all cases illustrated in the tables, the PDF and α s uncertainty grows marginally in the approximate NNLO predictions with respect to the NLO predictions. Another way to look at the PDF and α s uncertainty is shown in Figure 5 Table 6: Stop-pair production cross section for two different values of mt 1 at the LHC with √ S = 8 TeV. The numbers are obtained by using CT10 PDFs.
that both bands become larger for large stop masses. The approximate NNLO band in the left panel is almost everywhere inside the NLO band, while in the right panel the approximate NNLO band is larger than the NLO band. However, the bands obtained by using CT10 PDFs remain larger than the ones obtained when using MSTW2008 PDFs. Tables 3 to 8 also include the values for the K factors at NLO and approximate NNLO, which are both normalized to the LO cross section, i.e.
The NLO K factors tend to be slightly larger when CT10 PDFs rather than MSTW2008 PDFs are employed (roughly 1.5 vs. 1.3), but they are not very sensitive to the collider center of mass energy or to the mass of the top squark.   Finally, Figure 6 shows the cross section as a function of the top-squark mass up to mt 1 = 2 TeV for the LHC at 7, 8 and 14 TeV center of mass energy. In all cases, the bands represent the residual perturbative uncertainties, obtained as explained above. In Figure 6 we employ CT10 PDFs. One sees that, for large values of the stop mass, the approximate NNLO band is below the NLO band at 7 and 8 TeV center of mass energy, while, for the LHC at 14 TeV, the approximate NNLO band overlaps with the lower part of the NLO band. For comparison, in Figure 7 we repeat the same analysis by employing MSTW2008 PDFs for the LHC at 8 TeV. In this case, for large stop masses the approximate NNLO scale uncertainty bands tend to be slightly above the NLO bands.
We conclude this Section by comparing our results with the ones presented in [7] and [19], which have NLO+NLL accuracy. In Tables 9 and 10 we show the results obtained for the input parameters employed in [7], which coincide with the SPS1a' benchmark scenario of [50], as well as for two different values of the stop mass. Table 9 refers to the LHC with 7 TeV center of mass energy while Table 10 refers to the LHC with a center of mass energy of 14 TeV. We checked that, as expected, the NLO results in [7] coincide with the ones obtained with our modified version of Prospino. The central values of our approximate NNLO predictions for the total cross section are in very good agreement with the ones obtained by means of NLO+NLL formulas in [7]. However, the perturbative uncertainties at approximate NNLO accuracy appear to be smaller that the NLO+NLL scale uncertainties.   [49]. We report only the perturbative uncertainties.
In [19], the authors perform a simultaneous resummation of the production threshold logarithms and Coulomb singularities at NLL accuracy, including bound-state effects. Table  7 in [19] reports result for the production of top squark pairs at the CMSSM benchmark point 40.2.4 [43]. Results for higher values of the stop quark mass can be found in the files provided with the arXiv submission of [19]. Coulomb resummation and bound state effects increase the cross section, but the largest effect in the NLL results of [19] is due to soft resummation. The authors of [19] employ a private version of the MSTW2008 NLO PDFs in their calculations. In Table 11 we compare our findings with the ones of [19]. One can immediately see that our cross section prediction (inclusive of perturbative uncertainty) has a sizable overlap with the predictions quoted in [19] (inclusive of perturbative uncertainty). In all cases, the central value of our predictions falls within the perturbative uncertainty of the NLL calculations. The central values of the approximate NNLO calculations are lower than the ones of [19], as it can be seen from the fact that the quantitiesK i ≡ σ i /σ NLO are lower than the corresponding quantities in [19]. In analyzing this fact it is necessary to take into account that our result does not include the additional enhancements due to Coulomb corrections,   bound state effects and resummation of scale-dependent logarithms in the hard function 3 . In this sense, it would be more natural to compare our results with the predictions labeled NLL s and NLL s,fixed in [19]; these two approaches giveK factors which are roughly 0.03-0.08 smaller than the ones corresponding to the NLL predictions (see Fig. 15 in [19]), and fully compatible with our results. Additionally, experience with top quark production indicates that NNLL resummation gives cross sections which are larger than the ones obtained from the corresponding approximate NNLO formulas [21,22].
However, at this point in our discussion, the stage is set to comment on the reciprocal advantages and disadvantages of the various methods employed to implement soft gluon resummation and/or approximate NNLO formulas. This point can be better illustrated in the case of top quark pair production, where full NNLO results recently became available [51,52]. The various calculational schemes adopted by different groups can be classified with respect to two aspects: i) The kind of kinematics employed, and ii) the space in which the resummation   [43]. In particular, we set m t = 172.5 GeV, mg = 1386 GeV, mq = mt 2 = 1358 GeV and cos α = 0.39 as in [19]. The factorization scale is set equal to mt 1 . For the approximate NNLO results we used MSTW2008 NLO PDFs. It must be observed that the numbers in [19] were obtained with a private version of the MSTW NLO PDFs.
is carried out. Many predictions for the total top quark pair cross section are obtained by employing the traditional production threshold kinematics identified by the β → 0 limit [20,53], while PIM kinematics and 1PI kinematics, which additionally allow one to calculate differential distributions of phenomenological interest in top quark physics, are employed in [21] and [22,54], respectively. One can then decide to work in Mellin moment space [53,54], or in directly in momentum space [20][21][22]. All of these schemes allow one to obtain NNLL or approximate NNLO accuracy and they differ in the kind of formally subleading terms which are neglected. Without guidance from complete NNLO calculations there are no definitive arguments to prefer a priori one scheme to another. What one can do is to see how a given approximation does at NLO, where exact results are known, along the lines we followed at the beginning of this section. For the top quark production, all of the predictions obtained in the various schemes are compatible with each other when perturbative uncertainties are taken into account [55]. It is only a posteriori, when the results of a full NNLO calculation are available, that one can see which of the various approximate NNLO/NNLL results is closest to the exact NNLO result. In the case of the top quark total cross section, our method, based on PIM and 1PI kinematics and momentum space, provides results which are slightly lower than, but compatible with, the results of the other groups. It is important to stress that there is not a clear pattern indicating, even a posteriori, a definite methodological bias: the approximate NNLO prediction which best approximates the full NNLO total cross section is based upon 1PI kinematics and Mellin space [54], while Mellin space and momentum space calculations in the production threshold limit provide very similar results [20,53]. A detailed analysis of the various approaches was carried out in Section 5.2 of [22]; in particular, it was shown that the approach based on SCET and employed here significantly improves the agreement between the results obtained in 1PI and PIM kinematics. Furthermore, it is possible to argue that the good numerical agreement between the exact NLO result at the LHC and the corresponding approximation based on production threshold expansion is somewhat fortuitous: In fact, in spite of the fact that the β → 0 expansion fails to reproduce the correct shape of the (unphysical) distribution dσ/dβ in the gluon fusion channel, the integral of the approximate distribution is very close to the full NLO result (see for example the discussion in [56]). On the other hand it is fair to say that the difference between the full NNLO calculation of the top production cross section and the results of [21,22,44] is an a posteriori indication that the neglected formally subleading terms in those works turned out to be numerically larger than the neglected subleading terms in the other approaches.
In conclusion, we believe that, in absence of a full NNLO calculation, approaches leading to different kinds of approximate NNLO /NNLL predictions are well worth exploring. Furthermore, numerical differences among predictions based on different approaches should be conservatively taken as a measure of the uncertainty associated to the neglected subleading corrections. Precisely for this reason, following the procedure of [44], we decided to base our predictions for the stop production cross section on an average of the PIM and 1PI kinematics calculations and we estimated the perturbative error not only by means of the scale variation, but also by considering the maximal difference between the predictions obtained in the two kinematic schemes, as explained in this section.

Conclusions
Supersymmetry is certainly one of the best motivated scenarios for physics beyond the standard model. If supersymmetry is broken just above the electroweak scale, the supersymmetric partners of the known particles are expected to have masses of the order of 1 TeV, and they could soon be observed at the LHC. For these reasons it is important to have precise theoretical predictions for the production cross sections of supersymmetric particles at hadron colliders. Top squarks are expected to be among the lightest colored supersymmetric particles, and consequently they might be the supersymmetric partners which are most easily accessible at the LHC. In this paper we have employed effective field theory techniques in order to improve existing NLO calculations for stop-pair production by evaluating higher-order perturbative corrections arising from soft-gluon emissions.
In particular, we have adopted a framework which allows one to resum large logarithmic corrections at NNLL accuracy in the stop pair production process. The resummation can be carried out in two different kinematic schemes, referred to as PIM and 1PI. In principle, this fact enables one to obtain predictions for both the pair invariant mass distribution and the top-squark p T and rapidity distributions, as well as for the total cross section, which is at the moment the observable of primary interest in phenomenological studies. Furthermore, by re-expanding the NNLL resummation formulas, it is possible to obtain predictions for the cross section which have approximate NNLO accuracy in fixed-order perturbation theory. The evaluation of these approximate NNLO formulas represents the main goal of this work, where we obtained analytic expressions for all of the coefficients multiplying the plus distributions in the variables (1 − z) and s 4 in the NNLO hard-scattering kernels for the double differential distributions in the two kinematic schemes considered. This was made possible by a complete computation of the hard-function matrices in both production channels at NLO accuracy, the calculation of which represents the main technical result of the present work.
The impact of the approximate NNLO corrections on the stop pair production cross section has been examined. In order to obtain the best possible predictions and a better control on the neglected subleading terms, we have averaged the approximate NNLO results obtained in PIM and 1PI kinematics and matched them with exact fixed-order NLO calculations. As in the case of slepton-pair production studied in [17], we found that the total cross section depends strongly on the mass of the produced particles, while the dependence on the remaining SUSY parameters, such as the masses of other supersymmetric particles, is rather weak. We have found that including the approximate NNLO predictions for the pair production cross section reduces the perturbative uncertainty by more than a factor of two with respect to the NLO results. We stress that we explicitly accounted for the kinematic scheme uncertainty arising from the use of approximate NNLO formulas and we combined it with the scale uncertainty, as explained in Section 5. It was found that the approximate NNLO corrections have only a moderate impact on the central value of the NLO total cross section. We compared our result with analogous calculations at NLO+NLL accuracy [7,19] which were carried out by means of methods which are different from the ones employed in this work and within different kinematic schemes. We found very good agreement with the results of [7] and a substantial agreement within the quoted perturbative uncertainties with [19], although the cross section values we found tend to be slightly smaller than the ones quoted in [19]. We conclude by observing that, after the NNLO corrections are included, the main theoretical uncertainty on the total cross section comes from the PDFs. This is particular evident when large values of the stop mass are considered. Therefore, the results presented in this paper allow one to improve the precision of the stop-pair total cross section predictions by reducing their scale dependence to the extent that it becomes considerably smaller than PDF uncertainty. space is defined as (see Eq. (39) in [21]) where T (T ) indicates time ordering (anti-time ordering), while O s is the operator The four soft Wilson lines S are oriented along the directions of the incoming partons (n andn), and along the four-velocities of the outgoing top squarks (v 3 and v 4 ). The precise definition of the soft Wilson lines can be found in Eqs. (22) in [21]. The Laplace transformed soft functions can be obtained from the soft functions in position space through the relatioñ s (L, µ) = W x 0 = −2i e γ E µe L/2 , µ .
Here and below, we have omitted the dependence of the soft functions on the PIM or 1PI kinematic variables and on the heavy particle masses, as well as the subscripts qq or gg indicating the channel. The expansion ofs in powers of the strong coupling constant is At leading order the soft functions are the same both in PIM and 1PI kinematics: The bare soft function at NLO in position space can be written as The matrices w ij are related to the products of color generators and are the same for both kinematics. In the quark annihilation channel they are while for the gluon fusion channel one finds The functions I i ij are integrals over the soft-gluon phase-space. In PIM kinematics one finds I PIM where β t = 1 − 4m 2 t 1 /M 2 coincides with β in the z → 1 limit. Furthermore x s = (1 − β t )/(1 + β t ), y t = −t 1 /m 2 t 1 − 1, z u = −u 1 /m 2 t 1 − 1, and In 1PI kinematics one finds again that I 1PI 11 = I 1PI 22 = 0, while In the case of 1PI kinematics, the definition of β t should be changed to β t = 1 − 4m 2 t 1 /s .