Electroweakino pair production at the LHC: NLO SUSY-QCD corrections and parton-shower effects

We present a set of NLO SUSY-QCD calculations for the pair production of neutralinos and charginos at the LHC, and their matching to parton-shower programs in the framework of the POWHEG-BOX program package. The code we have developed provides a SUSY Les Houches Accord interface for setting supersymmetric input parameters. Decays of the neutralinos and charginos and parton-shower effects can be simulated with PYTHIA. To illustrate the capabilities of our program, we present phenomenological results for a representative SUSY parameter point. We find that NLO-QCD corrections increase the production rates for neutralinos and charginos significantly. The impact of parton-shower effects on distributions of the weakinos is small, but non-negligible for jet distributions.

While total production cross sections for electroweakino pair-production processes can be obtained for a large variety of MSSM parameter points via the public computer program PROSPINO or the RESUMMINO code package that additionally provides transverse-momentum and invariant-mass distributions of the gauginos, currently no dedicated Monte Carlo program exists for the calculation of differential distributions within arbitrary experimental selection cuts at NLO SUSY-QCD accuracy. Moreover, the afore-mentioned public programs do not provide information on the gaugino decay chains, and cannot be interfaced easily with Monte-Carlo programs such as PYTHIA [17] for the simulation of parton-shower, underlying-event, and multi-parton interaction effects. In principle, multi-purpose programs like MadGraph5_aMC@NLO [18] do provide building blocks for the computation of arbitrary processes at NLO-QCD accuracy. In the context of SUSY processes with a complex resonance structure, though, human interaction is required for the subtraction of on-shell resonances that currently cannot be accounted for automatically in the MadGraph5_aMC@NLO framework.
With the current work, we want to close existing gaps. We have developed a versatile code package for various electroweakino pair-production processes that provides NLO SUSY-QCD corrections to cross sections and differential distributions within arbitrary experimental selection cuts, and can be interfaced to parton shower programs via the POWHEG matching procedure [19,20]. We are using the framework of the POWHEG-BOX [21], a public repository for the simulation of scattering processes at hadron colliders at NLO-QCD accuracy matched with parton shower programs. For technical aspects related to features of the MSSM we build on experience gained in the implementation of slepton [22,23] and squark pair-production processes [24,25] in the framework of the POWHEG-BOX.

JHEP07(2016)083
In the following section, we will briefly describe technical aspects of our calculation that are specific to the implementation of electroweakino pair-production processes in the context of the POWHEG-BOX. In section 3 we will provide representative numerical results with a particular emphasis on the impact that NLO and parton-shower effects have on observables. Focusing on a specific SUSY benchmark point we will demonstrate how the code package we developed can be used for the calculation of experimentally accessible distributions including a simulation of supersymmetric decay chains. Our conclusions are given in section 4.

