Squark Production and Decay matched with Parton Showers at NLO

Extending previous work on the predictions for the production of supersymmetric (SUSY) particles at the LHC, we present the fully differential calculation of the next-to-leading order (NLO) SUSY-QCD corrections to the production of squark and squark-antisquark pairs of the first two generations. The NLO cross sections are combined with the subsequent decay of the final state (anti)squarks into the lightest neutralino and (anti)quark at NLO SUSY-QCD. No assumptions on the squark masses are made, and the various subchannels are taken into account independently. In order to obtain realistic predictions for differential distributions the fixed-order calculations have to be combined with parton showers. Making use of the Powheg method we have implemented our results in the Powheg-Box framework and interfaced the NLO calculation with the parton shower Monte Carlo programs Pythia6 and Herwig++. The code is publicly available and can be downloaded from the Powheg-Box webpage. The impact of the NLO corrections on the differential distributions is studied and parton shower effects are investigated for different benchmark scenarios.

GOLEM framework have been provided and compared to resummed predictions from jet merging. In [52][53][54], we have calculated the NLO corrections to the pair production of squarks of the first two generations and implemented the cross section in a fully flexible partonic Monte Carlo program without making any simplifying assumptions on the squark masses and treating the different subchannels individually. In the course of this calculation we have developed a new gauge-independent approach for the subtraction of on-shell intermediate gluinos at the fully differential level and compared our approach to several methods proposed in the literature. Moreover, we have extended the results [48,49,51], by matching our NLO calculation to parton showers using the Powheg-Box [55][56][57] framework.
These recent NLO calculations which take into account the full mass spectrum have shown that the K-factors of the individual subchannels can vary by up to 20%. Therefore, in order to improve the accuracy of the cross section predictions a proper NLO treatment of the individual subchannels is necessary, without relying on an averaged K-factor. Furthermore, it was found, that while the shapes of semi-inclusive distributions are only mildly affected by NLO corrections, this is not the case for more exclusive observables. Here the K-factors can vary by up to ±20% depending on the kinematics, both at the production level and after including squark decays supplemented by the clustering of partons to form jets. Irrespective of the use of fixed or dynamical scales, simply scaling LO distributions with a global K-factor is not a good approximation for exclusive observables.
In continuation of our effort to provide accurate predictions for SUSY production processes at the LHC we present in this work our results for the NLO SUSY-QCD corrections to squarkantisquark production of the first two generations. We furthermore combine our results both for squark pair production and for squark-antisquark production with the decay of the (anti)squark into the lightest neutralino and (anti)quark at NLO SUSY-QCD. All results are obtained at fully exclusive level and without making any simplifying assumptions on the squark mass spectrum. In order to obtain realistic predictions for exclusive observables we have combined our fixed-order NLO calculations with parton showers. To this end, the processes have been implemented in the Powheg-Box framework [54,57] and interfaced with different parton shower programs. The implementation has been made publicly available and can be obtained from [58].
The outline of the paper is as follows. In section 2 we present the NLO calculation of the squark-antisquark production process. Section 3 is devoted to the computation of the squark decays at NLO and the combination with the production processes. Here we study different approaches for the consistent combination at NLO. The implementation in the Powheg-Box as well as our results at fixed order and including parton shower effects are presented in section 4. Finally, we compare our results with results obtained with an approximate approach used in the SUSY searches by the LHC experiments. We summarize and conclude in section 5.

Squark-Antisquark Production at NLO
The calculation of the NLO corrections to squark-antisquark production is very similar to the one for squark pair production already presented in [52]. Therefore, the following discussion summarizes only the main steps and points out the most important differences between the two processes.

Contributing Channels
The production of a squark-antisquark pair at LO either proceeds via a pair of gluons or a quark-antiquark pair in the initial state: Here, the lower indices indicate the flavour of the particle, whereas the upper indices for the squarks denote the respective chirality. The contributing Feynman diagrams are depicted in Fig. 1. Due to the flavour conserving structure of the occurring vertices the gg initiated diagrams and the s-channel diagram contribute only to the production of squarks of the same flavour and chirality. The results for the individual matrix elements squared can be found in [54].
We consider in the following only the production of squarks of the first two generations mediated by the strong interaction. Correspondingly, the higher-order calculation comprises only SUSY-QCD corrections. In total, this leads to 64 possible final state combinations. This number can be reduced to 36 independent channels if the invariance under charge conjugation is taken into account. The number of independent channels can be reduced further if some of the squark masses are degenerate, as in this case the results for the qq initiated contributions differ only in the respective PDFs. However, we perform the calculation for a general mass spectrum and take advantage of this point only in the numerical analysis.

