Next-to-leading-order QCD corrections and parton-shower effects for weakino+squark production at the LHC

We present a calculation of the next-to-leading order QCD corrections to weakino+squark production processes at hadron colliders and their implementation in the framework of the POWHEG-BOX, a tool for the matching of fixed-order perturbative calculations with parton-shower programs. Particular care is taken in the subtraction of on-shell resonances in the real-emission corrections that have to be assigned to production processes of a different type. In order to illustrate the capabilities of our code, representative results are shown for selected SUSY parameter points in the pMSSM11. The perturbative stability of the calculation is assessed for the $ pp \rightarrow \tilde{\chi}^0_1 \tilde{d}_L $ process. For the squark+chargino production process $ pp \rightarrow \tilde{\chi}_1^- \tilde{u}_L $ distributions of the chargino's decay products are provided with the help of the decay feature of PYTHIA 8.


Introduction
After the discovery of the Higgs boson [1,2], with properties compatible with the Standard Model (SM) expectation [3], identifying signatures of physics beyond the Standard Model (BSM) remains a prime target of the CERN Large Hadron Collider (LHC). A strong motivation to search for new particles with the properties of Dark Matter (DM) candidates is provided by astrophysical observations that can best explained in the context of scenarios with massive, neutral particles that only weakly interact with constituents of the SM. A theoretically appealing class of models to account for such particles is provided by supersymmetric theories that predict a number of new particles differing from their SM counterparts in their spin, and receive larger masses through the soft breaking of supersymmetry (SUSY). The lightest stable SUSY particle (LSP) could constitute a viable DM candidate particle, with all the right properties to account for astrophysical observations. Searching for such particles in the very different environment of a hadron collider provides complementary means to establish the existence of these as of yet hypothetical entities.
SUSY particles can be produced in various modes at hadron colliders. The largest cross sections are expected for the production of color-charged particles via the strong interaction, e.g. QCD-induced pair-production of squarks or gluinos [4][5][6][7][8][9][10][11][12][13][14]. However, via their decays into SM particles and lighter, weakly charged SUSY particles squarks and gluinos give rise to rather complicated final states. Alternatively, one could consider the direct production of gaugino pairs [15], a gluino in association with a gaugino [16] or a squark in association with a neutralino or chargino (here generically referred to as weakinos) [17]. In particular in SUSY scenarios where the lightest neutralino, theχ 0 1 , is the LSP such production modes could result in a pronounced signature with a hard jet (due to the strong decay products of the squarks) and missing transverse energy stemming from a neutralino [17,18]. in several compressed scenarios [19], where the masses of SUSY particles are very close to each other. A similar search by the CMS collaboration resulted in limits on various simplified models of DM and models with large spatial extra dimensions [20]. The same search strategies could also be applied to interpret data in terms of SUSY scenarios where the lightest neutralino provides a DM candidate and is produced in association with a squark that ultimately gives rise to a hard jet.
In order to best exploit existing LHC data, accurate predictions for the considered production processes are essential. Ideally, they are provided in the form of Monte-Carlo programs that can easily accommodate experimental selection cuts and be combined with external tools for the simulation of decays, parton-shower, and detector effects. With the work presented in this article, we aim at providing such a program for weakinosquark production processes in the framework of the Minimal Supersymmetric Standard Model (MSSM). The fixed-order perturbative calculation at the heart of our endeavor in many aspects parallels the work of Ref. [17]. However, while that reference focused on providing next-to-leading order (NLO) QCD and SUSY QCD corrections to a wealth of weakino-squark production processes and an assessment of multi-jet merging effects at tree level we provide a fully differential Monte-Carlo program matching the NLO-(SUSY)-QCD calculation with parton showers using the POWHEG method [21,22] as implemented in the framework of the POWHEG-BOX [23]. Special care is taken in the subtraction of on-shell resonances that might spoil the convergence of the perturbative expansion if not treated properly. To that end, procedures developed and refined in Refs. [24][25][26] for the proper treatment of SUSY processes involving on-shell resonances in the context of the POWHEG-BOX are adapted. Details of the SUSY parameter spectrum can easily be set via an input file following the SUSY Les Houches Accord (SLHA) [27].
This article is organized as follows: We first describe details of our fixed-order calculation and, in particular, the treatment of on-shell resonances in the framework of the POWHEG-BOX. In Sec. 3 we illustrate the capabilities of our program by applying it to a representative phenomenological study. We select a particular point in SUSY parameter space and consider weakino-squark production at the LHC resulting in signatures with a hard jet and a large amount of missing transverse momentum before concluding in Sec. 4. We would like to stress that users are free to perform studies of their own design for arbitrary points in the MSSM parameter space with the public version of our code, available from the POWHEG-BOX site http://powhegbox.mib.infn.it/.