Framework of the calculation
Our implementation of weakino pair-production processes in the framework of the POWHEG-BOX builds on experience gained for related supersymmetric reactions, in particular slepton-and squark-pair production processes [22][23][24][25]. Rather than going into general features required for the implementation of a new process in the POWHEG-BOX repository, we here will focus on aspects that are specific to weakino pair production.
At the leading order the production of a pair of weakinos proceeds via the tree-level diagrams presented in figure 1 (a). In all channels the s-channel topology comprises Drell-Yan production, qq → V * , followed by the splitting V * →χ iχj , whereχ i stands for either a neutralinoχ 0 i (i = 1 · · · 4), or a charginoχ ± i (i = 1, 2), depending on the process of interest. The vector boson V denotes a Z boson in the case of neutralino pair production, V = W ± for the production of a neutralino and a chargino, and V = γ/Z for the production of a pair of charginos. In addition, diagrams with a squark being exchanged in the t-or u-channel occur. In the case of the production of a chargino and a neutralino only either tor u-channel contributions arise, while for the other considered production processes both types of topologies contribute. We work in a scheme with five massless quark flavors in the initial state, i.e. q/q = u, d, s, c, b. Numerically small bottom mass effects are disregarded throughout. This allows us to treat the scalar partners of these left-and right-chiral fermions as mass eigenstates. The CKM matrix is taken to be diagonal.
The NLO-QCD and SUSY-QCD corrections comprise virtual corrections to the qqinduced processes as well as real corrections with an extra parton in the final state. Only the sum of both corrections is infrared (IR) finite. Representative Feynman diagrams for the virtual corrections are shown in figure 1 (b). They include vertex and box corrections with gluon, gluino, quark, or squark exchange, as well as self-energy corrections in the case of the t-and u-channel diagrams with squark exchange. We use FeynArts 3.9 [26] to generate the virtual diagrams and FormCalc 8.4 [27] to calculate the amplitudes using the MSSM-CT model file of ref. [28]. The scalar loop integrals [29] are numerically calculated with LoopTools 2.12 [27,30]. In order to cancel the ultraviolet (UV) divergences, a renormalization procedure has to be conducted. In practice, this requires the calculation of suitable counterterms. We use the on-shell scheme for the renormalization of the wave functions of the massless incoming quarks as well as for the squark masses. Other fundamental parameters such as the electroweak coupling constant do not require renormalization at NLO in QCD.  In order to regularize the UV divergent loop integrals, in principle there are two possibilities in a supersymmetric theory. The standard procedure of the POWHEG-BOX is the dimensional regularization scheme (DREG), where the entire calculation is done in D = 4−2ε dimensions. The same regularization procedure is used in most publicly available sets of parton distribution functions (PDFs) that are needed for the computation of cross sections at a hadron collider. However, this scheme breaks supersymmetry at the level of the gauge interactions by introducing a mismatch in the (D − 2) transverse degrees of freedom of the gauge bosons and the two degrees of freedom of the gauginos. Hence, while SUSY invariance requires that the quark-squark-weakino Yukawa couplingĝ and the associated SU(2) gauge coupling g be equal at all perturbative orders, this relation is violated in DREG. In order to remedy this deficiency, a finite SUSY restoring counterterm is added at next-to-leading order in the strong coupling α s [8,31,32], The expansion in α s is done consistently to retain only the O(α s ) term that is induced by this finite SUSY restoring counterterm in the amplitude squared. An alternative way of regularization is the dimensional reduction scheme (DRED) for the calculation of the finite part of the virtual corrections. In DRED fields remain defined in four dimensions, while loop momenta are defined in D dimensions. Since this approach respects supersymmetry, the SUSY-restoring counterterm of eq. (2.1) is not needed anymore. However, to comply with the intrinsic treatment of the IR singularities in the POWHEG-BOX, a finite shift of the virtual amplitudes passed to the Monte-Carlo program is necessary [21],

