Top-squark pair production at the LHC: a complete analysis at next-to-leading order

We present a complete next-to-leading order study of top-squark pair production at the LHC, including QCD and EW corrections. The calculation is performed within the Minimal Supersymmetric Standard Model and numerical results are presented for parameter regions compatible with the observed Higgs boson. We employ the most recent parton distribution functions including QED corrections and we find NLO EW corrections to the inclusive stop-pair production cross section up to $25 - 30\%$ compared to the leading-order prediction. Besides corrections to inclusive cross sections also important kinematic distributions are investigated.


Introduction
The discovery of a signal in the Higgs boson searches at the LHC [1,2] has triggered extensive studies of the properties of this particle, in particular of its mass, couplings and spin [3,4]. The experimental results are in agreement with Standard Model (SM) predictions. Nevertheless, experimental evidence for beyond the Standard Model (BSM) physics may still be found in the future run, in particular for the minimal supersymmetric extension of the Standard Model (MSSM), where the value of the Higgs boson mass can naturally be explained.
Experimentally the production of third generation squarks is different from those of the other generations. Final-state bottom and top quarks from the decays of the squarks can be distinguished from their light-flavor counterparts. This leads to characteristic signatures which are extensively used in the experimental searches for SUSY at the LHC [62][63][64][65][66][67]. These searches assume a dominant decay channel and put limits on the mass of the light stop,t 1 , and of the decay products. For instance the analysis in Ref. [67] puts a lower bound on thet 1 one, mt 1 ≥ 750 GeV, assuming 100% branching ratio into a massless neutralinoχ 0 1 . 1 In the experimental analyses signal events are usually produced at LO and the total cross section is normalized to NLO QCD + NLL accuracy using the codes Prospino [68] and NLL-fast [33], along the lines of Ref. [69]. The theoretical uncertainty on the generated signal is obtained by varying the parton distribution function (PDF) sets and the renormalization/factorization scale. For a hadronic center-of-mass energy √ S = 14 TeV, the obtained uncertainty is depicted in Fig. 1. It is below 15% (20%) for masses below 1000 (1400) GeV. It is worth to mention that this procedure does not provide an accurate estimation of the impact of the higher order contributions. Indeed, the latter may depend strongly on the kinematics of the final state changing the shapes of differential distributions. Higher-order contributions to the decays of the produced particles can have substantial effects as well [57][58][59][60].
In this paper we analyze the numerical impact of the NLO EW corrections to the total cross section of the hadronic pair-production of the lightest top-squarkt 1 , by presenting the first phenomenological study combining NLO QCD and NLO EW corrections. A similar study of the NLO EW corrections has already been published in [19]. There, results were presented for a handful of Snowmass points [70] and for benchmark slopes passing along the SPS1a ′ scenario [71]. NLO EW corrections to inclusive cross sections were found to be small, below 5%, while differential distributions were found to be affected by up to 20% in their tails. In view of the developments over the recent years, updated and improved analyses are considered appropriate. Firstly, the previous study focuses purely on the EW corrections without including the (dominant) NLO QCD corrections. Only a consistent treatment of NLO QCD and NLO EW including common PDFs allows for a quantitative comparison. Secondly, in the MSSM the mass of the lightest (SMlike) Higgs boson is not a free parameter, but depends crucially on the top-squark masses and the mixing parameter X t = A t − µ/ tan β. Consequently, the experimental signal for a Higgs boson with a mass around 125 GeV yields a strong constraint for the scalar-top sector that has to be taken into account in any phenomenological studies. In general, the Higgs-mass constraint requires heavy stops and/or large mixing in the stop sector. Finally, a consistent computation at NLO EW accuracy requires a PDF set which includes both the photon structure function and the QED contributions to the evolution equations. When the study in Ref. [19] was performed, the MRST2004QED set [72] was the only one fulfilling this requirement. This set is affected by large uncertainties that have been significantly reduced in the newly available NNPDF2.3QED PDF set [73]. Here, latest data from deep inelastic scattering and from the LHC are included. It is worth to investigate how the new PDF alters the EW corrections to the production of the lightest top-squark.
Our study mainly focuses on the NLO corrections to the inclusive stop-anti-stop production cross section, including corrections of electroweak origin. This analysis can be directly used in experimental analyses to estimate the theoretical uncertainty related to the missing EW corrections. Besides corrections to inclusive cross sections we also consider the impact of the electroweak corrections on several kinematic distributions, which may be significant in specific SUSY scenarios [19].
This paper is structured as follows. Section 2 reviews the calculation of the total cross section of stop-anti-stop pair production, including both NLO QCD and NLO EW contributions. Section 3 is dedicated to the presentation of the numerical results, followed by our conclusions in Section 4.