Details of the calculation
Weakino-squark production at hadron colliders at leading order proceeds via two topologies, illustrated in Fig. 1: s-channel processes mediated by a quark, and t-channel processes mediated by a squark. In both cases, the initial state consists of a gluon and q g qχ ĩ q a qq aq a gχ i Figure 1: Representative tree-level Feynman diagrams for weakino-squark production through an s-channel (left) and a t-channel(right) process.
a quark. The weakino in the final state can be either a neutralino or a chargino. A neutralino-squark final state is characterized by a squark that has the same flavor as the quark in the initial state, e.g. gd →χ 0 id a , where the index a = 1, 2 represents the chirality index of the squark. The index i = 1, ..., 4 indicates the kind of neutralino produced. A chargino-squark final state will instead exhibit a squark of a different flavor, e.g. gd →χ − iũ a . In this case, i = 1, 2, as there are only two kind of charginos. We consider five massless quarks in the initial state, which allows us to treat the left-and right-handed squarks as mass eigenstates. We consider the CKM matrix to be diagonal.
For our calculation we are considering weakino-squark final states for all types of squarks and antisquarks of the first two generations and all types of weakinos, leading to a large number of independent channels. For neutralino-squark production, the four flavors of squarks and antisquarks can be produced both left-and right-handed, leading to 16 production channels for each neutralino and therefore to 64 channels in total, considering the four types of neutralinos. For chargino-squark production, only lefthanded squarks and anti-squarks are produced at the considered order of perturbation theory, as the quark-squark-chargino vertex, necessary to have a two-to-two process as seen in Fig. 1, only exists for left-handed squarks. Thus, considering the four flavors of squarks and the two types of charginos and their respective antiparticles, there are 32 channels for chargino-squark production.
The LO amplitudes for all channels have been generated using a tool within the POWHEG-BOX suite based on MadGraph 4 [28][29][30]. This tool also provides spin-and color-correlated amplitudes, that are necessary in order to use the automated version of the FKS subtraction algorithm [31] implemented in the POWHEG-BOX. To cross-check our results, we verified that these amplitudes were equivalent to those generated using FeynArts 3.9 [32] and FormCalc 9.4 [33] with the MSSM-CT model file from Ref. [34].
As a completely alternative approach we used an old (unpublished) calculation per-  of the t-channel LO diagram. Representative diagrams are shown in Fig. 2. Compared to the LO diagrams, the only new particle appearing in these diagrams is the gluino.
The virtual diagrams have been generated using the same tools used to generate the Born ones, FeynArts 3.9 and FormCalc 9.4, using the same MSSM-CT model file.
Since these diagrams are UV-divergent, a renormalization procedure has to be carried out. The divergences are first regularized using dimensional regularization and then removed using suitable counterterms, which can be generated by FeynArts. ance. This equality has to be restored introducing finite counterterms [15,36,37]. In our case, the relevant Yukawa coupling is the weakino-quark-squark one, for which the counterterm reads:ĝ The virtual diagrams also contain IR divergences, which are ultimately canceled by corresponding divergences in the real-emission corrections via the FKS subtraction method.
The real emission contributions have been generated using the aforementioned tool based on MadGraph 4. Differently from the LO case with quark-gluon induced channels only, the real corrections include q q and g g initial states, resulting in the following channels: Representative diagrams are shown in Fig. 4. Considering the various possibilities, the process that we are considering. They should instead be seen as the on-shell production of a different final state followed by a decay into two different particle, e.g. in the case of an on-shell resonant squark, as a di-squark production process followed by the decay of one of the squarks into a weakino and a quark. These contributions are therefore already taken into account in the respective production processes. Considering them as a part of our real corrections would lead to them being counted twice. They have therefore to be removed from our real-emission contributions to obtain a well-defined factorization of production and decay processes in the narrow-width approximation, which will then be again perturbatively well-behaved.
On-shell resonances rarely appear in SM processes (e.g. in W t production [38]), but are a relatively common feature of SUSY calculations. A similar subtraction procedure is necessary, for instance, in weakino-pair production [25,26] and in squark-pair production [24,39,40]. For the subtraction of on-shell resonances in our calculation we will resort to the procedure discussed in these references 1 .
The real-emission contributions to a process containing on-shell resonances can be divided in non-resonant (labelled as nr) and resonant (labelled as res) parts in the following way: where also an interference term between the non-resonant and resonant terms appears.
This separation is performed at the diagram level, meaning that only the diagrams that We remove the on-shell contributions from the resonant diagram, performing a pointwise subtraction of a counterterm. The first step is the regularization of the singularities present in the propagator of the on-shell particles, achieved by inserting a technical where m ij is the mass of the potentially on-shell particle ij and s ij = (p i + p j ) 2 , with p i and p j being the momenta of the final state particles that are daughter particles of the ij particle. The regulator Γ reg is not necessarily the physical width of the particle, but a technical parameter. The total cross section, after the removal of on-shell contributions, should not be dependent on its value, as the off-shell contributions are not. This also ensures that, while the applied method in general breaks gauge-invariance, this gaugeinvariance breaking effects are removed by going to the narrow-width approximation. This is achieved by approaching the plateau numerically, where the hadronic results do not depend on the regulator Γ reg anymore.
After the regularization, the removal of the on-shell contributions is performed by subtracting, locally from each resonant diagram, a counterterm that reproduces the behavior of the on-shell resonances. For a single on-shell resonance in a 2 → 3 process the general shape of this counterterm is: where again the particles originating from the potentially on-shell particle are labelled i and j, and the index k denotes the remaining particle, which is often referred to as the spectator particle. The first theta-function represents the condition that the centerof-mass energy squared s is high enough to produce the intermediate particle ij The Breit-Wigner factor BW is used to suppress the counterterm in correspondence of off-shell regions, to avoid the subtraction of off-shell contributions. It is defined as the ratio between the matrix element squared and the matrix element squared itself taken in the on-shell limit, i.e. s ij → m 2 ij and is equal to: In order for a BW factor to reproduce as closely as possible the behavior of the resonance, the regulator Γ reg is usually chosen such that it is as small as possible without causing numerical instabilities in the integration. In the limit Γ reg → 0, the BW factor approaches a Dirac delta, leading to fewer off-shell contributions being included in the counterterm.
Before integrating the counterterm over the phase space, however, another consideration is necessary: the matrix elements contained in the counterterm have been evaluated in a different phase space that meets the on-shell conditions, which is different from the phase space in which all the other real matrix elements have been calculated, i.e. the general three-particle phase space dΦ 3 . Therefore, the integration has to be performed over a separate phase space or, alternatively, a corrective factor reflecting the phasespace transformation, also called Jacobian factor, can be introduced before integrating the real contributions and the counterterms over the same phase space. If dΦ 3 is the onshell phase space element, it can be expressed in terms of the the regular three-particle phase space as: The Jacobian factor can be derived by imposing the on-shell condition to the integration over the dΦ 3 and reads: where the Källén λ function is defined as λ(x, y, z) = The cross section for the full real emission contributions can now be calculated as where σ OS is defined as There are two possibilities to perform the integration necessary to calculate these cross section. The simpler way is to actually not perform the integration in Eq. (9) before summing it into Eq. (8), but to perform only one collective integration for the whole real cross section. We will refer to this method as DSUBI. The second possibility is to perform separately the integrations, i.e. summing the contributions from Eq. (9) into Eq. (8) after they have been integrated. We will refer to this method as DSUBII. Requiring two separate integrations, the DSUBII method is obviously more time-consuming. However, it is the default option in our implementation, as for some scenarios, the DSUBI method was leading to numerical instabilities.
To better illustrate the on-shell subtraction scheme we employ let us specifically discuss the real-emission corrections to weakino-squark production. As already mentioned, squarks of every flavor and gluinos can become resonant in this class of processes, leading to 11 independent channels for resonances. In principle, a different regulator could be used for every channel. In our implementation, a distinction is made only between the regulator for squark resonances Γq and the one for gluino resonances Γg.
To check the stability of our on-shell subtraction method, we investigated the impact of the variation of the regulator on our results. We noticed that the DSUBI method was sufficient to give stable results when only gluino resonances were present but that when also squark resonances were possible, using the DSUBII method was necessary, because of the more complicated resonant structure. We therefore recommend the use of the DSUBII method. We also observed that the total cross section does not depend on the value of the regulator, as long as this value is chosen to be sufficiently small so that we have numerically reached the narrow-width limit. For both the gluino and the squark resonances, the ratio between the regulator and the mass of the resonant particle should not be larger than 10 −4 . Smaller values can be used, but they lead to larger numerical errors.
Our findings are shown in Fig. 6. We performed the calculation of the NLO cross section forχ 0 1d L using an artificial SUSY spectrum in which the mass values lead to the simultaneous appearance of gluino and squark resonances. The mass of the neutralino has been set to 314 GeV, the mass of the gluino has been set to 4.00 TeV and the squarks have been assumed to be mass degenerate with a mass of 2.31 TeV. The value of the two regulators has been changed simultaneously to explore the range 10 −8 ≤ Γ i m i ≤ 10 −1 , where i =q,g. In the range 10 −8 ≤ Γ i m i ≤ 10 −4 , the cross section appears to be fundamentally independent of the regulator, thus proving that our on-shell subtraction scheme is well-defined and that gauge-violating effects are numerically negligible once the the narrow-width approximation for the intermediate on-shell states has been reached.
For larger values of the regulator, the cross section is no longer constant.
As an additional check of our implemention, we reproduced the results presented in Ref. [17]  squarks and antisquarks for the individual production processes.