JHEP07(2016)083
where B is the Born amplitude for a specific partonic subprocess calculated in four dimensions. As a cross-check, we have employed both regularization methods in our calculation, and have in both schemes found the same value for the virtual amplitude V that enters the Monte-Carlo program.
In order to calculate the real emission contributions and provide the ingredients necessary for the construction of IR subtraction terms by the POWHEG-BOX, we make use of a build tool based on MadGraph 4 [33][34][35]. It can be used to generate the Born, the colorand spin-correlated Born and the real-emission amplitudes in a format that can be easily processed by the POWHEG-BOX. The IR divergences are canceled separately in the virtual and in the real parts by using the Frixione-Kunszt-Signer algorithm [36] that is implemented in the POWHEG-BOX. Representative real emission diagrams are displayed in figure 2.
While the calculation of the real-emission contributions for the qq -induced subprocesses of type qq →χ iχj g is straightforward, a subtlety arises in crossing-related partonic subprocesses with a quark in the final state. As noted in ref. [8], subprocesses of the type qg →χ iχj q include two types of contributions: first, we encounter one-parton emission diagrams being part of the genuine NLO-QCD corrections to weakino pair production (representative diagrams are displayed in the upper row of figure 2). Second, there occur contributions that can be interpreted as tree-level diagrams for the on-shell production process pp →q kχi , followed by the squark decayq k → q χ j , if the squark is sufficiently heavy for a (quasi) on-shell decay into the respective weakino, i.e. mq k > mχ j (representative diagrams for this class of contributions are displayed in the middle row of figure 2). This on-shell feature emerges not only in weakino pair production, but also in other supersymmetric production processes involving squarks or gluinos [24,25,37,38], or in the case of tW production in the framework of the SM [39]. While the genuine real-emission contributions of the first type clearly have to be taken into account in our NLO calculation, the resonant contributions that are to be considered part of the different production process pp →q kχi need to be removed consistently in order to avoid double-counting. To this end, we make use of a scheme that in a slightly different variant has first been applied in PROSPINO [37], and more recently has been adapted for the code structure of the POWHEG-BOX in refs. [24,39]. More specifically, we extend the procedure developed for the related case of squark pair production in the POWHEG-BOX [24,25] to the richer resonance structure of weakino pair-production processes, as more diagrams are involved in that case.
Each resonant diagram occurring in a subprocess of type qg →χ iχj q exhibits a propagator that diverges when the intermediate squark goes on shell. For instance, in the case of aq →χ j q decay, this implies (pχ j + p q ) 2 → m 2 q in terms of the momenta of the external particles. Such would-be divergencies can be tamed by assigning a finite width Γq to the respective propagator, as is done by default in our MadGraph-based real-emission amplitudes.
Having determined the resonance structure of the resonant diagrams that are considered as part of a squark-weakino rather than a weakino pair-production process, we  Figure 2. Representative Feynman diagrams for the partonic subprocesses qq →χ iχj g (lowest row), for diagrams including a squark resonance in the qg →χ iχj q channel (middle row), and for non-resonant diagrams in the qg →χ iχj q channel (upper row). In each case, a = 1, 2.
are now in a position to devise an on-shell (OS) subtraction procedure for isolating the genuine weakino pair-production process of interest. To this end, we split the full realemission amplitude into a purely resonant contribution, M res , and a regular remainder, M reg . Specifically, in the case of neutralino pair production, eight resonant diagrams M k res with k = 1, . . . , 8 contribute to M res , while for production processes involving charginos, only four resonant diagrams occur. In order to remove the resonant squark contributions, for each resonant diagram squared M k res 2 we introduce a counterterm of the form where p res = (pχ j +p q ), and the momenta entering M k res are to be remapped to the on-shell kinematics, cf. ref. [38]. The first step function in this equation ensures that the partonic center-of-mass energyŝ is sufficient to generate both an on-shell intermediate squark and an on-shell spectatorχ i . The second theta-function guarantees that the squark has a mass larger than the sum of the masses of the two particles into which it decays, so that it can become on-shell. Since we only consider massless quarks, m q = 0 GeV in our calculation. We also stress that the decay width Γq introduced in eq. (2.4) is to be viewed as a technical -6 -

JHEP07(2016)083
regulator in the on-shell subtraction procedure. It may, but not necessarily has to, be identified with the actual physical decay width for the resonant squark. Since after the on-shell subtraction results should not depend on the resonant contributions, final results must be independent of Γq.
After having identified the counterterms needed for the on-shell subtraction procedure, we are in a position to perform the phase-space integration separately for the regular and the on-shell subtracted resonant contributions, While the regular contributions are to be evaluated for the standard real-emission kinematics, the counter-term contributions have to be integrated over a remapped phase-space that is obtained from the original three-body phase space dΦ 3 via a Jacobian factor J k , where sq = p 2 res , and λ denotes the Kaellen-function, If the limits of the phase-space integration would not be adapted appropriately, an integration over the entire three-body phase-space would combine on-shell and off-shell contributions inconsistently. For further details on the rescaling of phase-space, see refs. [24,25]. For the actual evaluation of σ OS real in the POWHEG-BOX, we have devised a routine allowing for a mapping of phase space according to a specific resonance structure k. Following this procedure, we can in principle handle an arbitrary number of resonance structures. The routine we developed could thus be used for future POWHEG-BOX implementations of other processes requiring an on-shell resonance subtraction. We note that the on-shell subtraction procedure we are using has the advantage of numerical stability, but violates gauge invariance as Γq = 0. In order to overcome this drawback alternative methods have been explored in the literature [24], but were found to exhibit other disadvantages such as the requirement of artificial cuts, and also being quite involved in a MadGraph-based implementation of the real corrections. However, gauge invariance violating contributions in the approach we are using are numerically negligible, as demonstrated by the independence of our results on the technical parameter Γq discussed below.
In order to verify the validity of our implementation, we have performed a number of checks. First, we have tested that, after the subtraction of on-shell resonances, for collinear momentum configurations real-emission and IR subtraction terms approach each other. Second, we have found that the dependence of our predictions for weakino pair-production cross sections on the technical regulator Γq is negligible. Figure 3 illustrates the regulator dependence of the neutralino pair-production cross section for a SUSY benchmark point that features squarks heavy enough to on-shell decay into a neutralino and a quark. We consider the mSUGRA spectrum SPS 1a [40] with m 0 = 100 GeV, [41], we use it in order to illustrate the technical details of the regulator dependence as it easily provides a spectrum for which the squark masses induce resonances to be regulated. We do not use this benchmark point for phenomenological studies. In the range Γq/mq = 10 −5 to 10 −1 , where mq = 556 GeV is the average of the four squark masses of the first generation, the dependence of the cross section on the regulator is entirely negligible, thus confirming the stability of the applied on-shell subtraction procedure. Finally, we have computed total cross sections at LO and NLO accuracy in the setup of ref. [8] and found agreement with the published results.