Computational details
The leading-order contribution to the hadronic cross section for the production of stopanti-stop pairs is of the O α 2 s , described as follows, indicating the perturbative order O α a s α b of the partonic processes X contributing to the total hadronic cross section as σ a,b X . The sum in Eq. (2.1) runs over the four lightest quark flavors, q = u, d, c, s. We neglect the bottom-initiated partonic processes, since they are suppressed by the bottom-quark parton distribution function and turn out to be small [28,32]. 2 The partonic cross sections for the gg and qq channel can be found in Eq. (3) and Eq. (4) of [15], respectively. They depend only on the mass of the produced top-squarks.  They have been computed in Ref. [15] and are implemented in the public code Prospino [68].
The NLO QCD corrections depend on the masses of the top squarks and of the gluino, and on the stop mixing angle.
The EW contributions arise at O α 2 , O (α s α) and O α 2 s α ; they can be separated into four different channels according to the partonic initial states, In principle these electroweak terms depend on the full set of MSSM parameters. We compute them by using FeynArts [75,76] and FormCalc [76,77], together with LoopTools [76] for the numerical evaluation of the one-loop integrals. Ultraviolet divergences are cancelled by renormalization at the electroweak one-loop level, along the lines of Ref. [28]. Infrared and collinear singularities are handled by using mass regularization and are computed by using the double cut-off phase-space-slicing method [78][79][80], as described in Ref. [24]. The initial-state collinear singularities of gluonic (photonic) origin are factorized and absorbed in the parton distribution functions by using the MS (DIS) scheme. Our computation has been numerically checked against the results presented in [19]. The EW contributions have been implemented in the code SusyHell, a (to be public) Monte Carlo integrator for the production of colored SUSY particles at the LHC.
The inclusive cross section for stop-anti-stop production, complete at NLO, is obtained by summing the LO cross section, Eq. (2.1), and the NLO contributions, Eqs. (2.2) and (2.3). For a discussion of the QCD and EW effects separately, it is convenient to define the individual and summed cross-section parts as follows,  Theoretical uncertainty (%) Figure 1: Theoretical uncertainty on the total cross section fort 1t * 1 production at the LHC. It has been computed by using the code NLL-fast [33] and the procedure described in Ref. [69].

Numerical analysis
The numerical values of the Standard Model input parameters are chosen according to The results presented in this section are computed for a hadronic center-of-mass energy of √ S = 14 TeV. For the numerical evaluation of the hadroinc cross sections, we use the NNPDF2.3QED PDF set [73]. The latter has been implemented in Prospino and SusyHell through the LHAPDF interface [81], which also automatically accounts for the consistent evolution of the strong coupling constant. The factorization scale µ F and the renormalization scale µ R are set to a common value, equal to the mass of the lightest top-squark, In the following we consider various phenomenological MSSM scenarios characterized by ten TeV-scale free parameters, 2) from which the individual soft-breaking parameters are obtained in the following simplified way: thereby, for M 1 gaugino-mass unification at the GUT scale is assumed. From these soft parameters we calculate the physical spectrum using tree-level relations. The only exception is the physical mass of the Higgs bosons, computed by using the code FeynHiggs 2.10 [82][83][84][85][86].

