NNLL momentum-space resummation for stop-pair production at the LHC

If supersymmetry near the TeV scale is realized in Nature, the pair production of scalar top squarks is expected to be observable at the Large Hadron Collider. Recently, effective field-theory methods were employed to obtain approximate predictions for the cross section for this process, which include soft-gluon emission effects up to next-to-next-to-leading order (NNLO) in perturbation theory. In this work we employ the same techniques to resum soft-gluon emission effects to all orders in perturbation theory and with next-to-next-to-logarithmic (NNLL) accuracy. We analyze the effects of NNLL resummation on the stop-pair production cross section by obtaining NLO + NNLL predictions in pair invariant mass and one-particle inclusive kinematics. We compare the results of these calculations to the approximate NNLO predictions for the cross sections.


Introduction
After the discovery of a Higgs boson in 2012, the search for supersymmetric partners of the Standard Model particles is one of the most important goals of the experimental program at the Large Hadron Collider (LHC). This is especially true in the light of the fact that the LHC is expected to run at a center-of-mass energy of 14 TeV when operations will resume in 2015. Consequently, the LHC will be able to further investigate the existence of supersymmetric particles with masses in the TeV range. The Minimal Supersymmetric Standard Model (MSSM) with R-parity conservation predicts the production of supersymmetric particles in pairs. At a hadron collider such as the LHC, one expects to observe primarily the production of supersymmetric particles carrying color charge, such as squarks and gluinos. In unified supersymmetric theories the third generation of squarks can have masses which are significantly lighter than the masses of the first two generations of squarks, as a consequence of large Yukawa and soft couplings entering the evolution of the mass parameters from the unification scale down to low energies. Consequently, the lightest of the two supersymmetric partners of the top quark could be the lightest squark in the spectrum and the first supersymmetric particle to be observed at the LHC.
This fact motivated the work of several groups, who in the last fifteen years obtained predictions for the top-squark pair production cross section with an accuracy beyond the JHEP03(2014)066 leading order in supersymmetric quantum chromodynamics (SUSY-QCD). The calculation of the next-to-leading order (NLO) corrections to the production of top-squark pairs was carried out in [1]. The NLO corrections enhance the production cross section if the renormalization and factorization scales are chosen close to the top-squark mass. The NLO corrections to the stop production cross section, as well as the corrections to the production cross section for several other supersymmetric particles, are implemented in the public codes Prospino and Prospino 2 [2]. Electroweak corrections to stop pair production have a quite sizable effect on the tails of the invariant-mass and transverse-momentum distributions, but they only have a moderate impact on the total cross section. These corrections were evaluated in [3,4]. Recently, phenomenological analyses which consider squark pair production and decay at NLO in QCD were carried out in [5,6].
Corrections from soft gluon emissions account for a large fraction of the NLO SUSY-QCD corrections. For this reason the resummation of these corrections to next-to-leading logarithmic (NLL) accuracy was studied in [7] by means of a standard technique based upon the resummation of threshold logarithms in Mellin space. Within the same approach, the resummation of the soft gluon corrections in the production of pairs of gluinos and squarks of the first two generations was carried out up to next-to-next-to-leading logarithmic (NNLL) accuracy in [8][9][10]. Furthermore, in [11][12][13] approximate next-tonext-to-leading order (NNLO) formulas for squark, stop, and gluino pair production were derived by means of resummation techniques. These formulas include threshold corrections and Coulomb corrections. Recently, the fully differential NLO cross section for the production of squark pairs, including squark decays, was evaluated and matched with parton showers [14].
Over the last few years, an alternative approach to resummation, which makes use of soft-collinear effective theory (SCET) methods, has been developed [15,16]. This approach, which allows one to work directly in momentum space, was applied to several processes of interest in collider physics, such as Drell-Yan production [17], Higgs production [18,19], direct photon production [20], top-quark pair production [21][22][23], and slepton pair production [24]. A similar method, combined with non-relativistic QCD techniques employed to resum Coulomb corrections, was independently developed in [25] in order to study top-quark pair production at NNLL accuracy, and it was applied to the study of squark-pair, gluino-pair and stop-pair production at NLL accuracy [26,27].
In particular, the studies of the top-pair differential distributions performed in [21,22] can be repeated in a straightforward (but laborious) way for the production of top squarks. The method adopted in those works relies on the factorization of the partonic cross section, which takes place in the soft limit. The partonic cross section can in fact be expressed as the convolution of two different factors. Schematically, these factors are a hard function, which includes the effects of virtual corrections, and a soft function, which describes the emission of soft gluons from the external particles involved in the process. Two different kinds of soft limits were considered. In [21], where the goal was the calculation of the invariant-mass distribution of the top-quark pair, the soft limit was defined as z = M 2 /s → 1, where M is the pair invariant mass and s is the square of the partonic center-of-mass energy. This framework is conventionally referred to as "pair invariant mass" (PIM) kinematics. In [22] JHEP03(2014)066 instead, the goal was the calculation of the top-quark transverse-momentum and rapidity distributions, therefore the soft limit was defined as s 4 → 0 with s 4 = (p 4 + k) 2 − m 2 t , where p 4 is the momentum of the unobserved anti-top quark, k is the momentum of the additional real radiation in the final state, and m t is the top-quark mass. The latter set up is know in the literature as "one particle inclusive" (1PI) kinematics. One can obtain predictions for the total cross section by integrating either one of the two kinds of distributions over the whole available phase space. The two calculations of the total cross section differ numerically because they neglect different sets of terms that are formally subleading in the soft limit. Total cross-section predictions which account for this source of uncertainty are obtained by averaging results from calculations carried out in the two kinematic schemes [28]. The calculations in PIM and 1PI kinematics share the same hard functions but require different soft functions. Furthermore, the soft functions in the two schemes do not depend on the spin of the particles involved in the process, so that the NLO soft functions employed in [21,22] are the same soft functions that one needs in order to study top-squark pair production.
In [29], the hard functions for the production of top-squark pairs was evaluated up to NLO in SUSY-QCD. By combining the hard functions with the soft functions of [21,22], it was possible to obtain approximate NNLO predictions for the pair invariant-mass and stop transverse-momentum distributions. At the moment, the most relevant observable in the study of top-squark production is the total cross section, which is employed in order to set lower bounds on the allowed values of the top-squark mass. For these reasons, in [29] we integrated the differential distributions in order to obtain approximate NNLO formulas for the total top-squark production cross section. A systematic comparison showed that the approximate NNLO results of [29] are compatible with the NLL results of [7] and [27].
In this work we study the resummation of the soft-gluon corrections to the top-squark pair production cross section by solving the renormalization-group equations (RGEs) satisfied by the soft and hard functions. The known expressions for the relevant anomalous dimensions [30,31] along with the NLO hard and soft functions calculated in [29] and [21,22], respectively, are sufficient in order to carry out the resummation up to NNLL accuracy. The paper is organized as follows: in section 2 we summarize our notation and conventions. In section 3 we review the resummation procedure, which uses the same scheme adopted in [21,22] for the study of top-quark pair production. In section 4 we discuss the matching of the NNLL resummation considered here to fixed-order NLO calculations; the matching is carried out in order to obtain NLO+NNLL predictions for the total cross section. The phenomenological impact of these predictions and their relations to other studies found in the literature are presented in section 5. Finally, we collect our conclusions in section 6.