Phenomenological results
A collection of electroweakino pair-production processes will be made publicly available in the framework of the POWHEG-BOX via the project website http://powhegbox.mib.infn.it/. Since the public code can be used for specific user applications, we refrain from presenting an extensive numerical analysis here, but only intend to highlight some representative phenomenological results.
For our numerical studies, we consider proton-proton collisions at the LHC with a center-of-mass energy of √ s = 14 TeV. For the parton distribution functions (PDFs) of the -8 -

JHEP07(2016)083
proton we use the PDF4LHC15 NLO set [42] as implemented in the LHAPDF library [43]. Since no LO set is provided by the PDF4LHC15 working group, we use the NLO set also for the computation of LO results. Unless explicitly specified otherwise, we choose fixed values for the renormalization and factorization scales, µ R and µ F , proportional to the sum of the masses of the weakinosχ A andχ B produced in the specific process under consideration, µ R = µ F = ξµ 0 with µ 0 = mχ A +mχ B . The scale parameter ξ is set to one by default. When combining fixed-order results with a parton-shower program, we use PYTHIA 6.4.25 [17]. QED radiation, underlying event, and hadronization effects are switched off throughout. Partons arising from the real-emission contributions of the NLO-QCD calculation or from the parton shower are recombined into jets according to the anti-k T algorithm [44] as implemented in the FASTJET package [45] with R = 0.4 and η jet < 4.5.
As electroweak input parameters we choose the Z boson mass, m Z = 91.1876 GeV, the electromagnetic coupling, α −1 (m Z ) = 127.934, and the Fermi constant, G F = 1.16638 × 10 −5 GeV −2 . The other SM and MSSM parameters required for our calculations are provided in the form of a file complying with the SUSY Les Houches Accord (SLHA) [46,47] that can be computed with an independent external spectrum calculator. We have used the SuSpect 2.43 program [48] for the calculation of the spectrum and the SDECAY program [49] for the decay widths and branching fractions to obtain such an SLHA file. Specifically, we consider a minimal supergravity (mSUGRA) benchmark point suggested in ref. [50] that is consistent with a Higgs mass of about 126 GeV as well as further collider and dark matter constraints. This benchmark point is characterized by the following SUSY input parameters: For this benchmark point, we first consider the neutralino pair-production process pp →χ 0 1χ 0 1 . We find a total cross section of σ LO = 4.780 ab at LO and of σ NLO = 5.595 ab at NLO. The NLO SUSY-QCD corrections thus enhance the production rate by more than 15%. In order to quantify the dependence of these results on the unphysical renormalization and factorization scales, we have varied µ R and µ F in the range 0.1µ 0 to 10µ 0 around our default choice µ R = µ F = µ 0 = 2 mχ0 1 , cf. figure 4. At LO, neutralino-pair production is a purely electroweak process and thus only depends on µ F via the parton distribution -9 -JHEP07(2016)083  functions of the scattering protons. The scale behavior of the LO results directly reflects the µ F dependence of the (anti-)quark distribution functions in the probed kinematic regime. At NLO, additionally µ R enters and, in contrast to µ F being effectively accounted for only at lowest order, dominates the scale uncertainty of σ NLO . However, in the typically considered range µ 0 /2 to 2µ 0 the NLO cross section changes by only about 3%, indicating that the perturbative expansion is rather stable, and the scale uncertainty is reduced compared to the LO predictions. In order to assess the impact of the higher-order corrections and parton shower effects on kinematic features of weakino pair production, we consider the representative chargino pair-production process pp →χ + 1χ − 1 . Numerical uncertainties are at the permille level and not shown in the plots that follow. Figure 5 illustrates the transverse-momentum and pseudo-rapidity distributions of theχ + 1 at fixed order, and after the matching of the NLO result to the parton shower (NLO+PS). Analogous results are obtained for the other chargino,χ − 1 . In figure 6 (left) we depict the invariant-mass distribution of the chargino pair. As expected from the above discussion of total cross sections for the related case of neutralino-pair production, we notice that the normalization of these distributions changes significantly when going from LO to NLO. On the other hand, their shapes are only slightly affected by the NLO corrections, as illustrated by the dynamical K-factors,  More pronounced effects of the parton shower emerge in jet observables, such as the transverse-momentum distribution of the hardest jet shown in figure 6 (right). For the reaction pp →χ + 1χ − 1 , at NLO, jets can only result from a hard parton of the real-emission contributions. After matching with a parton shower, additional jets can occur that will, however, be mostly soft or collinear. From the displayed figure it is apparent that while in the fixed-order calculation the transverse-momentum distribution of the jet diverges towards small values of p jet T , the Sudakov form factor of the NLO+PS calculation tames this would-be divergence. We note, however, that a precise description of jet observables in weakino pair-production processes would require considering the related reactions with an associated jet being present at LO already. Only a full NLO calculation for theχ + production process would yield NLO-accurate predictions for jet distributions. In our calculation of pp →χ + 1χ − 1 , jet observables are described effectively only at LO accuracy and thus associated with significant theoretical uncertainties. Our results confirm the findings obtained in the context of a jet veto resummation formalism [51] for the related case of slepton pair production, which revealed that theoretical uncertainties at the lowest resummation order are large enough to weaken current exclusion limits relying on searches making use of jet vetoes.
In many SUSY scenarios, theχ 0 1 represents the LSP that, due to the requirement of R-parity conservation, does not decay. Being electrically neutral, such an LSP cannot be observed directly in the detector, but only via its imprint on the missing transverseenergy spectrum. Depending on the mass hierarchy of a SUSY parameter point, heavier neutralinos and charginos decay via chains into a combination of stable particles, such as partons, leptons, neutrinos, and the LSP. Particularly clean experimental signatures emerge from final states with charged leptons that are rare in the context of the Standard Model. A prime example is provided by the leptonic decay chain of the process pp →χ 0 2χ + 1 that gives rise to a three-lepton final state as depicted in figure 7. Having full access to supersymmetric decay chains in a Monte-Carlo simulation is thus of great phenomenological relevance. The codes we developed for weakino pair-production processes offer such an option by an interface to the SUSY decay feature of PYTHIA. We can thus provide predictions that are at the same time NLO accurate for the hard weakino pair-production process, include parton-shower emission effects, and give full access to the kinematic properties of the stable particles in specific decay chains using the narrow-width approximation.
To illustrate this feature, we focus on final states with three charged leptons plus missing transverse energy arising from theχ 0 2χ + 1 production process. For the setup of this simulation, we follow closely the strategy of the ATLAS analysis reported in ref. [5]. We only consider events with an electron, a positron, a muon, and a large amount of -12 -JHEP07(2016)083 missing transverse energy in the final state. Each charged lepton is required to exhibit nonvanishing transverse momentum, be located in the central-rapidity region and sufficiently well separated from each other in the rapidity-azimuthal angle plane, p T > 10 GeV , |η | < 2.5 , ∆R( , ) > 0.05 . (3.5) In addition, the missing transverse momentum is required to be large, This latter observable is computed from the negative sum of the final-state particles that are detected, i.e. the electron, positron, muon, and jets with a transverse momentum p j T ≥ 20 GeV, similar to what is done in the experimental analyses. As the sum of the transverse momenta of the final-state particles should add to zero, this is effectively similar to the sum over the non-detected particles, i.e. the LSP, the neutrinos emerging in the decay chain, and the softer jets. Figure 8 (left) shows the missing transverse momentum distribution obtained with our POWHEG+PYTHIA simulation after the cuts listed above are applied. Here and in the following, results are presented for the default NLO+PS setup obtained by matching the NLO result via the POWHEG formalism with PYTHIA, and for reference also for a LO sample matched with PYTHIA using the same parton-shower settings, referred to as LO+PS. The ratios help to quantify the impact of the NLO corrections in the presence of parton-shower effects on distributions of the decay products encountered in the considered reaction. We find that the general features of the NLO corrections are very similar for distributions of the decay particles as for the weakinos produced in the primary hard scattering process, pp →χ 0 2χ − 1 . In particular, the R ratio is flat over the entire range of missing transverse momentum, with a size of about 1.2 resembling the ratio of the integrated NLO and LO cross sections. Figure 8 (right) illustrates the invariant mass distribution of the e + e − system in the considered process. Apparently, the decay of theχ 0 2 into a lepton pair and theχ 0 1 LSP is dominated by e + e − pairs with an invariant mass close to the Z pole. Similarly, the decay of theχ + 1 into a µ + ν µ pair and an LSP features a lepton-neutrino pair dominated by the W resonance. Since the invariant mass of the µ + ν µ pair cannot be fully reconstructed because of the non-detectable neutrino we refrain from showing that distribution here. Similar to the case of missing transverse momentum, the R ratio turns out to be flat for the invariant mass distribution of the e + e − system, with slightly more statistical fluctuations far away from the resonance region at around M e + e − ∼ M Z than in the peak region.
The transverse-momentum distribution of the electron is depicted in figure 9. Because of the selection cuts of eq. (3.5) that we impose, no events with a transverse momentum smaller than 10 GeV occur. Over the entire plot range, the R ratio amounts to about 1.2, i.e. the NLO corrections are distributed rather uniformly for this distribution. The r.h.s. of figure 9 shows the azimuthal-angle separation ∆Φ e + µ + of the two positively charged leptons occurring in the e + e − µ + + E miss T final state. We note that the azimuthal-angle separation of the positron and the muon peaks at ±π. Also for this distribution, the impact of NLO corrections is flat over the entire range considered.