Total cross section
For our numerical evaluation we consider several benchmark scenarios as defined below. Within these scenarios we perform various one-dimensional scans over the parameters listed in (3.2). In these scans we focus on regions allowed by limits from the Higgs sector, i.e. we require a Higgs-state close to the observed one, with couplings compatible with the observed rates in the Higgs search channels. The compatibility between the considered models and the experimental results has been checked by using the codes HiggsBounds [87][88][89][90] and HiggsSignals [91]. They compute a χ 2 measure from the predictions of the model and the measured Higgs rates and masses. From this measure a p-value is estimated, testing the consistency between the model and the data. As a practical rejection criterion, we discard a model point if the corresponding p-value is below 0.0027, i.e. we require consistency at the three-sigma level.
The results of the scans are shown in Figs Light-Higgs scenario This scenario is inspired by the best-fit point of Ref. [92]. In that study a fit of phenomenological MSSM scenarios compatible with electroweak precision observables and with experimental searches at the LHC and the Tevatron was performed. The values of the parameters characterizing this scenario are shown in Table 1(a). The left (right) panels of Fig. 2 show the results of the scan over X t (Mq 3 ), and Fig. 3 displays the results of two scans over µ. Additional scans over M 2 and mg are very insensitive to the values of these parameters and we do not include them here.
The dependence of the total cross section on X t is mild and basically related to the dependence of mt 1 on X t . On the other hand, the total cross section depends strongly on Mq 3 and varies over four orders of magnitude in the considered parameter range, cfr. Fig. 2(b). For Mq 3 ≥ 1400 GeV (mt 1 ≥ 1300 GeV) the NLO cross section drops below 1 fb. Correspondingly high luminosities are required to exclude this scenario. In the high Mq 3 region the relative impact of the EW corrections become more important, up to 20 − 30% of the LO cross section. The relative impact of the EW corrections is enhanced by the behavior of the QCD corrections whose relative yield decreases from 50 to 40% above thet 1 →gt threshold, Mq 3 ≃ 1300 GeV. As can be inferred from Fig. 2(d), the EW corrections from the qq and qg channels are negligible, while the ones from the gg-channel  Table 1: Benchmark scenarios within our ten-parameter phenomenological MSSM.
are of the order of 5% of the LO cross section, irrespective to the value of Mq 3 . In this scan the gγ-channel dominates the EW corrections for Mq 3 ≥ 1000 GeV. Its inclusion in experimental studies would significantly reduce the theoretical uncertainties, i.e. the systematic uncertainty originated from neglecting the EW corrections. The LO and NLO QCD cross section is independent of µ, see Fig. 3(a) and Fig. 3(b). However, the NLO EW and thus also the full NLO cross section do depend on µ. The NLO EW corrections double their impact on the LO cross section when µ varies from 500 to 2500 GeV. In the left part of Fig. 3, where we take Mq 3 = 1000 GeV, their contributions are below 15%, i.e. four times smaller than those of the QCD corrections. The right scan, where we assume Mq 3 = 1250 GeV, exhibits larger EW corrections: they amount to more than 15% of the LO cross section for µ ≥ 1500 GeV. The relative impact of the various channels of the EW corrections for Mq 3 = 1000 GeV and Mq 3 = 1260 GeV is shown in Fig. 3(c) and Fig. 3(d) respectively. Again, the qq and gq-channel are negligible, while the contribution of the gγ-channel is constant and larger for Mq 3 = 1260 GeV. In the Mq 3 = 1000 (1250) GeV scan this channel contributes about 7% (11%) of the LO cross section. Instead, the contributions from the gg-channel strongly depend on µ and increase as this parameter increases. For Mq 3 = 1000 GeV it dominates the EW corrections for µ ≥ 2700 GeV. It is worth to note that in both cases the inclusion of the tree-level, model-independent gγ-channel would significantly reduce the EW theoretical uncertainty. In these scans we observe at 1 →χ 0 3 t threshold at µ ≃ 650 (950) GeV for Mq 3 = 1000 (1250) GeV.
Light-stop scenario This scenario has been defined in Ref. [93] in correspondence to the values of the parameters summarized in Table 1(b). The value of Mq 3 is significantly below the TeV-scale, while X t is chosen to maximize the mass of the lightest CP-even Higgs boson h, leading to a phenomenologically-viable scenario with a light stop.
The dependence of the total cross section on Mq 3 is shown in Fig. 4(a). Both the LO and the NLO cross section depend strongly on this parameter: their values span three orders of magnitude when Mq 3 varies from 400 to 1000 GeV. The relative impact of the QCD corrections is large, while its dependence on Mq 3 is mild: over the entire Mq 3 interval they are 50 − 55% of the LO cross section. Outside thet 1 →χ 0 4 t threshold region, the EW corrections are positive and small, i.e. they are 8% of the LO cross section at most. The relative impact of the EW corrections of the different channels as a function of Mq 3 is shown in Fig. 4(c). Outside of the threshold region the bulk of the EW contributions originate from the gγ-channel, owing to mutual cancellations between the gg-and the qq-channels. Fig. 4(b) shows the dependence of the total cross section on µ. The LO and NLO total cross sections are independent of the value of µ in a vast portion of the scanned region. The total NLO cross section decreases in the high µ region, i.e. µ ≥ 1200 GeV, as a consequence of the decrease of the EW corrections. The variation is below 10% and it takes place when approaching thet 1 →b 1 W threshold located in the phenomenologically excluded region, Mq 3 ≥ 1780 GeV.
As can be inferred from Fig. 4(d), the gγ-channel dominates in the µ = 300 − 600 GeV range. For higher values of this parameter, µ = 700 − 1200 GeV, the gg-and the gγchannel dominate and result in an overall negative NLO EW correction. The inclusion of the gγ-channel only would thus overestimate the impact of the EW corrections.
Light-stau scenario This scenario exhibits a sizable mixing in theτ sector and a rather low value of Ml 3 , allowing for a lightτ state and a possible enhancement of the h → γγ decay rate compared to its SM prediction [93]. The value of the parameters characterizing this scenario are listed in Table 1(c).
In the left (right) plots of Fig. 5 we investigate the dependence on Mq 3 (µ) within this scenario. The NLO QCD corrections to the total cross section do only marginally depend on Mq 3 and negligibly on µ. They amount to 50% of the LO cross section, cfr. Fig. 5(a) and 5(b).
The dependence of the EW contributions on Mq 3 is less trivial. Below thet 1 →χ 0 3 t threshold at Mq 3 ≃ 800 GeV the EW corrections are flat and small, below 5% of the LO cross section. Above the threshold the corrections grow as Mq 3 grows reaching 10% (15%) of the LO cross section for Mq 3 ≃ 1150 (1200) GeV, i.e. for mt 1 ≃ 1050 (1200) GeV. In the high Mq 3 region the cross section prediction is thus substantially modified by the inclusion of the EW corrections. As shown in Fig. 5(c), outside the threshold region the EW contributions are decently approximated by the gγ-channel only; indeed, the other channels contribute only up to 2% of the LO cross section.
Here, the EW corrections are only mildly affected by the variation of µ, cfr. Fig. 5(b). Far from the threshold region, µ ≃ 700 GeV, they are flat and of the order of 7− 10% of the LO cross section. The leading contribution to the EW corrections is the µ−independent contribution of the gγ-channel, which below (above) the threshold is slightly (considerably) enhanced by the other channels.
Tau-phobic scenario This scenario belongs to the two-parameter family of scenarios also introduced in [ The relative size of the QCD corrections, about 50% of the LO the cross section, is not affected by the variation of either Mq 3 or µ, see Figs. 6(a) and 6(b). The EW contributions, on the other hand, do depend on Mq 3 ; their relative size is tripled when Mq 3 varies from 1000 to 1800 GeV. Hence, the EW corrections are important and not negligible; they are more than 10% of the LO cross section over almost the entire region of Mq 3 considered, i.e. Mq 3 ≥ 1100 GeV or, equivalently, mt 1 ≥ 900 GeV. Moreover they can be as large as 25% of the LO cross section in the large Mq 3 region, Mq 3 ≥ 1700 GeV.
The dependence of the EW corrections on µ is milder but nevertheless significant. The EW contributions relative to the LO cross section increase for increasing µ and double their value as µ varies from 100 to 3000 GeV reaching 10% at µ ≃ 2500 GeV, cfr. Fig. 6(b).
In the considered parameter region mutual cancellations between the gg and qq channels render the gγ-channel the dominant one. Its inclusion would considerably lower the remaining EW theoretical uncertainty down to 5% of the LO cross section. The detailed behavior, however, strongly depends on µ, as can be inferred from Fig. 6(d). In particular approximating the EW contributions by just the gγ-channel becomes less and less accurate for increasing µ: the gγ-channel yields less than half of the EW corrections for µ ≥ 2000 GeV. m mod± h scenarios These scenarios are defined in Ref. [93] as modifications of the m max h scenario [94]. They lead to a smaller mass of the lightest CP-even Higgs boson by reducing the ratio The input parameters of the m mod + h and of the m mod − h scenario are collected in Table 1(e) and Table 1(f), respectively. They differ in the value of r t and in the sign of X t .
The dependence of the NLO corrections on Mq 3 in both scenarios is similar to the scenarios investigated before, see Fig. 7(a) and Fig. 8(a). The QCD corrections stay at 50% of the LO cross section over the entire Mq 3 interval. The relative impact of the EW corrections increases with Mq 3 . They are larger than 10% of the LO cross section for Mq 3 ≥ 1100 GeV, that is for mt 1 larger than 950 GeV. As shown in Fig. 7(c) and Fig. 8(c) mutual cancellations between the qq and gg-channel effects make the EW corrections effectively equal to the contribution of the gγ-channel. Fig. 7(b) shows the µ-dependence of the NLO contributions in the m mod + h scenario. Outside thet 1 →χ 0 4 t and thet 1 →b 1 W threshold region, located at µ ≃ 710 GeV and µ ≃ 2620 GeV respectively, the EW corrections relative to the LO cross section are almost flat and vary from 8% to 12% as µ varies in the considered interval. As can be inferred from Fig. 7(d), the largest EW contribution originates from the gγ-channel. Outside the threshold regions the combined effect of the other channels are 2% of the LO cross section at most. In the m mod − h scenario the relative yield of the EW corrections doubles its value from 7% to 15% as µ varies from 100 to 3000 GeV. In the low µ region, i.e. below thet 1 →χ 0 4 t threshold at µ ≃ 700 GeV, the gγ-channel is a good approximation of the EW corrections. Above the threshold its importance is reduced by the increase of the corrections of the gg-channel, which become dominant in the region µ ≥ 2500 GeV.

Differential distributions
In the previous subsection we studied the numerical impact of the next-to-leading order contributions to the total cross section for stop-anti-stop production at the LHC combining contributions of EW and QCD origin. Although the current experimental studies account for the higher order corrections by re-weighting the events generated at LO with the NLO(+NLL) predictions for the total rate, it is well known that higher-order contributions can significantly alter the shape of kinematic distributions.
In this subsection we study the impact of the EW corrections on differential distributions with respect to the invariant mass, the transverse momentum and the pseudo-rapidity, defined as respectively. The quantities p j , p T j and η j are the four-momentum, the transverse momentum and the pseudo-rapidity of the particle j, respectively. We have considered the six SUSY scenarios described in Table 1, but here we present results for the light-stop and for the tau-phobic scenario only, since the other scenarios exhibit similar features of the latter one. The results are collected in Fig. 9  Invariant mass distribution In scenarios with a heavy stop exemplified by the tauphobic scenario in Fig. 10, the relative impact of the EW corrections to the invariant mass distributions decreases as M inv increases. As can be inferred from Fig. 10(a) the EW corrections are positive and sizable in the small M inv region while in the large M inv region the EW corrections become negative. In this regime the relative impact of the corrections is almost flat and their size is small. The behavior of the EW corrections can be explained looking at the contributions of the different EW corrections collected in Fig. 10(b). The O(α s α) and the O(α 2 s α) corrections are the numerically dominant ones. In the small M inv region they are both positive and sum up enhancing the EW corrections. In the large M inv region, instead, the positive O(α s α) corrections partially compensate the large negative contributions of the O(α 2 s α) corrections. Also in the light stop-scenario the relative yield of the EW corrections decreases as M inv increases but, immediately above the production threshold, they become negative, see Fig. 9(a). As one can see in Fig. 9(b), this behavior is again related to mutual cancellation between the O(α s α) and the O(α 2 s α) corrections. In the entire M inv region the O(α 2 ) corrections are irrelevant for all considered scenarios.
Transverse momentum distribution The EW corrections to the transverse momentum distribution exhibit features similar to those for the invariant mass distribution.
In the low p T region the EW corrections are positive and sizable in all scenarios but the light-stop one, i.e. their relative yield is above 10%. For high values of p T , instead, the absolute value of the EW corrections decreases and they become negative, see Fig. 10(c). As shown in Fig. 10(d), this behavior is originated by the same interplay between the O(α s α) and O(α 2 s α) corrections present in the M inv distribution. As shown in Fig. 9(c), in the light-stop scenario the relative impact of the EW corrections is below 10%. This feature can be understood by looking at the behavior of the various EW contributions, Fig. 9(d). In the low p T region the O(α s α) and O(α 2 s α) contributions are both positive but small, i.e. below 5%, while in the high p T region the two contributions cancel against each other.
Pseudo-rapidity distribution As illustrated in Fig. 10(e), in all heavy-stop scenarios the EW corrections are important for large pseudo-rapidities, i.e. they are above 10% for |η| ≥ 2. In the low pseudo-rapidity region the EW corrections are positive but small, and almost negligible in the central region. The behavior of the EW corrections can be understood by looking at Fig. 10(f). The large and positive O(α s α) corrections originate from the gγ partonic process. This process proceeds via a t-and u-channel tree-level diagram, thus it produces preferably final particles in the forward region. In this region its large contribution is further enhanced by the positive O(α 2 s α) corrections, cfr. Fig. 10(f). In the low pseudo-rapidity region, instead, the positive contributions of O(α s α) are canceled by negative corrections of O(α 2 s α). As shown in Fig. 9(e), in the light-stop scenario the EW corrections are generally small, i.e. below 10% for |η| ≤ 4.
The importance of the EW corrections considerably depends on the kinematics of the producedt 1t * 1 . Thus, the yield on experimental event rates can be substantially affected by the application of kinematical cuts. Furthermore, NLO QCD corrections to the decay of the produced top-squarks may further alter significantly the shape of the kinematical distributions [57,58]. Therefore, such NLO corrections to the decay should also be investigated systematically for the EW corrections, including EW corrections also in the decay.

Conclusions
In this paper we have presented the first phenomenological study fort 1t * 1 production at the LHC including both the complete NLO QCD and NLO EW contributions. We have used the most recent PDF sets including QED effects, presenting a thorough study of parameter regions compatible with a SM-like Higgs boson observed at the LHC. The allowed MSSM regions are characterized by a rather heavyt 1 and/or by a large mixing in the stop sector, i.e. large X t .
Our analysis has shown that NLO EW contributions tot 1t * 1 production are not always negligible even on an inclusive level; they can be sizable, as large as 15 − 20% of the LO cross section, particularly in the parameter regions with heavy stops. This is mainly due to the contribution of the gγ-channel which increases with the mass of the produced topsquarks. This contribution can be easily included in the MSSM predictions for experimental investigations. Moreover, it could reduce the residual theoretical uncertainty related to the EW corrections below 5% of the LO cross section, at least in the parameter configurations where µ is not too large. The dependence on the remaining parameters of the model is found to be rather weak.
The presented study was performed using the Monte Carlo integrator SusyHell. This code includes the NLO EW corrections to all squark and gluino production channels and it will be made public in the future. It will allow for detailed phenomenological studies of the EW corrections to colored SUSY particle production, on an inclusive and also a fully differential level. Using this tool also the electroweak contributions should eventually be combined with higher order corrections to the decay processes.         Figure 6: Same as Fig. 4, but for the tau-phobic scenario. The value of the parameters not involved in the scans are collected in Table 1(d).    Table 1(e).      Figure 10: Same as Fig. 9 but for the tau-phobic scenario defined in Table 1(d).