Phenomenological analysis
In order to demonstrate the capabilities of our code, we present numerical results for two selected scenarios: First, we consider the on-shell production of a neutralino and a squark for a realistic SUSY parameter point in the phenomenological MSSM model with eleven parameters, the pMSSM11, suggested in Ref. [41], which we call scenario a. This Our POWHEG-BOX code can be used to produce event files in the format of the Les Houches Accord (LHA) [45] for the on-shell production of a squark and a weakino.
These event files can in turn be processed by a multi-purpose Monte-Carlo program like  Table 2: Cross sections at LO and NLO accuracy and K factors for different final states in scenario a. The quoted numbers include the sum of channels for weakino production in association with a squark and an anti-squark of the given flavor. All cross-section numbers are given in units of [ab]. The numerical uncertainties do not affect the digits reported.
PYTHIA [46,47] that provides a parton shower (PS) to obtain predictions at NLO+PS accuracy. PYTHIA furthermore provides the means for the simulation of tree-level decays of unstable SUSY particles. To illustrate that feature, we consider the squark+chargino production channel pp →χ − 1ũ L for scenario b and simulate the decays of the squarks, u L → uχ 0 1 , and the charginos,χ − 1 → e −ν eχ 0 1 with PYTHIA 8 [47], thus providing predictions for pp → e −ν eχ 0 1 uχ 0 1 in the narrow-width approximation for the squark and chargino decays. We note that QCD corrections do not affect the purely weak chargino decay. QCD corrections in principle relevant for the squark decay are not taken into account. When labeling the perturbative accuracy of our results, we will only refer to the production process, implicitly assuming that no QCD corrections are provided for the squark decays.
Throughout our analysis, the renormalization and factorization scales, µ R and µ F are set to with The scale variation parameters ξ R , ξ F are used to vary the scales around their central value. If not specified otherwise, they are set to ξ R = ξ F = 1. For the parton distribution functions (PDFs) of the protons we use the CT14LO and CT14NLO set [48] as provided by the LHAPDF library [49] Figure 7: Dependence of the total cross section for the process pp →d Lχ 0 1 on the renormalization (µ R ) and factorization(µ F ) scales. Three curves are shown for the LO and NLO cases, respectively. The solid lines correspond to µ R = µ F = ξµ 0 , the dotted lines to µ R = ξµ 0 with µ F = µ 0 and the dash-dotted lines to µ F = ξµ 0 with µ R = µ 0 for scenario a.
In Tab. 2, we report cross sections and K factors, defined as the ratio of NLO to LO cross sections, for various final states in scenario a. We quote results for the production of squarks and anti-squarks of the first generation in association with either the lightest neutralinoχ 0 1 or charginoχ ± 1 . From our results the relevance of the NLO corrections is apparent. An exceptionally large K-factor of 4.16 is observed for theχ 0 1d R channel, due to the suppression of the LO matrix elements and the large numbers of quark-and antiquarkinduced channels of the real corrections.
In Fig. 7 we present the dependence of the cross section for the process pp →d Lχ can be traced back, to a large amount, to the sizable number of quark-and antiquarkinduced channels of the real corrections. The overall variation of the cross section in the range 0.5µ 0 to 2µ 0 is approximately 19%.
For our phenomenological study of the representative squark+neutralino production channel pp →χ 0 1d L in scenario a we assume the neutralino gives rise to a signature with a large amount of missing transverse momentum. We do not consider squark decays for this study, as it is mainly intended to demonstrate the perturbative stability of our results and the applicability of our code to searches for DM and other types of new physics.
In this scenario with stable squarks, jets can only arise from real-emission contributions or parton-shower effects. We reconstruct jets from partons using the anti-k T jet algorithm [50] with an R parameter of 0.4. While we do study the properties of such jets to assess the impact of the parton-shower matching on the NLO-QCD results, we do not require the presence of any jets in the event selection, but impose cuts only on the missing transverse momentum of the produced system. The missing transverse momentum of an event, p miss T , is reconstructed from the negative of the vectorial sum of the transverse momenta of all objects assumed to be visible to the detector (i.e. jets with a transverse momentum larger than 20 GeV in the pseudorapidity range |η jet | < 4.9 and squarks). The absolute value of p miss T is sometimes referred to as "missing transverse energy", E miss T = | p miss T |. In Fig. 8 we show the transverse-momentum distribution of the hardest jet, before any cuts are imposed. In our setup, with a stable squark, such a jet can only result from the real-emission contributions or from the parton shower. This distribution is thus particularly sensitive to the NLO+PS matching. Indeed, the figure shows the typical Sudakov behavior expected for this distribution: Towards low values of p jet T , the fixedorder result becomes very large. The Sudakov factor supplied by the NLO+PS matching procedure dampens that increase. At higher transverse momenta the NLO+PS results are slightly larger than the fixed-order NLO results. Since our analysis is fully inclusive with respect to jets, distributions related to the directly produced squark and neutralino do not exhibit strong sensitivity to this Sudakov damping. Figure 9 illustrates the angular separation of the squark from the neutralino, ∆φ(χ 0 1 ,d L ), and the missing transverse energy at NLO and NLO+PS accuracy before any cuts are imposed. The ∆φ(χ 0 1 ,d L ) distribution shows that squark and neutralino tend to be produced with large angular separation, and this feature is not significantly altered by parton-shower effects. For the missing transverse-momentum distribution, in the bulk the NLO and NLO+PS results are very similar. Only towards very large values of | p miss T | where we expect soft radiation effects to become important, the fixed-order NLO result is slightly larger than the corresponding NLO+PS result, which is related to the behavior of the jet at low p jet T discussed above. Using a cut on the missing transverse energy thus can be considered to be a perturbatively safe choice. In the following, in order to consider an event in our analysis, we require it to be characterized by a large amount of missing transverse energy, Next, let us consider the squark+chargino production channel pp →χ − 1ũ L at NLO+PS accuracy, combined with tree-level decays of the squark,ũ L → uχ 0 1 , and the chargino, χ − 1 → e −ν eχ 0 1 , as provided by PYTHIA 8. For the considered parameter point, the quark resulting from the squark gives rise to a very hard jet, shown in Fig. 11. Further jets can be generated by real-emission corrections and the parton-shower. As demonstrated by Fig. 11 (r.h.s), such jets exhibit an entirely different shape.
In addition to the jet, the squark decay gives rise to a very hard neutralino. It contributes to the missing momentum of the final-state system together with the neutralino and the neutrino stemming from the chargino decay and jets that are too soft to be identified. The total missing transverse energy of the e −ν eχ 0 1 uχ 0 1 final state is depicted in Fig. 12. From this distribution it is clear that we can afford to impose a hard cut on E miss T to suppress background processes with typically much smaller amounts of missing transverse momentum, without significantly reducing the signal cross section. In the  following we therefore impose a cut of which reduces the cross section for e −ν eχ 0 1 uχ 0 1 production by less than one percent. Because of the mass pattern of mother and decay particles, the momentum balance of the chargino-decay system is quite uneven: The heavy neutralino acquires a large amount of transverse momentum, peaked at about 450 GeV, while the lepton exhibits a much softer distribution, see Fig. 13. Even softer leptons are expected in SUSY scenarios where the masses of the chargino and its decay neutralino are yet closer than for the considered pMSSM11 point. This feature has to be taken into account in searches aiming to use the lepton to tag a particular signal.

Conclusions
In this article we presented a calculation of the NLO-QCD corrections to the entire range of weakino+squark production processes, and their matching to parton-shower programs. We discussed in detail the subtraction of on-shell resonances that appear in the real-emission corrections, but de facto should be considered as part of a different production process. Two completely independent implementations of the fixed-order calculation were devised and compared to ensure the reliability of our work. One of these constitutes the central element of a new POWHEG-BOX implementation that will be made public at http://powhegbox.mib.infn.it/.
In order to illustrate the capabilities of this implementation, we presented results for two parameter points of the pMSSM11 for the representative pp →χ 0 1d L and pp →χ − 1ũ L processes, considering tree-level decays of the squark and chargino in the latter. We PYTHIA decays of the unstable SUSY particles of the squark+weakino production process can be simulated.
applications of the program we developed, and hope that henceforth the tool will find ample use in customized applications by the phenomenological and experimental highenergy communities.