Conclusions
In this work, we have presented a new set of implementations for weakino pair-production processes in the framework of the POWHEG-BOX. The newly developed code allows for the calculation of the NLO SUSY-QCD corrections for the hard production process, and provides an interface to parton-shower programs such as PYTHIA via the POWHEG method. The program can process SLHA files obtained with an external spectrum calculator for the computation of a specific SUSY parameter point in the context of the MSSM. If desired, decay chains of the weakinos can be simulated with a dedicated option in PYTHIA.
We have described the technical aspects of the implementation specific to weakino pair-production processes. To illustrate the capabilities of the code package we developed, -14 -

JHEP07(2016)083
we have discussed phenomenological features of a few selected weakino pair-production processes focusing on theoretical uncertainties and the impact of parton shower effects on experimentally accessible observables. We have found that, in accordance with previous results reported in the literature, generally NLO corrections have a significant impact on production rates. Scale uncertainties of the NLO results are moderate, however. Partonshower effects are very small for weakino distributions, but can be significant for jet observables that are effectively described only at lowest order in perturbation theory. Getting better control on these would require a full NLO calculation for weakino pair-production in association with a jet. Finally, we have illustrated the capability of the code to account for the kinematic distributions of observables related to SUSY decay chains for a specific mSUGRA benchmark point. Any user of the code is free, however, to consider a SUSY spectrum of her own choice and obtain NLO+PS results for any set of observables within arbitrary selection cuts.