Notation
The production of top-squark pairs is described by the scattering process

JHEP03(2014)066
We focus on the stop production at the LHC, so that N 1 and N 2 indicate the incoming protons, while X is an inclusive hadronic final state. In this work, we treat the top squarks as on-shell particles and neglect their decay; this approximation introduces an uncertainty of order Γt 1 /mt 1 , where mt 1 is the stop mass and Γt 1 represents its width. The two partonic subprocesses contributing to the stop pair production at lowest order in perturbation theory are The momenta of the incoming partons p i (i = 1, 2) are related to the hadronic momenta through the relation p i = x i P i . The relevant invariants for the hadronic scattering process are In order to describe the partonic scattering, we employ the Mandelstam invariants In Born approximation s + t 1 + u 1 = 0, and consequently M 2 = s and s 4 = 0.
Following the procedure employed in [29] and in the papers devoted to the calculation of differential distributions for top-quark pair production [21,22,32,33], we consider two different kinematic schemes, each of which has its own threshold limit. In PIM kinematics the threshold region is defined by the limit s → M 2 , while in 1PI kinematics the threshold region is approached by taking the limit s 4 → 0. The two different kinematics are suitable for the calculation of different differential distributions: PIM kinematics is used in order to calculate the pair invariant-mass distribution, while 1PI kinematics is employed in order to evaluate the stop transverse-momentum and rapidity distributions. In contrast to the production threshold region, which is defined by the limit β = 1 − 4m 2 t 1 /s → 0 and is often employed in the calculation of the total cross section in the soft limit, in the PIM and 1PI threshold regions top squarks are not necessarily produced nearly at rest. For instance, if we require to observe a stop pair with an invariant mass M , the squared partonic centerof-mass energy should be larger than M 2 , which can be much larger than the production threshold s ≥ 4m 2 t 1 . In both kinematic schemes, the partonic cross section in the threshold region is numerically dominated by the contribution of soft gluon emission.
A fact which is particularly relevant for resummation purposes is that in the soft limit the partonic cross section factors into products of hard and soft functions. Each of these two factors satisfies known RGEs. The anomalous dimensions entering these equations are know up to NNLO [30,31], while the matching coefficients are known up to NLO. This allows one to solve the RGE in Laplace space [15,16] and obtain resummed formulas which are valid up to NNLL accuracy.