Virtual and Real Corrections
At NLO the squark-antisquark production processes receive contributions from virtual and real corrections. For the calculation of the virtual corrections we use the Mathematica packages FeynArts 3.8 [59][60][61] and FormCalc 6.1 [62,63]. The numerical evaluation of the loop integrals is performed with Looptools 2.7 [62,64].
In order to regularise the occurring ultraviolet (UV) divergences we apply Dimensional Regularisation (DR) [65][66][67][68][69]. The UV divergences are absorbed into the fields and parameters of the theory by introducing renormalization constants. For the renormalization of the strong coupling constant we use the MS scheme and decouple the heavy particles, i.e. the gluino, the top-quark and the squarks, from the running of the strong coupling constant α s . In the numerical analysis the 2-loop results for the determination of α s at the scale of the process are used, hence we require the 1-loop decoupling coefficient, which can be found e.g. in [21,[70][71][72]. Dimensional Regularisation violates SUSY explicitly by changing the number of degrees of freedom of the gluon field, inducing a mismatch between the gauge and the Yukawa couplings beyond LO. At NLO this effect can be cured by adding a finite SUSY restoring counterterm to the counterterm of the Yukawa coupling, see [73]. With these steps it is possible to use the five-flavour α (5),MS s in the numerical analysis. The occurring fields and masses are renormalized using on-shell renormalization conditions. As the relevant counterterms are not included in the FormCalc version we use, they had to be implemented by hand in the MSSM model file.
The actual calculation of the corrections is performed such that the full mass dependence is preserved. In principle this requires the generation of all possible production modes with FeynArts, which is obviously a very inefficient procedure. Instead, we generated only the virtual contributions for uū →ũ LūL , uū →ũ LūR , ud →ũ LdL , dd →ũ LūL and gg →ũ LūL , where the indices L and R refer to the left-and right-handed chirality of the squarks. All other combinations of squarks in the final state can be traced back to one of these cases. However, this procedure requires a generalization of the masses of the internal squarks, if the corresponding propagators are connected to an external squark or quark line. In case of squark pair production this step amounted to simply replacing all internal squark masses in the vertex and box corrections with the masses of the external squarks, while the self-energy corrections could be left unchanged. For squark-antisquark production this generalization is more involved and requires a dedicated consideration of the individual diagrams. Some sample graphs are depicted in Fig. 2. The first two diagrams in the upper row are examples for the case where all internal masses have to be kept, i.e. here no changes are necessary. In the next two graphs the masses of the squarks in the loop have to be replaced case by case according to the flavour of the initial state quarks. Note that both chiralities have to be taken into account. The diagrams depicted in the lower row of the figure are examples for the case where one or more internal squarks are connected directly or indirectly to the final state squarks. The masses in the corresponding propagators and loop integrals have to be generalized accordingly.
The real corrections consist of the contributions with one additional gluon in the final state: q iqj →q c1 kq c2 l g , g g →q c iq c i g . Moreover, at NLO a new channel occurs with a gluon and an (anti)quark in the initial state: These channels are related to each other by invariance under charge conjugation.
In order to calculate the q iqj , q i g and gq j channels it is sufficient to perform the calculation for one of them and construct the other combinations by either crossing the gluon or by charge conjugating the respective process. Here, the calculation is performed analytically for the q i g → q c1 kq c2 l q j subprocesses. The occurring traces are evaluated with FeynCalc 8.2 [74]. The calculation is performed using two gauges for the external gluon, the Feynman gauge and a light-cone gauge.
The gg-channels are obtained from MadGraph 5.1.3.1 [75] by generating the HELAS calls [76] for the specific process g g →ũ LūL g, generalizing the masses of the occurring squarks and removing the widths of the intermediate particles.
All these contributions exhibit infrared (IR) divergences, which cancel by virtue of the Kinoshita-Lee-Nauenberg theorem [77,78] against the corresponding divergences in the virtual contributions. As apt for a Monte Carlo event generator this cancellation is achieved by means of a subtraction formalism. We employ the FKS method [79], which is automated in the Powheg-Box.
In the qg-initiated channels q i g →q iqj q j a second type of singularity occurs for scenarios with mg > mq j . 3 For these mass configurations the intermediate gluino in the diagrams depicted in Fig. 3 can be produced on-shell, causing a resonant behaviour. A similar problem has already been encountered in the calculation of squark pair production [21,52]. Being formally equivalent to the Born contribution of on-shell squark-gluino production with the gluino decaying subsequently into a quark and an antisquark these contributions are large and require a proper definition of the process of interest. Keeping these terms would cause a double counting if the predictions for squark-antisquark production were combined with the ones for squark-gluino production. Hence, in order to obtain a meaningful result these on-shell contributions have to be subtracted consistently.
There exist several methods to cope with this type of singularities, which have been developed in the context of tW production [80], squark pair production [52,48] and squark/gluino production [21]. These approaches can be categorized as follows: • Diagram Removal (DR): In this approach the resonant contributions are removed by either completely neglecting the Feynman diagrams in Fig. 3 (DR-I) or by keeping the interference terms with the non-resonant contributions, but removing the amplitude squared of the depicted graphs (DR-II). Both approaches are rather easy to implement in a Monte Carlo program, but break gauge invariance.
• Diagram Subtraction (DS): These methods aim at a pointwise subtraction of the onshell contributions by constructing a counterterm and performing a suitable reshuffling of the momenta. Hence both the interference terms and the off-shell contributions are kept by construction. In order to regularise the singular behaviour for (pq j + p q j ) 2 → m 2 g a finite width Γg for the resonant gluino has to be introduced (in fact this is also required in the DR-II scheme in order to regularise the integrable singularity in the interference terms). In the original proposal for tW production [80] this is achieved by replacing the corresponding propagator: However, this approach is only gauge invariant in the limit Γg → 0. A fully gauge invariant modification of the DS scheme (denoted DS * in the following) has been proposed in the context of squark pair production [52]. In this approach the analytic expression for the amplitude squared is expanded in the poles (pq j + p q j ) 2 − m 2 g ≡ s jg before introducing the regularising width: The coefficients f k (k = 0, 1, 2) are gauge invariant quantities, i.e. introducing a regulator Γg at this point preserves gauge invariance and leads to The differences between the expressions obtained with the DS * and with the 'usual' DS method vanish for Γg → 0 as expected, see [52]. The counterterm for the subtraction of the on-shell contributions in this method is given by f 0 and reproduces the one used in the DS scheme in the limit (pq j + p q j ) 2 → m 2 g . For more details on the momentum reshuffling and the construction of the subtraction term see [52].
The comparison of these different subtraction methods for squark pair production revealed for the scenario considered in [52] only discrepancies in the total cross section at the per mille level. Repeating this study for squark-antisquark production, however, leads to larger differences, as the contributions of the qg initiated channels are larger in this case. To illustrate this point the predictions for the total production cross sections of squarks of the first generation as obtained with the DR-II (using the light-cone gauge) and the DS * scheme are summarized in the second and third column of Tab. 1. The scenario considered here corresponds to the mSUGRA point 10.4.5 [81] specified in Sec. 4. For the regularising width we choose Γg = 1 GeV. As can be inferred from the percental difference between the respective numbers given in the fourth column the predictions obtained with these two methods differ by up to 11%, leading to a discrepancy of 2.7% after summing these channels. Taking into account the contributions of the squarks of the second generation, too, increases this discrepancy further: The charge conjugate channels are combined. The last four columns contain the numerical values for the quantity σqg as defined in the text and the respective contribution to the full NLO cross section, again for both the DS * and the DR-II method.
corresponding to a discrepancy of 3.6%.
The fifth and seventh column of Tab. 1 contain the respective predictions σ qg for the qg contribution to each channel. This (unphysical) quantity comprises the 2 → 3 parts of the respective channel, i.e. the real amplitudes squared and the corresponding FKS counterterm and hence allows for a direct estimation of the effects of the applied subtraction scheme. As can be inferred from the table these contributions make up several percent of the individual cross sections. Hence the large discrepancies observed between the two subtraction methods have significant effects on the predictions for the total cross sections.
Even larger effects of the chosen subtraction scheme can be observed in differential distributions which are sensitive to the emitted parton of the real corrections. As an example, the p T distribution of the radiated parton obtained with the DR-II scheme and the DS * method is shown in Fig. 4 (left). For p j T > 200 GeV the two predictions differ by about 30%. In contrast, the shape of the mqq distribution (right plot in Fig. 4), which is supposed to be less sensitive to additional radiation, is not affected by the chosen method. Solely the normalization reflects the 3.6% discrepancy already encountered in the total cross section

Tests and Comparison
The calculation presented in the last section has undergone numerous checks and comparisons. An obvious test for the correctness of the calculation consists in a comparison with the public program Prospino2 for the limit of a mass degenerate spectrum. Unfortunately, a direct comparison of the results obtained with this public code is not straightforward, as it implicitly takes into account the sbottom production processes gg →bb and qq →bb, while the contributions bb →bb are neglected. Moreover, at NLO the contributions qg →qbb and the charge conjugate processes are taken into account. Instead of mimicking the way the total K-factor is calculated in Prospino2 we have compared the numerical results of our calculation with a non-public implementation of the original results from [21], denoted Prospino * in the following. Besides testing our calculation for the special case of degenerate squarks we have intensively checked the individual building blocks: • The Born expressions have been compared with results given in the literature [16][17][18]. In addition, the numerical comparison of the total cross section with Prospino * provides a simple cross check for the correctness of the nontrivial combinatorics of the contributing channels.
• The UV finiteness of the virtual corrections has been checked both analytically and numerically. The correct structure of the IR poles has been verified by comparison with the known structure for the case of massive coloured particles in the final state, see e.g. [82]. The correctness of the modifications performed in the virtual routines in order to generalize them to an arbitrary mass spectrum has been tested by performing this generalization for both gg →ũ LūL and gg →d LdL and comparing the outcome numerically. Likewise, the other cases mentioned in Sec. 2.2 have been checked.
• The analytic results for the real matrix elements squared have been compared numerically for a multitude of arbitrary phase space points with the routines generated with Mad-Graph 5. The cancellation of the IR poles against the FKS counterterms has been tested using the automatic procedure provided by the Powheg-Box. The gauge invariance of the DS * scheme has been explicitly checked by comparing the outcome of the two different gauges used in the calculation. Furthermore, the equivalence of the DS and the DS * scheme in the limit Γg → 0 has been verified numerically.
• The individual results for the three production channels gg, qq and qg have been compared for degenerate mass spectra with Prospino * .

Squark Decays at NLO and Combination with Production Processes
The calculation of NLO SUSY-QCD corrections to production processes is only a first step towards a realistic prediction of possible events at the LHC. The next step requires the inclusion of the decays of the produced particles. Here, only the decay mode into a quark and the lightest neutralino for squarks,q → qχ 0 1 , or into an antiquark and the lightest neutralino for antisquarks, q →q +χ 0 1 , will be taken into account. In many SUSY scenarios, in particular the ones studied in this paper, the lightest neutralino is the lightest supersymmetric particle and stable (if R-parity conservation is assumed). The SUSY-QCD corrections to this decay have been known for several years. However, in the original calculation [83] only results for the partial width have been given, so that a differential description is not possible. Recently, a fully differential calculation in the context of squark pair production and decay has been presented [48] where the radiative corrections to the decay have been included by using the phase-space-slicing technique. In the next subsection we present a recalculation of the decay at the fully differential level by applying the subtraction method developed for single top production and decay [84] to our process. The second part of this section deals with the consistent calculation of the total squark width, which is required for the combination of the production and decay processes described in the last part of this section.

Decay Width forq
The LO contribution to the decay width of the process comprises only one Feynman diagram, which is depicted in Fig. 5 (a). As in the production process, at NLO virtual and real corrections have to be taken into account. The two diagrams contributing to the virtual corrections are shown in Fig. 5 (b). The calculation is performed in DR and the external fields are renormalized on-shell. All integrals have been evaluated analytically. The package HypExp [85,86] has been used to expand hypergeometric functions. Using DR requires, again, the introduction of a finite counterterm. Here, the squark-quark-neutralino Yukawa couplingĝ is affected. The SUSY-restoring counterterm leads to the following relation to the gauge coupling g [73],ĝ The real corrections involve an additional gluon, emitted either from the squark in the initial state or the quark in the final state, as displayed in Fig. 5 (c). The IR divergences arising from the soft and/or collinear emission of the gluon cancel against the corresponding ones in the virtual corrections. This cancellation is achieved on the differential level by applying the subtraction method developed in [84] for single top production and decay to our decay process. The divergences in the real radiation process are cancelled by a local counterterm, which is constructed such that it has the same singular behaviour as the full matrix element. It takes the form of the LO matrix element squared multiplied by a function D, which describes the emission of the soft or collinear radiation: In the limit of soft emission, when p g → 0, or where the momenta of the quark p q and the gluon p g are collinear, the counterterm on the right-hand side of Eq. (11) has the same singular structure as the full matrix element squared on the left-hand side. The LO matrix element M 0 is evaluated with modified momenta p q and p χ 0 1 which absorb the momentum carried away by the gluon. They are subject to momentum conservation pq i = p q + p χ 0 1 , as well as to the on-shell conditions p 2 q = 0 and p 2 In the process at hand advantage can be taken from the fact that the LO matrix element squared can be easily factorized from the divergent part of the real matrix element squared, with the function f (y, z), calculated in d = 4 − 2 dimensions, defined as In this function the following substitutions have been made with r denoting the squared ratio of the neutralino over the squark mass, The coefficient of the LO matrix element squared can then be chosen to serve as the divergent part, denoted by D, in the counterterm in Eq. (11), In order to cancel the IR divergences in the virtual corrections this counterterm needs to be integrated analytically over the one-particle phase space of the emitted gluon. The results for the necessary integrals can be found in Table 1 of [84]. The integrated counterterm then reads . (18) All steps of the analytical calculation have been checked against [87]. The results for the partial widths at LO and NLO have been compared to the result obtained from Sdecay 1.3 [88]. Moreover, we have compared our results to the results presented in the independent calculation of squark pair production and decay of [48], in particular to the results given in Table 6 for the benchmark point 10.1.5 and the corresponding distributions, and have found agreement. In addition, this decay has been implemented in the Powheg-Box. The virtual corrections for this independent calculation have been calculated with FeynArts/FormCalc and the loop integrals have been evaluated with LoopTools. The real matrix elements squared, calculated by hand, have been tested numerically for a multitude of phase space points against the corresponding routines obtained with MadGraph. In the Powheg-Box the cancellation of the divergences is achieved automatically via the implemented FKS subtraction method. We have found perfect agreement between the calculation presented in this section and the implementation in the Powheg-Box.

Total Squark Width at NLO
For the calculation of the squark branching ratios we also need the total decay width Γq tot , both at LO and NLO. Furthermore, the NLO total decay width will be necessary to normalize the expressions for the combination of the production and decay processes, as we will see in the next subsection. Since we only consider the decay into a quark and the lightest neutralino as possible 'decay chain' for the produced squarks, it is not necessary to describe all other partial decay widths differentially. Therefore, they can be extracted from the literature or their implementation in Sdecay. In order to implement the various decay routines from Sdecay in our code the following adaptions had to be made for the individual decay modes: The decays into neutralinosχ 0 k or charginosχ ± l are mediated by electroweak interactions. The decay into charginos is only possible for left-chiral squarks. In the routines for the (N)LO results [83] taken from Sdecay only the conventions for the parameters, especially those entering the calculation of the squark-quark-gaugino vertex, had to be adapted. In our calculation the weak mixing angle θ W is determined according to Eq. (10.11) from [89], yielding All other parameters needed for the numerical evaluation of the decay widths can be found in Sec. 4.2.
• Strong decay:q i → q ig The NLO corrections to this strong decay mode have been calculated in [90]. However, this calculation has been performed for degenerate squark masses and implemented in the same way in Sdecay. In order to incorporate the full mass dependence, we have calculated the gluino self energy using the corresponding function from the calculation of the stop decays [91]. In these decays the correctt 1,2 masses, the top quark mass and thet mixing angles have been used. For each squark of the first two generations this function is called by replacing the appropriate squark mass and setting the quark mass and mixing angles to zero.
Also in the calculation of α s , where the heavy particles are decoupled from the running, the squark masses are assumed to be degenerate. To restore the full mass dependence in α s , the logarithms of the masses of the heavy, decoupled particles have been modified to obtain the logarithms given in Eq.
(3) of [52]. In Sdecay the strong coupling constant is converted from the MS scheme, used in the original calculation, to the DR scheme. In order to use α s as implemented in the Powheg-Box we calculate α s in the MS scheme by omitting the conversion factor introduced in Sdecay.

Combination with the Production Processes
A consistent combination of the production processes at NLO with the subsequent decays of the squarks,q → q +χ 0 1 , or antisquarks,q →q +χ 0 1 , at NLO is the next necessary step. In this combination we take into account only those contributions to the process pp → 2q +2χ 0 1 that lead to two on-shell intermediate squarks. In the narrow width approximation, which is valid in the scenarios analysed here since the widths of the squarks fulfil Γq i /mq i 1, the differential cross section factorizes into the production cross section times the branching ratios of both squark decays By applying the narrow width approximation we not only neglect contributions with off-shell squarks, which are known to be suppressed by Γq i /mq i , but also non-factorizable higher-order contributions. The latter ones comprise interactions between particles of the production and decay stage or between final-state particles of the two decays. These contributions are expected to be suppressed by Γq i /mq i as well [92,93]. Only long-range interactions induced by the exchange of soft gluons could still affect the results of exclusive observables. However, an analysis of these effects is beyond the scope of this publication. Aiming at a combination of the decays at NLO with the production process at NLO the factors in Eq. (20) have to be replaced by the NLO quantities: This expression obviously includes beyond-NLO contributions. In order to strictly consider NLO accuracy it has to be expanded to NLO in α s . There exist two approaches for this problem, both developed in the context of single and pair production of top quarks [84,94].
In the first approach a Taylor expansion of the full expression is performed. This leads to a formula which is normalized to the LO total widths and subtracts the ratios of the NLO corrections to the total widths over the LO total widths from the first term: This subtracted term might lead to negative contributions, if the NLO corrections to the total width are positive and large while the corrections to the partial widths are small. However, this expansion has the advantage that the sum over all possible decay channels reproduces the production cross section, i.e. the branching ratios of all subchannels add up to one.
In the second approach only the numerator is expanded in α s while the NLO total widths are kept in the denominator. This expansion avoids the problem of potentially negative contributions and leads to the expression: In this approach summing over all possible decay channels does not reproduce the production cross section, as the branching ratios do not add up to one, and in this sense unitarity is violated.
Since both expansions to NLO accuracy may cause problems the complete expression in Eq. (21) can be used as an alternative approach. On the one hand, in this approach the branching ratios sum up to one, but on the other hand only parts of the possible beyond-NLO corrections are included. Given a good convergence of the perturbative series we expect these terms to be small, however.
In chapter 4 results for all three possible combinations of production and decays at NLO, according to Eqs. (21)-(23), will be presented and compared. However, in all other results the Taylor expansion of the cross section, Eq. (22), will be used since in this approach the unitarity of branching ratios is preserved. In the scenarios analysed here the subtracted terms in Eq. (22) are unproblematic, i.e. the NLO corrections to the total decay widths are small (see Tab. 7).

Implementation and Results
After a brief discussion of the steps required for the implementation of squark production and decay in the Powheg-Box this section summarizes our main findings, including both numerical results at fixed order perturbation theory and after application of different parton showers. Moreover, we present some results for total rates after applying realistic experimental search cuts.

Implementation in the Powheg-Box
The implementation of squark-antisquark production in the Powheg-Box is essentially identical to the case of squark pair production, which has been extensively discussed in [52]. Besides several changes in the main code required for the consideration of processes with strongly interacting SUSY particles the process-dependent parts have to be provided, which comprise • all independent flavour structures contributing to the Born and real channels, • the Born and the spin/colour-correlated matrix elements squared, • the finite part of the virtual contributions, • the real matrix elements squared • and the colour flows for the Born configurations.
The implementation of the various subtraction schemes discussed in Sec. 2.2 is rather involved and has been described in detail in [52], too. In essence, we have implemented (besides the two DR schemes) several versions of the DS scheme by splitting the real matrix element squared into a part containing the resonant gluino contributions and the corresponding subtraction terms and a part containing all other terms. The resonant parts do not contain any IR singularities and can therefore be treated independently from the Powheg-like event generation, similar to the 'hard' part R h of the real matrix elements squared introduced below.
We have implemented these building blocks for squark-antisquark production into the version 2 (V2) of the Powheg-Box and ported our previous implementation of squark pair production to the V2. This newer version of the Powheg-Box allows for the consideration of NLO corrections to the decays of the on-shell produced squarks. We use this new option to combine our results for the NLO production processes with the corrections to the specific decayq → qχ 0 1 described in the previous section. Besides taking into account the decay products in the flavour structures as described in the manual of the Powheg-Box V2 this requires the combination of the production and decay matrix elements according to the combination formula in Eq. (22). Moreover, the FKS subtraction of the IR divergences related to the gluon emission off either the squark or the quark in the NLO corrections to the decay process requires the specification of the colour correlated Born matrix elements squared. These are trivial in the case at hand and read in the convention of the Powheg-Box Bq q = C F B.
In order to check the correctness of the implemented results the same tests as described in [52] have been performed. These comprise a comparison of numerous differential distributions evaluated at NLO with the corresponding results after generation of the hardest emission according to the Powheg method, both at the level of the production processes and after including the decays. While we find an excellent agreement for inclusive quantities, a strong enhancement of the Powheg results compared to the respective NLO distributions is observed for exclusive quantities like the transverse momentum of the squark-antisquark system, pqq T . The same artificial enhancement has already been observed in case of squark pair production and can be cured by using the soft/collinear limits of the real matrix elements squared R instead of the full expressions for the generation of the hardest radiation. In the Powheg-Box this is achieved by introducing a function F which separates the soft/collinear part R s and the hard part R h of the real matrix elements squared: This function F has to fulfil F → 1 in the soft/collinear limit and should vanish far away from the corresponding phase space regions. In the Powheg-Box the functional form 10.3.6 * 825 GeV 550 GeV 0 GeV 10 +1 10.4.5 1150 GeV 690 GeV 0 GeV 10 +1 is used, with the transverse momentum p T of the emitted parton with respect to the emitter and a damping parameter h (see [55] and [57] for further details). Similar to our earlier studies on squark pair production we use h = 50 GeV throughout. This choice was found to damp the artificial enhancement in the pqq T distribution and reproduces the NLO prediction for pqq T 200 GeV, while maintaining the Sudakov damping for small transverse momenta inherent in the Powheg method.

Setup
For the numerical analysis we consider two mSUGRA scenarios which are not yet excluded by data, see e.g. [95,96]. The scenarios are based on the CMSSM points 10.3.6 * 4 and 10.4.5 from [81]. The input parameters of these scenarios are summarized in Tab. 2. The mass spectrum of the SUSY particles has been generated with Softsusy 3.3.4 [97], the resulting on-shell masses are then used as input parameters. For the SM parameters the following values are used [89]: As Softsusy implements non-vanishing Yukawa corrections, there is a small difference between the masses of the second-generation squarks and the corresponding first-generation ones, i.e. mũ L = mc L etc. To simplify the analysis and save computing time these masses are replaced by the mean of the mass pairs, i.e. mũ L and mc L are replaced by (mũ L + mc L )/2 and so on. The obtained masses for the squarks of the first two generations and the gluino masses are summarized in Tab. 3. Note that for the point 10.3.6 * the mass hierarchy is mq > mg, while for 10.4.5 mq < mg, the latter point requiring the subtraction of contributions with on-shell intermediate gluinos as described in Sec. 2.2. Here, the DS * method is used, with a default value for the regulator Γg = 1 GeV (recall that this regulator is only needed if a subtraction is required, thus in all other cases it is set to zero).  Furthermore, the partial and total decay widths of the squarks depend on the masses of the charginos and neutralinos and the respective mixing matrices. The masses of the neutralinos  and charginos for the two scenarios are given in Tab In order to diagonalize the chargino mass matrix two matrices are needed, one for the left-handed components (denoted U ) and one for the right-handed ones (denoted V ). These 2 × 2 mixing matrices are parametrized as (i = U, V ) The mixing angles are given by cos θ U = 0.97213 and cos θ V = 0.98594 for the parameter point The PDFs are taken from the LHAPDF package [98]. For the LO results shown in the following the LO set CTEQ6L1 [99] with α s (m Z ) = 0.130 is used, while the NLO results are calculated with the NLO set CT10NLO with α s (m Z ) = 0.118 [100]. The strong coupling constant for the LO results is correspondingly computed using the one-loop renormalization group equations (RGEs), while the value used in the NLO results is obtained from the two-loop equations.
All results are calculated for the LHC with √ s = 14 TeV. The error bars shown in the following represent the statistical errors of the Monte Carlo integration.
Taking into account the decays of the produced squarks into qχ 0 1 or applying a parton shower algorithm leads to a potentially large number of partons in the final state. These partons are clustered into jets with Fastjet 3.0.3 [101,102]. To this end the anti-k T algorithm [  adopted, using R = 0.4. In the following only minimal cuts are applied on the transverse momentum and the pseudorapidity of the resulting jets: p j T > 20 GeV and |η j | < 2.8 .
Except for the results shown in Sec. 4.3.3 no event selection cuts are imposed.

Results at Fixed Order
The first part of this section is devoted to a discussion of the NLO corrections to squarkantisquark production. In the second part we present some results for the combination of production and decay, both for squark-antisquark and squark pair production. Hence this part extends our previous results for the squark pair production processes in [52] by also including the NLO corrections to the decay.

Squark-Antisquark Production
The results for the total squark-antisquark production cross sections determined at LO and NLO for the two benchmark scenarios defined in Sec. 4.2 are summarized in Tab. 5. In order to assess the theoretical uncertainties we vary the renormalization and factorization scales by a factor two around the central value µ = mq. The resulting percental uncertainties are also given in the table. Considering the resulting K-factors we note that in both cases the SUSY-QCD NLO corrections are positive and large, resulting in K ≡ σ NLO /σ LO ≈ 1.4. The scale uncertainties are strongly reduced by taking into account the NLO corrections, as expected. Turning next to the individual K-factors for the subchannels contributing to squark-antisquark production we observe that they differ significantly from the total K-factor obtained after summing the cross sections for all individual channels. To illustrate this point, the LO/NLO cross sections and the resulting K-factors for the production channels involving only squarks of the first generation are given in Tab. 6 for the CMSSM point 10.3.6 * . Note, that the channels with squarks of the same flavour and chirality in the final state, displayed in the first four rows of the table, have contributions from gg initial states and therefore larger K-factors than channels with squarks of different flavour or chirality. Hence the assumption that the individual K-factors can be approximated by the total K-factor obtained from Prospino is in general not valid. Determining the individual corrections consistently is especially important if the decays are taken into account and the branching ratios of the different squarks differ significantly for the specific decay channel under consideration. In order to assess the possible numerical impact of this approximation we consider the decayq → qχ 0 1 at LO at the level of total cross sections, i.e. we multiply the production cross sections for the individual squark-antisquark production channels with the respective LO  branching ratios. In this step we take into account the contributions of the second generation squarks as well.
We first consider the benchmark scenario 10.3.6 * . Using the correctly calculated NLO results for the individual production channels, multiplying them with the corresponding branching ratios and summing all channels, we obtain channels σ NLO · BR LO q →χ 0 1 q · BR LO q →χ 0 1q = 0.139 fb.
To mimic the way Prospino obtains the individual NLO results a common K-factor has to be calculated, using an averaged squark mass mq = 1779.31 GeV. In the case at hand this leads to σ avg LO = 2.315 fb , σ avg NLO = 3.218 fb ⇒ K avg = 1.39 .
Note that the difference compared to the full calculation given in Tab. 5 is marginal and not visible when rounding to the second decimal place. This is due to the fact that the spread in the squark masses is rather small. Multiplying the LO result for each subchannel with this common K-factor and the corresponding branching ratios gives Thus the rate obtained with the approximation relying on a constant K-factor for all subchannels is roughly 10% smaller for this special case. Repeating this procedure for the benchmark scenario 10.4.5 one obtains for the Prospinolike K-factor σ avg LO = 3.090 fb , σ avg NLO = 4.356 fb ⇒ K avg = 1.41 .
Again, comparing this result to the full calculation given in Tab. 5 the discrepancy is only marginal.

Squark Production and Decay at NLO
As discussed in Sec. 3.3 we have used three different approaches to combine the production and decay processes at NLO, differing in the way the combined expression is expanded in α s . All approaches require the calculation of the total squark width, either at LO or NLO accuracy. The results for the two considered benchmark scenarios are summarized in Tab  In a first step we compare the numerical results obtained with these approaches, both for differential distributions and total cross sections. In Fig. 6 the distributions for the transverse momenta of the hardest and the second hardest jet, p j 1 /j 2 T , their invariant mass m j 1 j 2 and the missing transverse energy / E T are depicted for squark pair production using the benchmark scenario 10.3.6 * . Here, App. 1 corresponds to the Taylor expansion according to Eq. (22), whereas in App. 2 only the numerator in the combination formula is expanded, see Eq. (23). Approach 3 is the result obtained without any expansion, i.e. these distributions include contributions which are formally of beyond-NLO.
The discrepancies between the approaches 1 and 2 can amount to up to O(15%) for the jet distributions and are largest close to threshold, while the results for / E T reflect only the overall discrepancy in the total cross sections, which amounts for this scenario to approximately 4%. The distributions obtained with the third approach do not show any large deviations from the results obtained with the other two approaches, but suffer from large statistical fluctuations. These result from the more complicated structure of the phase space integrations.
The total cross sections for the combined production and decay processes as obtained with the approaches 1 and 2 are summarized in Tab. 8, both for the scenario 10.3.6 * and 10.4.5. Note that the predictions for the LO cross sections are identical in both approaches and have been calculated according to Eq. (20) using the LO quantities. Comparing the results for the different predictions at NLO reveals only rather small discrepancies < 4% for the total rates for the scenarios considered here. In the rest of this chapter we use exclusively the first approach to combine production and decay processes. In order to assess the influence of the NLO corrections on differential cross sections we consider in the following the differential K-factors for several observables. In Fig. 7 the LO and the NLO distributions for the transverse momentum of the hardest jet, p j 1 T , its rapidity, y j 1 , the missing transverse energy / E T and the effective mass m eff ≡ p j 1 T + p j 2 T + / E T are depicted for squark-antisquark production, using the benchmark scenario 10.3.6 * . The results for the scenario 10.4.5 are qualitatively the same. Considering the p T distribution of the hardest jet one observes a strong enhancement of the NLO corrections for small values of p T , while they turn even negative for large values. The result for the second hardest jet, which is not shown here, is qualitatively the same. A similar observation holds for the effective mass: the NLO curve is dragged to smaller values of m eff and the differential K-factor depicted in the lower panel is far from being flat over the whole region. For the / E T predictions, in contrast, the deviation of the differential K-factor from the total one is rather small, of O(5%), except for events with very small or very large missing transverse energy. Likewise, the shape of the rapidity distribution of the hardest jet is hardly affected by the NLO corrections. Next we consider the same set of qq -App. 1 9.29 · 10 −2 1.03 · 10 −1 1.11 5.73 · 10 −1 9.15 · 10 −1 1.60 qq -App. 2 9.29 · 10 −2 9.88 · 10 −2 1.06 5.73 · 10 −1 9.32 · 10 −1 1.63 distributions for squark pair production with subsequent decays, this time for the scenario 10.4.5, depicted in Fig. 8. Again the results for 10.3.6 * are qualitatively identical and not shown here.
In essence, the behaviour is very much the same as for squark-antisquark production in Fig. 7 and differs only in details. For example the differential K-factor of the rapidity distribution y j 1 shows slightly larger deviations from the total K-factor, whereas the one for / E T is a bit flatter in the range considered here.

Parton Shower Effects
In order to investigate parton shower effects we have combined our implementations of the squark production and decay processes with different parton shower programs. To this end five million events have been generated for squark-antisquark and squark pair production for each of the two benchmark scenarios defined in Sec. 4.2. The results shown in the following have been obtained by setting the folding parameters of the Powheg-Box to the values n ξ = 5, n y = 5, n φ = 1 , reducing such the number of events with negative weights. However, in the context of squark production and decay processes two further sources of negative events can occur. The first one originates from the way production and decay are combined in Eq. (22), see the discussion in Sec. 3.3. It is not possible to apply the folding procedure described above in this case, since the negative contributions to B are directly related to the (modified) Born contribution. Using a different expansion of the combination formula, e.g. Eq. (23), would remedy this point, however this approach violates unitarity and should therefore be avoided. The implemented subtraction schemes described in Sec. 2.2 present another source of contributions with negative weights. While these are completely absent for the DR-I method and their number can be reduced again by means of folding for the DS * and the DR-II method, they inevitably occur for the methods relying on a splitting of R.
All in all, using the (for conceptual reasons preferable) DS * subtraction scheme with split real matrix elements squared and Eq. (22) for the combination of production and decay leads unavoidably to events with negative weights, which cannot be neglected. Therefore, they are kept in the generated event files by setting the Powheg-Box flag withnegweights=1.
The generated event samples have been showered with two Monte Carlo event generators, using three different parton shower algorithms implemented in these programs: • Pythia 6: We use the version 6.4.28 [104]. All results have been obtained with the Perugia 0 tune [105], invoked by setting MSTP(5) = 320. A comparison with the Perugia  (22). In all plots the lower panels depict the respective differential K-factor (full) and the total K-factor from Tab. 8 (dashed). 11 tune (MSTP(5) = 350) yields only tiny discrepancies. 5 In order to study only effects of the parton shower, hadronization and multi-parton interaction (MPI) effects have been turned off by setting MSTP(111) = 0 and MSTP(81) = 20, invoking thus the use of the p T -ordered shower.
However, in the simulation of the full process, including NLO corrections to the production and the decays, a further subtle difficulty arises when using Pythia, which is related to the way the starting scales for the shower are chosen. The Powheg approach relies on the assumption that the p T of the emitted final-state parton is larger than the transverse momentum of any subsequent splitting generated by the parton shower. This requires the application of a p T veto in the parton shower, with the maximal scale being read for each event from the event file. However, if final-state resonances are present the mass of those has to be preserved by the reshuffling operations performed in the shower algorithm. Therefore, the showering of partons originating from the decays of these resonances, i.e. in the processes considered here the produced squarks, is performed separately in Pythia. The starting scale for these shower contributions is set to the invariant mass of all decay products, hence in the case at hand to the mass of the respective squark. In the scenarios considered here this scale is typically an order of magnitude larger than the upper scale written to the event file, leading to much more radiation and thus to a strong bias of the results. In order to correct for this effect, the Pythia routines had to be adapted to use the scale specified in the event file as starting scale in all individual contributions to the parton shower.
• Herwig++: The default shower of Herwig++ [106] is ordered in the angles of the branchings. Applying this shower to an event sample generated according to the Powheg method requires again the use of a p T veto. However, this combination lacks the emission of soft wide-angle partons, as the first emission in an angular-ordered shower is not necessarily the hardest one. In principle these missing parts have to be simulated in an extra step via a so-called vetoed truncated shower, which is not provided by Herwig++ and thus not taken into account in the following. The effect of this missing part will be estimated by comparing the results to those obtained with the p T -ordered Dipole-Shower [107,108], which is also part of the Herwig++ framework. The results presented in the following sections have been obtained using the version 2.6.1 [109]. In the following Herwig++ refers only to the default shower, while the results labeled Dipole-Shower or, for the sake of brevity, Dipole refer to the Dipole-Shower included in the Herwig++ framework.
The showered results for squark-antisquark production are shown in Fig. 9, using the scenario 10.3.6 * . Likewise, Fig. 10 depicts the results for squark pair production, obtained with the scenario 10.4.5. All plots show the outcome of the three parton showers described above and the NLO prediction, which serves as normalization in the ratio plots shown in the lower panels. The results for squark pair production using scenario 10.3.6 * and squark-antisquark production with scenario 10.4.5 do not reveal any new features compared to the depicted combinations and are therefore not shown here.
Comparing the predictions for the individual observables shown in the two figures we note that in all cases considered here the p j 1 T result obtained with Herwig++ is in the low-p T region slightly (O(10%)) enhanced compared to the other parton showers, whereas the Dipole-Shower and Pythia essentially agree here. At the other end of the spectrum, however, both the Herwig++ and the Dipole-Shower predict O(10%) smaller rates than Pythia, which is almost in accordance with the NLO result for large values of p j 1 T . The outcome of Herwig++ and the Dipole-Shower is identical in this kinematic regime. Similar conclusions can be drawn from the p j 2 T , the m j 1 j 2 and the m eff distributions not shown here. In contrast, the distributions describing the third hardest jet show more pronounced differences. Comparing first the results for p j 3 T obtained with Herwig++ and the Dipole-Shower one notices that they agree within O(5 − 10%). The result for the third jet obtained with Pythia is in all cases smaller compared to the other two parton showers. While the discrepancy using the benchmark scenario 10.3.6 * is for both squark-antisquark and squark pair production smaller than 10%, it amounts to 10-15% for the scenario 10.4.5 in both cases. The largest differences in the three shower predictions emerge in the results for the pseudorapidity of the third hardest jet, η j 3 . While Pythia and the Dipole-Shower agree within 5% for all cases and differ in case of squark pair production only in the overall normalization, but not in the shape of the distributions, Herwig++ predicts evidently more jets in the central region η j 3 1. Comparing the Herwig++ result and the Pythia outcome for squark-antisquark production, this enhancement amounts to a 20% higher rate in the centre and a reduction of the same magnitude for η j 3 ≈ 2.8. In case of squark pair production, this effect is smaller, of O(10%), but still clearly visible. The predictions for the missing transverse energy / E T agree very well and essentially reproduce the NLO result. Tiny deviations are only visible in the tails of the distributions, however they are smaller than 5% in all cases.
All in all, the predictions of the different parton showers for the observables depending solely on the two hardest jets agree within O(10%) or better. Comparing the showered results with the outcome of a pure NLO simulation the effects of the parton showers on these observables are at most of O(10−20%), except for the threshold region. By and large, the two Herwig++ showers yield larger deviations from the NLO outcome for these observables, whereas Pythia reproduces the NLO curves within O(10%). The / E T distribution is in all cases hardly affected by parton shower effects. Larger deviations between the different parton showers emerge in the predictions for the third hardest jet, which is formally described only at LO in the hard process. Especially the Herwig++ prediction differs significantly from the other two showers and predicts more jets in the central region of the detector. At this point it is not possible to decide ultimately if this discrepancy is an effect of the missing truncated shower or simply a relict of the way the phase space is populated in the different shower algorithms. This would require the actual implementation of a vetoed truncated shower, which is beyond the scope of this work. However, comparing the outcomes of the Dipole-Shower and Herwig++ reveals only very small discrepancies in  Table 9: Total cross sections after applying the event selection cuts defined in Eq. (37) for the different production modes in the two benchmark scenarios. The decays of the squarks (antisquarks) to qχ 0 1 (qχ 0 1 ) are included at NLO. The given results have been obtained at the level of a pure NLO simulation and including parton shower effects with Pythia and Herwig++, respectively. The last two rows have been obtained by rescaling LO events after application of the Pythia/Herwig++ shower with the constant K-factor and the individual NLO branching ratios.
other observables. Hence the overall effect of the neglected truncated shower seems to be small.

Total Rates
In the last step the created event samples are analysed using a realistic set of event selection cuts, which corresponds to the definition of the signal region 'A-loose' for the SUSY searches in two-jet events performed by the ATLAS collaboration [95].
The event selection cuts used in this analysis are Here, the effective mass m eff is defined as the sum of the p T of the two hardest jets and / E T , whereas the inclusive definition of this observable includes all jets with p j T > 40 GeV, Moreover, ∆φ(j i , / E T ) denotes the minimal azimuthal separation between the direction of the missing transverse energy, / E T , and the i th jet. The additional cut ∆φ(j 3 , / E T ) > 0.4 is only applied if a third jet with p j 3 T > 40 GeV is present. Applying these cuts at the level of a pure NLO simulation yields for the total cross sections for squark (anti)squark production combined with the subsequent decays in the two benchmark scenarios 10.3.6 * and 10.4.5 the results given in the first row of Tab. 9. Matching these NLO results with a parton shower hardly affects the outcome after using the cuts defined in Eq. (37), as can be inferred from the results obtained with Pythia and the Herwig++ default shower listed in the second and third row, respectively. Note that due to the mixture of cuts on inclusive and exclusive quantities the rates predicted by the two showers are slightly larger compared to the NLO case. Moreover, the two parton showers yield identical rates within 1-2%.
In order to compare these results obtained with our new calculations and implementations with the values determined according to the setup used so far for the simulation of these processes we proceed as follows: first the production and decays of the squarks are simulated with LO accuracy. The resulting events are reweighted with a common K-factor for squark-antisquark or squark pair production, which is obtained from Prospino, i.e. assuming degenerate squark masses and averaging over all channels. Each individual production channel is then multiplied with the corresponding NLO branching ratios for the produced squarks. The rescaled events are subsequently processed with the Pythia and the Herwig++ default shower, neglecting again effects of hadronization, MPI, etc. The results obtained with this approximate setup after applying the event selection cuts defined in Eq. (37) are summarized in the last two rows of Tab. 9. Comparing these total rates with those obtained in the full simulation one notes that the discrepancy is almost negligible in case of squark pair production, but amounts to 15-20% for squark-antisquark production. This discrepancy is mainly caused by assuming a common K-factor for all subchannels instead of using the exact results with individual K-factors when combining production with decay. This effect in squark-antisquark production has already been demonstrated in section 4.3.1 for the case of LO decays. In squark pair production, however, subchannels with K-factors close to the global K-factor have large branching ratios and therefore the exact and the approximate method give similar results. These examples illustrate that in order to obtain precise predictions it is not in all cases sufficient to use the approximate approach.

Conclusions
One of the main tasks of the LHC is the search for beyond the SM physics, in particular supersymmetry. At the high-energy run of the LHC coloured SUSY particles can be produced with masses up to the multi-TeV range. In order to find these particles and be able to measure their properties, reliable predictions for the production cross sections both at the inclusive and at the exclusive level are mandatory. In this paper we continue our effort in providing accurate theoretical predictions by presenting results for the squark-antisquark production of the first two generations at NLO SUSY-QCD without making any simplifying assumptions on the sparticle masses and by treating the different subchannels individually. As developed in our previous paper [52] we have performed the subtraction of possible on-shell intermediate gluinos in a gauge-invariant approach and compared to several methods proposed in the literature. While in squark pair production for the investigated scenarios the differences in the total rates turned out to be negligible and quite small for distributions, in squark-antisquark production, where the contributions of the qg initiated channels are more important, larger differences were found. They amount to about 4% for the inclusive NLO cross section in the investigated scenario. Even larger effects are found in the distributions, where the discrepancies between the investigated methods can be up to 30% in the p T distribution of the radiated parton. The invariant mass distribution of the squark-antisquark pair is not affected by the chosen method, however, and only reflects the discrepancy in the total cross section.
The K-factor for squark-antisquark production has been found to be sizeable and positive with K ≡ σ NLO /σ LO ≈ 1.4 and the scale uncertainty is strongly reduced by taking into account the NLO corrections. The comparison of the results for individual K-factors for the subchannels contributing to squark-antisquark production and the K-factor obtained after summing the cross sections differ significantly, so that the use of a global K-factor in general does not lead to accurate predictions. Combining the NLO production cross section with LO decays of the (anti)squark into the lightest neutralino and (anti)quark leads to discrepancies of about 10% between the exact result and the one assuming a common K-factor. The more the branching ratios of the squarks for the specific decay channel under consideration differ, the more important the consistent treatment of the individual corrections becomes.
In a next step we have combined our results for squark pair and for squark-antisquark production with the decays of the final state (anti)squarks into the lightest neutralino and (anti)quark at NLO SUSY-QCD at fully differential level. In this context we have discussed two methods for the combination of production and decay with NLO accuracy in the kinematics. One is based on a Taylor expansion respecting unitarity, but suffering from possibly negative contributions. The second approach, which does not expand the total decay width entering the branching ratios of the decays, avoids this problem, however violates unitarity. The results for these two approximations and for the case where no expansion in the strong coupling is performed at all, differ by at most 4% for the total cross sections. In the jet distributions the discrepancies between the two approximations can be up to 15%, whereas in the / E T distribution they are purely given by the discrepancy in the total cross sections. In view of these small deviations, in particular for the inclusive quantities, we have adopted the unitarity preserving approach in the remaining numerical analysis.
The influence of the NLO corrections on the distributions has been investigated for several observables. While in the / E T distribution the deviation of the differential K-factor and the total one is of O(5%), the K-factor for the p T distributions of the two hardest jets can vary in a range of ±40%, hence the assumption of a constant K-factor clearly is not valid here any more.
In order to obtain realistic predictions for exclusive observables we have matched the NLO cross sections with parton showers using the Powheg-Box framework. The implementation is publicly available and can be downloaded from [58].
The matched NLO results have been interfaced with the p T ordered shower of Pythia6 as well as the default shower and the Dipole shower of Herwig++. To allow for a consistent comparison of the three showers, in Pythia the starting scale for the radiation off the decay products had to be modified. The largest differences in the three shower predictions is found in the pseudorapidity distribution of the third hardest jet. Thus Herwig++ predicts more jets in the central region, which is particularly pronounced for squark-antisquark production in the investigated scenario. The comparison of the showered result with a pure NLO simulation shows small differences for more inclusive quantities. In more exclusive distributions, in particular Herwig++ shows large deviations from the pure NLO result, as e.g. in the predictions for the third hardest jet. To decide if this is an effect of the missing truncated shower or a relict of the way the phase space is populated would require further investigations and is beyond the scope of this work.
Finally, we performed a cut-based analysis of the total cross sections in two benchmark scenarios using realistic event selection cuts taken from an ATLAS analysis. Comparing our results with the approximate approach used by the experiments revealed small discrepancies for squark pair production, but up to 20% differences for squark-antisquark production. This effect could be traced back to assuming a common K-factor for the production cross sections of all subchannels instead of using the exact results. These examples show that the effects can be sizeable and precise theoretical predictions should take into account the full NLO calculation for the production processes, consistently combined with the squark decays at NLO.
The reliable exploitation and interpretation of the LHC data in the search for new physics requires accurate theoretical predictions for production and decay of SUSY particles including higher order corrections not only for inclusive quantities but also for distributions. Our results for the fully differential calculation of the SUSY-QCD corrections to squark pair and squarkantisquark production combined with their subsequent decay at NLO SUSY-QCD and matched with parton showers show that the independent treatment of the contributing subchannels is essential and that differential K-factors can not be assumed to be flat. The results presented here are the next step in our program of providing a fully differential description of SUSY particle production and decay at the LHC.