PIM kinematics
In order to deal with PIM kinematics, it is useful to introduce the quantities the threshold region is defined by the limit z → 1. Because of the QCD factorization theorem [34], the double-differential cross section in M and θ (the stop scattering angle in the partonic rest frame) can be factorized as where µ f is the factorization scale, and the sum runs over the incoming partons. 1 As usual, parton luminosities ff ij are defined as the convolutions of the non-perturbative parton distribution functions (PDFs) for the incoming partons: The functions C ij in eq. (2.6) are the hard-scattering kernels, which are related to the partonic cross sections and can be calculated in perturbation theory. The hard-scattering kernels depend on the top-squark masses mt 1 and mt 2 (where we assume mt 1 < 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 thet 1 -t 2 mixing angle α. However, in order to avoid the use of an unnecessarily heavy notation, we drop these quantities from the list of arguments of the hard-scattering kernels.
At lowest order in α s , only the quark annihilation and gluon fusion channels contribute to the hard-scattering kernels, therefore ij ∈ {qq, gg}. In order to go beyond leading order, one needs to consider virtual and real emission corrections to the Born approximation, so that new production channels such as qg →t 1t * 1 q open up. However, it is a well-known fact that 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 partonicthreshold limit. Therefore, eq. (2.6) can be rewritten as where we omit terms of O(1 − z). In eq. (2.8) 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 JHEP03(2014)066 of eq. (2.8) differ in the fact that in the first term the quark (antiquark) comes from the hadron N 1 (N 2 ) in eq. (2.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 limit z → 1, the hard-scattering kernels C ij factor into a product of hard and soft functions [21]: 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 . 2 Throughout this paper, we work in the s-channel singlet-octet basis already employed in [29]. A factorization formula analogous to eq. (2.9) for the top-quark pair production was derived by employing SCET and heavy-quark effective theory in [21]. A completely analogous procedure can be followed in order to derive eq. (2.9), which is valid in the case of top-squark pair production. The hard functions, computed in [29], are obtained from 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. The soft functions are identical to the ones needed for the case of top-quark pair production, which were evaluated up to NLO in [21]. The hard functions were evaluated up to NLO in [29]. The RGEs satisfied by the hard and soft functions are identical to the ones satisfied by the corresponding quantities in the top quark production case and are discussed in detail in [21]. The anomalous dimensions regulating these RGEs are know up to NNLO. As discussed in section 3, by solving these RGEs it is possible to implement the resummation of soft gluon emission corrections up to NNLL accuracy.

1PI kinematics
1PI kinematics is used whenever one needs to consider kinetic properties of a single particle, rather than of the pair. One can then write the double-differential distribution in the topsquark transverse momentum and rapidity as (2.10) Obviously, only the quark-annihilation and gluon-fusion channels contribute to the hardscattering kernels C ij at the lowest order in α s . The hadronic Mandelstam variables T 1 and U 1 can be expressed in terms of the stop rapidity and transverse momentum as . Therefore, the variables s, s 4 , t 1 , u 1 , which are arguments of the 1PI scattering kernels, can be expressed in terms of p T , y, x 1 , x 2 . The lower integration JHEP03(2014)066 limits in eq. (2.10) are In order to obtain the total cross section, it is necessary to integrate the double-differential distribution in eq. (2.10) 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 limit s 4 → 0 factor into a product of hard and soft functions, in analogy to eq. (2.9): (2.14) 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 . We deal with this ambiguity following the methods described in section 4 of [22]. As in the case of PIM kinematics, the hard and soft functions are matrices in color space, arising from virtual and soft-emission corrections, respectively. The 1PI hard functions are identical to the ones encountered in the PIM kinematics. The 1PI soft functions, which differ from those derived in PIM kinematics, depend on plus distributions which are singular in the limit s 4 → 0. They were originally computed up to NLO in [22] for the top-quark pair production cross section. The RGEs satisfied by the hard and soft functions are identical to the ones discussed in [22], therefore all of the elements are in place to implement the resummation up to NNLL accuracy.

Resummation
Our main goal is to resum the leading singular terms in (1 − z) (PIM kinematics) or s 4 (1P1 kinematics) in the region of (partonic) phase-space where the stop production cross section is dominated by the threshold terms. This is accomplished by deriving and solving RGEs for the hard and soft functions. The RGEs for the hard functions do not depend on the virtual particles running in the loops or on the spin of the final state particles, and therefore they are precisely the same equations that have been discussed and solved in [21,22] up to the order appropriate for NNLL resummation. The RGE satisfied by the PIM soft functions and its solution can be found in section 5.1 of [21], while the solution of the RGE satisfied by the 1PI soft functions can be found in section 3.2 of [22]. Here we limit ourselves to collect the resummation formulas for the hard-scattering kernels appearing in eqs. (2.6) and (2.10).

JHEP03(2014)066
The resummed expression for the hard-scattering kernels in PIM kinematics is where we dropped the indices indicating the partonic channel. The channel-dependent hard matrices H are described in section 3.1 of [29], where they were evaluated up to NLO. The Laplace transform of the soft matrices,s, was defined in section 4.2 of [21]. The introduction of the Laplace transform of the soft matrices is motivated by the fact that, in Laplace space, soft functions are regular polynomials of their first argument, which satisfy ordinary first-order differential equations [15]. The PIM evolution matrices U and the exponential factor a γ φ are defined in section 5 of [21]. The parameter η arises from the solution of the RGE for the Laplace-transformed soft functions. The notation is such that one must first take the derivatives with respect to η appearing in the first argument ofs and then set η = 2a Γ (µ s , µ f ), as discussed in section 5 of [21]. For values µ s < µ f one finds that η < 0 and consequently one must use a subtraction at z = 1 and analytic continuation to express integrals in terms of plus distributions [35]. For example, for a smooth function g(z) that is not singular for z → 1 one can analytically continue the integrals from the region η > 0 to the region η > −1/2 by means of the relation If necessary, it is possible to analytically continue the integral on the left-hand side of this equation to the region η > −n/2 for an arbitrary positive integer n. This can be done by subtracting an increasing number of terms from the Taylor expansion of g(z) at z = 1.
Although the all-order hard-scattering coefficients C depend on the factorization scale µ f but do not depend on the soft and hard scales µ s and µ h , any practical implementation of the resummation formula eq. (3.1) will have a residual dependence on these two scales. This is due to the fact that the anomalous dimensions appearing in the evolution factors in eq. (3.1) are evaluated up to a given finite order in perturbation theory. The order at which this truncation takes places, together with the order at which the hard and soft functions are evaluated, defines the accuracy at which the resummation formula is implemented. The anomalous dimensions and the hard and soft functions are known at an order which is sufficient to carry out the resummation with NNLL accuracy. The choice of the numerical values for the hard and soft scales is discussed in section 4.
The resummation formula for the hard-scattering kernels in 1PI kinematics is (see section 3.2 in [22]) The evolution factors and the hard functions in eq. (3.3) are the same as in the PIM case (see [21]). The Laplace transform of the 1PI soft functions was evaluated up to NLO and can be found in section 3.1 of [22]. As for the PIM case, for values of the scale such that η < 0 one must use analytic continuation to interpret the formula in terms of plus distributions. Also in the case of 1PI kinematics, the resummation of the top-squark pair production cross section can be carried out at NNLL accuracy.

Matching and scale choices
Although the method employed allows us to obtain predictions for the pair invariant-mass distribution of the stop pair and for the transverse-momentum and rapidity distribution of a single top squark, we will limit ourselves to the calculation of the observable of phenomenological interest at the moment, i.e. the total stop-pair production cross section. The total cross section can be obtained by integrating the double-differential distributions in PIM and 1PI kinematics over the complete phase space, as explained in section 2.
Obviously, one wants to combine NNLL resummation with the most accurate fixedorder calculations of the total cross section available to date. Currently, the total stop-pair production cross section is known at NLO [1]. The NLO calculations can be matched to NNLL calculations of the total cross section as follows: where the subscript i ∈ {PIM, 1PI} indicates the kinematic scheme employed. Furthermore, the subscripts in eq. (4.1) indicate the scales (µ f , µ h , µ s ) on which each term depends. In eq. (4.1), σ NLO i is the exact result in fixed-order perturbation theory, while captures the leading singular terms in the threshold limit. If the various scales are set equal to each other, the resummed expressions for the cross section automatically reduce to fixedorder perturbative expansions. Consequently, the second term in the first line of eq. (4.1) includes the set of NLO terms which are not included in the resummed formulas, and it can be added to the first term, which includes the NNLL corrections, without introducing any double counting. The issue of the choice of numerical default value for the scales in the first term on the right-hand side of eq. (4.1) is addressed below. NLO predictions for the stop pair cross section can be conveniently obtained from the programs Prospino and

JHEP03(2014)066
Prospino2 [2]. The matching procedure of eq. (4.1) can be carried out separately for each of the two kinematic schemes considered. Since the total cross section can be obtained starting from either of the two kinematics, but each kinematics neglects different sets of subleading corrections, we follow the procedure already adopted in [28,29] and average the two results. Schematically, our resummed prediction for the total cross section is then obtained as Similarly, in evaluating the perturbative error associated with our result, we want to reflect also the uncertainty associated to the choice of the kinematic scheme. In order to achieve this goal, we start by varying separately each scale where µ 0,i denotes the default choice for the scale µ i , which is discussed in the next two sections. We then evaluate the quantities where we neglected the subscript NLO + NNLL for each of the cross sections appearing on the left-hand side of eqs. (4.4). In complete analogy, we also evaluate the quantities ∆σ ± h and ∆σ ± s by varying the hard or soft scales, while keeping the other two scales equal to their default values. Finally, the perturbative uncertainty on the cross section is obtained by combining the quantities ∆σ ± i in quadrature, i.e.
At this stage we turn our attention to the choice of the default values for the soft, hard and factorization scales.

Choice of the hard and factorization scales
The hard scale µ h should be set to the characteristic scale of the underlying partonic subprocesses shown in eq. (2.2). An obvious possibility would be the invariant mass M of the stop pair, which is the lower bound on the partonic center-of-mass energy √ s. However, the observable M is only defined in PIM kinematics, whereas the pair invariant mass is not observed in 1PI kinematics. We will therefore use the other obvious possibility, the production threshold µ 0,h = 2mt 1 , as the default value for the hard scale in both kinematic JHEP03(2014)066 schemes. For the factorization scale, we follow the standard choice made in fixed-order perturbation theory calculations, namely we set µ 0,f = mt 1 . As is common practice, we will vary the scales µ h and µ f independently by factors of 2 about the default values.

Choice of the soft scale
Contrary to the hard matching scale, the soft matching scale is not associated with a parameter entering the partonic cross sections. Rather, it is generated dynamically when the partonic cross sections are convoluted with the steeply falling PDFs [16]. Our procedure for fixing the value of the soft scale is similar to the one employed in the case of top-quark pair production in [21,22]. In the case of the top-squark pair production considered here, the problem is slightly more complicated because the stop mass is not known, and it becomes a parameter in the determination of µ s . In general, one expects to find that the soft function has a well-behaved perturbative expansion when µ s is set equal to a scale characteristic of the energy of the real soft radiation, which is expected to be smaller than the hard scales mt 1 and √ s. In order to find this scale for a given kinematic scheme and fixed center-of-mass energy and mt 1 , we look for the minimum of the α s corrections to the total cross section arising from the soft function as a function of µ s . In order to isolate these corrections, we select the part of the NNLL resummed formula for the hardscattering kernels which arises froms (1) (i. e. the NLO contribution to the soft function), evaluate the contribution of these terms to the total cross section, and divide what we find by the NLL cross section. We furthermore set µ s = µ f = µ h , which is equivalent to considering the fixed-order corrections at NLO accuracy. When plotting these corrections as a function of µ s /mt 1 for fixed s and mt 1 , one finds that they show a minimum. We further plot the location of the minimum as a function of mt 1 . The curve which emerges is that of a smooth, monotonically decreasing function, which for fixed kinematics and collider energy can be well approximated by a quadratic polynomial. We employ such fits in order to determine the default value of the soft scale for fixed mt 1  A similar curve is found in the case of 1PI kinematics. The resulting functions are shown in figure 1. In order to account for the uncertainty introduced by the scale choice, in phenomenological predictions we allow the chosen soft scale to vary in the range [µ 0,s /2, 2µ 0,s ], as explained above.

Phenomenology
In this section we analyze the numerical predictions for the stop pair production cross section at NLO+NNLL accuracy. In particular, i) we compare the results obtained in PIM and 1PI kinematics and their average, ii) we investigate the dependence of the predictions on the variation of the hard, soft and factorization scales, iii) we provide numerical tables for different values of the stop mass and for different choices of the PDF sets, and iv) we compare the predictions with NLO+NNLL accuracy to the approximate NNLO cross section studied in [29]. In order to keep the presentation concise, we consider two values for the LHC center-of-mass energy: √ S = 8 TeV, which is the energy at which the machine was running before the shutdown in 2013-2014, and √ S = 14 TeV, which is the targeted energy when operations resume in 2015. Furthermore, as in [29] we fix the SUSY parameters other than the light stop mass to the value characterizing the benchmark point 40.2.5 in [36]. As it was shown in [29], the total cross section shows little sensitivity to the SUSY parameters other than mt 1 . Table 1 collects the values of the input parameters entering the hard functions employed throughout this section. The benchmark point 40.2.5 uses a value of 1087.15 GeV for the top-squark mass. We employ this value in the tables below. 3 However, in the same tables we also consider mt 1 = 500 GeV, which is representative of the current experimental lower bounds on this quantity. In addition, we plot mass scans for the total cross section in the range mt 1 ∈ [500, 2000] GeV. In the following, unless we explicitly write that we do otherwise, it is understood that we employ NNLO PDFs in NLO+NNLL calculations and approximate NNLO calculations, while we employ NLO PDFs in NLO and NLL calculations. In each plot or table, we explicitly indicate the use of either CT10 [37,38] or MSTW2008 [39] PDFs.

Comparison between 1PI and PIM kinematics
Calculations which rely on the use of PIM and 1PI kinematics neglect different sets of power-suppressed terms and therefore lead to numerically different predictions. In order to account for the scheme uncertainty, we combine the NLO+NNLL predictions in the two kinematic schemes as explained in section 4. The differences of the predictions obtained using PIM and 1PI kinematics can be inferred from figure 2, where the two dark solid lines are obtained by considering, for each value of mt 1 , the quantities where σ without superscript indicates the average between the 1PI and PIM predictions, obtained according to eq. (4.3). To obtain these lines all scales (soft, hard, and factorization) are set at their default values discussed in section 4. In both panels, the 1PI prediction for the resummed cross section is slightly larger than the PIM prediction in the entire range of values for mt 1 considered in the figure. However, in both cases the spread between the 1PI and PIM predictions is significantly smaller than the perturbative uncertainty of the combined result, represented by the light brown band and determined as discussed in section 4. The slight dent in the bands at mt 1 ≈ 1660 GeV, which is particularly evident in the right panel of the figure, coincides with the gluino-top-quark production threshold.

Scale dependence of the resummed cross section
An anticipated effect of the resummation at NNLL order is that phenomenological predictions should be less sensitive to the choice of the soft, hard, and factorization scales when compared to calculations at NLL accuracy. We study this aspect in figure 3. In all panels the top-squark mass is set equal to 1087 GeV. The plots in the left column refer to a LHC center-of-mass energy of 8 TeV, while the ones on the right column refer to 14 TeV.
The two panels in the first row shows the effect of varying the factorization scale about its standard value µ f = mt   In these two plots, the hard and factorization scales are set to their default values. Finally, all panels in figure 3 indicate that the NNLL corrections increase the cross section with respect to NLL calculations. We remind the reader that in order to asses the total scale uncertainty of the 1PI and PIM predictions both at NLL and at NNLL, we first vary each scale in the range [µ 0,i /2, 2µ 0,i ] (i = f, s, h) and then add the three uncertainties obtained in this way in quadrature. In view of the behavior shown in the plots, one can expect a slightly smaller scale uncertainty at NNLL than at NLL. At this stage we proceed to discuss the effect of the NNLL corrections on the total stop production cross section.

Total cross section
We now proceed to study the effect of the resummation at NNLL accuracy on the total stop pair production cross section. We study this aspect by comparing the NNLL predictions for the cross section (matched to the fixed-order NLO cross section) with NLO, NLL, and approximate NNLO predictions for the same observable. In this section, LO   We remind the reader that PDFs at different orders (and, consequently, the cross-section predictions) employ different values of α s . We start by discussing the stop production cross section at two different values of the top-squark mass: mt 1 = 500 GeV and mt 1 = 1087 GeV. Tables 2 and 3 show the predictions for the LHC operating at a center-of-mass energy of 8 TeV, obtained using the PDF sets   MSTW2008 and CT10, respectively. Tables 4 and 5 show the corresponding results for the case in which √ S = 14 TeV. The K-factors are defined as Inspecting the tables, we observe that the NLL calculations predict a smaller cross section than the NLO calculations. This effect is particularly pronounced at mt 1 = 500 GeV. This means that NLO contributions not included in NLL soft gluon emission corrections are JHEP03(2014)066 numerically sizable, in particular for smaller values of the stop mass, where hard gluon emission is less suppressed by phase-space constraints. A similar behavior was already encountered in the study of the top-quark pair production cross section (see for example table 4 in [21]).
This situation should be compared to the relation between approximate NNLO and NNLL +NLO predictions. While the total stop production cross section at NLO+NNLL accuracy is slightly smaller that the approximate NNLO cross section for all cases considered in the tables, the two predictions are well within the respective perturbative uncertainties, which are indicated by the first error next to the central values reported in the tables. The second error in the tables accounts for the PDF and α s uncertainty. Both approximate NNLO and NLO+NNLL cross sections agree within perturbative uncertainties with the NLO calculations. The relative size of the approximate NNLO and NLO+NNLL corrections in the mass range mt 1 ∈ [500, 2000] GeV is shown in figure 4. As already observed in the tables, the approximate NNLO cross section is slightly larger than the NLO+NNLL one except in the case of very large mt 1 masses (mt These considerations serve as a posteriori self-consistency check of our calculational framework and indicate that the approximate NNLO and the NLO+NNLL predictions, which are in good agreement with each other, are robust. Of course, a full calculation of the NNLO corrections to the stop-pair production process would be the only way of assessing with certainty to which extent approximate NNLO calculations reproduce the exact NNLO results. Furthermore, NNLO calculations in fixed-order perturbation theory could be easily matched to the NNLL resummed cross section discussed in this work. Unfortunately to date the large number of mass scales involved makes a full evaluation of the NNLO corrections an extremely challenging task.

Comparison with other results in the literature
We conclude our phenomenological analysis by comparing our NLO+NNLL predictions for the total cross section with the results obtained in [7] and [27], which have NLO+NLL JHEP03(2014)066  Table 6.
Comparison between the NLO+NLL cross section of [7], the approximate NNLO calculation of [29], and the NLO+NNLL result of the present work. The table refers to the LHC with √ S = 14 TeV and to mt 1 = 400 GeV, the remaining input parameters are set at the values characterizing the SPS1a' benchmark point in [40]. The PDFs employed are the CTEQ6.6 set. We report only the perturbative uncertainty. accuracy and are obtained with calculational methods different from the ones employed here. We focus our attention on values of the stop mass which are close to or higher than the current lower bounds on this parameter in the MSSM.
In table 6 we show the results obtained for the input parameters employed in [7], which coincide with the SPS1a' benchmark point in [40]. In the table we consider a collider energy of 14 TeV and set mt 1 = 400 GeV. The PDF set employed is CTEQ6.6. We checked that, as expected, the NLO results in [7] coincide with the ones obtained with the Prospino version JHEP03(2014)066 we employ. Our central value for the NLO+NNLL cross section is in very good agreement with the NLO+NLL value obtained in [7] and with the approximate NNLO prediction obtained in [29]. The perturbative uncertainty of the NLO+NNLL result is essentially identical to the one affecting the approximate NNLO calculation. Both are smaller than the NLO+NLL scale uncertainty.
Reference [27] presents results obtained by resumming simultaneously production threshold logarithms and Coulomb singularities with NLL accuracy. Bound-state effects are also included in that calculation. Results for the top-squark pair production at the CMSSM benchmark point 40 Table 7. Comparison between the NLO+NLL results of [27], the approximate NNLO calculation of [29], and the NLO+NNLL results of the present work. The numbers refer to the benchmark point 40.2.4 [36]. In particular, we set m t = 172.5 GeV, mg = 1386 GeV, mq = mt 2 = 1358 GeV, and cos α = 0.39 as in [27]. The numbers refer to the LHC operating at √ S = 7 TeV (upper portion) and 8 TeV (lower portion). The factorization scale is set equal to mt 1 . For the approximate NNLO results and the NLO+NNLL results we used MSTW2008 NLO PDFs. The errors indicate only the perturbative uncertainty.
The central values at NLO+NNLL accuracy are marginally smaller than in approximate NNLO calculations. The perturbative uncertainty is slightly larger than the one found at approximate NNLO, but smaller than the one quoted in [27]. The lower portion of table 7 shows that the same observations apply to the case of the LHC at √ S = 8 TeV, for which the authors of [27] provide predictions in an ancillary file included in the arXiv submission of their work.
Finally, we briefly comment on a few papers which were recently published. Ref. [41] analyzes the impact of finite-width effects on threshold corrections to squark and gluino production, finding them to be negligible for a moderate decay width, Γ/m ≤ 5%, which corresponds to the case of interest for present searches. This result confirms the validity of the analysis presented here, which neglects these effects. Refs. [9,42] present the first results in threshold resummation (in the β → 0 limit) for squark and gluino production at NNLL accuracy (the latter including Coulomb gluon effects as well). Since these papers focus on squark and gluino production and do not consider stop pair production, a direct comparison is not possible. It will be interesting to compare the different approaches when a comprehensive phenomenological analysis for stop pair production will be available.

Conclusions
In this paper, we have completed the analysis of the soft-emission corrections to the production of top-squark pairs started in [29]. In particular, we have considered the resummation of partonic threshold logarithms at NNLL order. Our method relies on the factorization JHEP03(2014)066 of the partonic cross section in a trace of the product of two matrices, the hard and soft functions, in color space. This factorization is valid in the soft limit. The hard function accounts for virtual corrections, while the soft function accounts for the emission of soft gluons. In [29], it was shown that the use of the threshold limit of the partonic cross section allows one to obtain reliable predictions for hadronic observables in stop pair production, at least for the range of values of mt 1 considered in that work and in the present one. This happens because of the mechanism of dynamical threshold enhancement [17], which essentially amounts to the fact that PDFs enhance the relative weight of the soft-emission region in the partonic phase-space integrals appearing in the calculation of hadron-initiated production process. Furthermore, in [29] we presented the calculation of the hard and soft functions up to NLO and derived the anomalous dimensions of the hard and soft functions required in order to obtain approximate NNLO formulas for stop-pair production observables. These approximate formulas include all of the plus distributions appearing in the partonic cross section at NNLO, which capture the leading singular terms in the soft limit.
In the present work, we have solved the RGEs satisfied by the hard and soft functions in order to carry out the resummation of threshold logarithms directly in momentum space, with NNLL accuracy. The relevant anomalous dimensions are identical to the ones employed in the study of top-quark pair production considered in [21,22]. We have carried out the analysis in two different kinematic schemes, PIM and 1PI, which in principle allow us to obtain different differential distributions, such as the stop-pair invariant-mass spectrum or the top-squark transverse-moment and rapidity distributions.
However, top squarks have not been discovered yet. Consequently, the most interesting observable in top-squark pair production is the total cross section, which must be evaluated as a function of the mass of the hypothetical top squark. Our technique allows us to obtain the total cross section by carrying out the resummation in either PIM or 1PI kinematics and subsequently integrating the double-differential distribution over the available phasespace. Furthermore, the difference between the predictions obtained in the two kinematic schemes provides a handle for how to estimate subleading corrections neglected in the soft limit. In fact, a different set of formally subleading corrections is neglected in the two different schemes. This scheme uncertainty is combined with the usual scale uncertainties in order to estimate the total perturbative uncertainty affecting our predictions.
The phenomenological predictions for the total cross section as a function of the topsquark mass have been obtained by matching the resummation formulas at NNLL order with the complete NLO calculation obtained from the code Prospino [2]. The NLO+NNLL calculations lead to values of the total cross section which are very close to the approximate NNLO calculations of the same observable, first presented in [29]. The perturbative uncertainty affecting the NLO+NNLL calculations is essentially the same we found in approximate NNLO calculations; moreover, it is smaller than both the perturbative uncertainty affecting NLL calculations and the residual PDF and α s uncertainty. We consider the good agreement between approximate NNLO and NLO+NNLL calculations an indication of the fact that our calculational framework is self consistent and robust, a priori equivalent to other schemes employed to carry out the resummation of soft gluon emission in this process. We note, however, that the resummation of higher-order Coulomb corrections, studied for example in [27] at NLL accuracy, is not considered in the present work.

JHEP03(2014)066
Finally, we emphasize that the procedure described here and in [29] can be adapted and applied to the study of other production processes involving colored supersymmetric particles, such as gluino pairs, sbottom pairs, and pairs of squarks of the first and second generation.