Automated simulations beyond the Standard Model: supersymmetry

The MadGraph5_aMC@NLO framework aims to automate all types of leading- and next-to-leading-order-accurate simulations for any user-defined model that stems from a renormalisable Lagrangian. In this paper, we present all of the key ingredients of such models in the context of supersymmetric theories. In order to do so, we extend the FeynRules package by giving it the possibility of dealing with different renormalisation options that are relevant to supersymmetric models. We also show how to deal with the problem posed by the presence of narrow resonances, thus generalising the so-called on-shell subtraction approaches. We extensively compare our total rate results with those of both Prospino2 and Resummino, and present illustrative applications relevant to the 13 TeV LHC, both at the total-rate and differential levels. The computer programmes that we have used to obtain the predictions presented here are all publicly available.


Introduction
After more than fifty years since its proposal, the Standard Model (SM) has been proven to be an extremely successful theory of Nature: its predictions agree well with the vast majority of the data collected so far, sometimes at an astonishing level of precision. Despite its success, however, the SM leaves some deep questions unanswered, and suffers from various conceptual issues and limitations; thus, it is widely understood as a low-energy effective theory that is supposed to emerge from a suitable UV-complete theory. Among the candidates for the latter, weak-scale supersymmetry (SUSY) constitutes one of the best-motivated options from a theoretical viewpoint. Naturally extending the Poincaré algebra by linking the fermionic and bosonic degrees of freedom of the theory [1][2][3][4][5][6][7][8][9], SUSY and in particular its minimal incarnation, the so-called Minimal Supersymmetric Standard Model (MSSM) [10,11], addresses several of the shortcomings of the SM. For example, SUSY can tackle some aspects of the hierarchy problem inherent to the SM by stabilising all scalar masses relatively to quantum corrections [12], it can ensure the unification of the three SM gauge couplings at high energies [13][14][15], and many SUSY realisations include a candidate explaining the presence of dark matter in the Universe [16,17].
The existence of a superpartner with a mass equal to the one of its SM counterpart is experimentally excluded, so SUSY must be broken. For theoretical and phenomenological reasons, this breaking must be soft and is expected to shift the SUSY particle masses in the TeV regime. Thus, also owing to the solid theoretical motivations of supersymmetric theories, the quest for SUSY particles still plays a major role in the searches carried out by the LHC experiments, and features prominently in the strategies that are currently being laid out for future-collider projects. Both the ATLAS and CMS collaborations have begun to release full LHC run 2 results; the persistent absence of beyond-the-SM (BSM) signals implies that the space of parameters of SUSY theories is constrained in an increasingly severe manner. For example, squarks and gluinos are bounded to have masses well above 1-2 TeV [18][19][20][21], whilst the bounds on the electroweak superpartners now reach several hundreds of GeV [22,23]. However, such experimental limits have been extracted either in the framework of specific MSSM benchmark scenarios, or in the context of simplified models of new physics that are inspired by the MSSM, and that feature only a small number of new particles and new interactions with respect to those of the SM [24,25]. Therefore, these limits can be evaded both in various non-minimal supersymmetric theories and with specific, fine-tuned, parameter configurations of the MSSM.
Experimental analyses for SUSY searches are currently mostly based on Monte Carlo simulations of signals in which tree-level matrix elements of different partonic multiplicities are consistently combined and interfaced with parton-shower (PS) Monte Carlos; the absence of double counting between matrix elements and parton showers is guaranteed by the use of a merging prescription -CKKW [26,27], MLM [28], and CKKW-L [29,30] (see also refs. [31][32][33]). For any given jet multiplicity, such predictions are therefore leading-order (LO) plus leading-logarithm (LL) accurate. They are typically further improved by normalising them, at the fully inclusive level, to the best available (in the sense of perturbative information) total cross sections for the production of the relevant SUSY particles. Typically, these cross sections combine fixed-order predictions at the next-to-leading-order (NLO) accuracy in α S (which we shall call "NLO QCD" henceforth, understanding that quarks, gluons, squark, and gluinos run in the loops), with calculations that carry out the resummation of threshold logarithms at the next-to-leading-logarithmic (NLL) accuracy. In some cases, the results include some of the next-to-next-to-leading-order (NNLO) and next-to-next-to-leading-logarithmic (NNLL) contributions as well.
As the LHC moves steadily towards a precision-physics phase, BSM search strategies, in view of their null results thus far, must evolve too. Among other things, finding elusive signals requires improving the control over systematics. At the theoretical level, this demands a better use of higher-order results than the simple rescaling of tree-level merged simulations, without losing the realistic description of multi-jet final states provided by parton showers. Fortunately, solutions originally devised for SM physics, which have been thoroughly and successfully tested since the start of the LHC operations, can be applied to BSM scenarios as well, pretty much as their tree-level counterparts. More specifically, NLO computations can be matched to parton showers (NLO+PS henceforth) by means of the MC@NLO [34] and Powheg [35] methods (or of any of their variants and less-used alternatives [36][37][38][39][40][41][42][43]) for any underlying SUSY process. Likewise, the extensions of tree-level merging techniques to NLO [44][45][46][47][48][49][50][51][52][53][54][55] work equally well in the SM and BSM theories. These facts have stimulated recent theoretical work (which we shall briefly review in section 2.4), that is characterised by being essentially on a process-by-process basis and in simplified models.
The main goal of this paper is to render systematically feasible, for any user-defined process and in the fully-fledged MSSM, fixed-order, matched, and merged simulations that include NLO QCD effects. In order to do this, we shall rely on automated techniques, which are by now extremely well established; in particular, we shall work in the MadGraph5 aMC@NLO framework [56] (MG5 aMC henceforth). We remind the reader that MG5 aMC is a meta-code, which uses hard-wired information on Quantum Field Theories (such as general rules of Feynman diagrammatics, the structures of matrix elements, phase spaces, and cross sections, and so forth), and external information equivalent to the Lagrangian of the theory one is interested in (these are a set of rules, called a UFO model [57], obtained automatically from the Lagrangian by means of codes such as FeynRules [58], possibly used jointly with NLOCT [59]), to construct on the fly a computer code specific to the production process that one is interested to simulate.
Thanks to the significant amount of development and validation activity performed on MG5 aMC in the past few years, most of the work for the current paper has gone into the construction of a UFO model for the MSSM (see section 2). A notable exception, at the level of the meta-code, is the following. In the context of theories with rich particle spectra, that include relatively narrow resonances, the computation of contributions beyond the LO might give results spoiling the "convergence" of the perturbative series. This happens when real-emission corrections include partonic sub-processes that feature those resonances in s-channels, whose integration over the phase space is either divergent (if the resonance propagator is not Dyson-summed), or grows like an inverse power of the resonance width (and is therefore numerically dominant). A familiar example in the SM is that of tW − production, for which a subset of real corrections features a tW −b final state, to which diagrams with an s-channelt quark contribute. We point out that the lack of convergence of the perturbative series in these cases is perfectly justified, since one is trying to compute higher-order corrections to the cross section of a process which is simply ill-defined. In the SM example just mentioned, tW − production does not exist as such: its definition requires the simultaneous presence of a particle (the top) and of one of the decay products of its antiparticle (the W − ). In spite of this, one can try and give an operative meaning to these ill-defined processes: they are conceptually useful, since they correspond to an intuitive physical picture which is easy to understand, and they can constitute, within reasonable approximations, valuable perturbative tools. In this sense, BSM theories in general, and SUSY in particular, provide one with many possibilities to test different strategies. In this paper we shall discuss the definitions of such strategies, that we collectively call Simplified Treatments of Resonances (STR henceforth), and their implementations in MG5 aMC in a way suited to NLO+PS simulations. Thus, STR encompass the procedures known as Diagram Removal (DR) and Diagram Subtraction (DS), and generalise the so-called On-Shell Subtractions (OSS) schemes.
This paper is organised as follows. In section 2 we briefly review the status of simulations in SUSY theories, and the basic characteristics of MG5 aMC and of the models it uses for BSM physics. We also introduce the theoretical framework of the MSSM, and its various renormalisation options. We validate our implementation against Prospino2 and Resummino in section 3. The total cross sections of several SUSY benchmark processes in simplified scenarios are presented in section 4. The general algorithms for the treatments of resonances in perturbative computations, and their implementations in MG5 aMC, are discussed in section 5. A case study of jets plus missing energy at NLO+PS accuracy at the 13 TeV LHC is considered in section 6, in the context of a non-simplified benchmark point. We draw our conclusions in section 7. Some details about one-and two-point loop integrals, the MoGRe package (the acronym standing for More General Renormalisation in FeynRules), and the settings relevant to the decoupling mass limit in MG5 aMC can be found in the appendices A, B, and C, respectively. The MSSM, the simplest SUSY extension of the SM, results from the direct supersymmetrisation of the SM [10,11]. The gauge symmetry group of the theory is the SM one, and it relies on the three vector supermultiplets V B , V W and V G , where the SM left-handed (q L ) and right-handed (u R and d R ) quarks, as well as left-handed ( L ) and right-handed (e R ) leptons, are complemented by scalar counterparts, the left-handed ( q L ) and right-handed ( u R and d R ) squarks and left-handed ( L ) and right-handed ( e R ) sleptons. We have moreover indicated, in the equation above, the representation of the various supermultiplets under the MSSM gauge group SU (3) c × SU (2) L × U (1) Y . The MSSM Higgs sector is constituted of two chiral supermultiplets H D and H U , which allows for the cancellation of chiral anomalies and to generate masses for all up-type and down-type particles.
All kinetic and gauge-interaction Lagrangian terms are fixed by gauge and SUSY invariance, and can be casted very compactly within the superspace formalism [6][7][8], Gα θ·θ + h.c. , (2.4) where the notation [ . ] X indicates that after expanding the superfield inside the bracket in terms of the Grassmanian variables θ andθ, only the X-component is retained. The first line of the Lagrangian refers to the chiral sector of the theory and includes a sum upon the superfields associated with all the previously-introduced chiral supermultiplets. The vector superfields appearing in the exponents are considered contracted with the relevant representation matrices, and g y , g w and g s denote the hypercharge, weak and strong coupling constants. The last line of this Lagrangian describes the gauge sector, and involves squares of the superfield strength tensors associated with the three gauge subgroups, the summed spin indices α being indicated explicitly.
Assuming R-parity conservation to avoid the presence of baryon-and lepton-number violating interactions, which will challenge the experimental observations [60], the superpotential interactions include Yukawa couplings, generating in particular the SM quark and lepton masses, as well as the off-diagonal Higgs mass-mixing µ term. The corresponding Lagrangian in superspace is written as where the superpotential W MSSM (Φ) is given, in the flavour space and with all flavour indices understood for clarity, by (2. 6) In this expression, the matrices y u , y d and y e are the usual 3 × 3 Yukawa matrices. The expansion of the above Lagrangian in terms of the component fields relies on standard techniques detailed, e.g., in refs. [61,62], and allows for the extraction of the SUSY-conserving part of the MSSM Lagrangian. The results read where all fermionic fields are two-component left-handed (ψ and V ) or right-handed (ψ and V ) spinors, the dot products are invariant products in spin space, σ µ = (1, σ i ) consists of one of the possible four-vectors built upon the Pauli matrices and we used the usual gauge field strength tensor V µν and covariant derivatives D µ taken in the appropriate representation of the gauge group. The summation over k refers to the three MSSM gauge groups (the generic notation T k and g k being used for the gauge-group representation matrices and coupling constants), whereas we use the indices i and j and the generic (φ, ψ) notation for the sum over the scalar and fermionic component of the model supermutiplets. In addition, the gaugino-scalar-fermion interactions are included in the last term of the second line, the third line is constituted of the Yukawa interactions deduced from the superpotential, and the scalar potential is shown in the last line of the above Lagrangian.

Soft supersymmetry breaking and particle mixings
As for any phenomenologically realistic SUSY theory, the MSSM exhibits soft SUSY breaking to introduce a mass splitting between a SM particle and its superpartner. We supplement the Lagrangian of eq. (2.7) by all possible soft terms breaking SUSY explicitly [63], The soft SUSY-breaking Lagrangian first includes mass terms for the gauginos, the mass parameters of the U (1) Y , SU (2) L and SU (3) c gauginos being denoted by M 1 , M 2 and M 3 . The following seven terms consist of mass terms for all scalar fields, the parameters mQ, mL, mŨ, mD and mẼ being 3×3 Hermitian matrices in the flavour space and m Hu and m H d are the two Higgs mass parameters. In addition, bilinear and trilinear soft multiscalar interactions can be deduced from the form of the superpotential, the corresponding coupling strengths being organised in the three 3 × 3 T u , T d and T e matrices (in the flavour space) and the complex number b related to the soft SUSY-breaking Higgs mixing term. At the minimum of the scalar potential, the neutral components of both Higgs doublets H 0 u and H 0 d get non-vanishing vacuum expectation values and the electroweak symmetry is broken. As a result, the electroweak vector bosons B µ and W 3 µ mix into the massless photon A µ and massive Z µ states, where we have introduced the electroweak mixing angle θ w . As in the SM, the charged weak boson physical states are defined by diagonalising the third generator of SU (2) in the adjoint representation, which gives rise to The eight degrees of freedom included within the two Higgs doublets give rise to three Goldstone bosons G ± and G 0 that become the longitudinal modes of the weak gauge bosons and five physical Higgs bosons h 0 , H 0 , A 0 and H ± defined by diagonalising the Higgs sector, where we have introduced the Higgs mixing angles α and β, the tangent of the latter being given by the ratio of the vacuum expectation values of the neutral Higgs gauge eigenstates, tan β = v u /v d . Mixings in the fermionic electroweak sector are also induced by the breaking of the electroweak symmetry. In the neutral sector, the gaugino and Higgsino gauge eigenstates mix into the physical neutralino eigenstates χ 0 i (with i = 1, 2, 3, 4), whereas in the charged sector, they mix into the two chargino states χ ± i (with i = 1, 2), (2.12) In those relations, the mixing matrices N , U and V are unitary and allow to diagonalise the neutral and charged electroweakino mass matrices. As in the SM, the diagonalisation of the quark sector requires four unitary matrices V u , V d , U u and U d and the one of the lepton sector relies on three unitary rotation matrices V e , V ν and U e as we have omitted the right-handed neutrino (super)fields in the model definition. This leads to the following field redefinitions in the flavour space, and we follow the traditional approach of casting these rotations through a redefinition of the left-handed down-type quark field only via the CKM matrix V CKM , (2.14) In agreement with SUSY, these field redefinitions are promoted to the superfield level so that one must consider an extra rotation acting on the basis of left-handed down-type squarks similar to the one of eq. (2.14),d The two field redefinitions of eq. (2.14) and eq. (2.15) define the so-called super-CKM basis [64] in which the resulting 6 × 6 squark mass matrices are non-diagonal. Following the SUSY Les Houches Accord (SLHA) conventions [65], the superpotential and soft parameters are redefined according to where the hatted quantities refer to the new free parameters of the theory. The matricesŷ u ,ŷ d and y e are diagonal and real 3 × 3 matrices in the flavour space, whereas all the other matrices are in principle possibly flavour-violating and CP -violating. We however assume a constrained realisation of the MSSM in which organising principles of the soft terms forbid any source of flavour and CP violation on top of those inherent to the CKM matrix. In this case, all subsequent flavour-violating effects in the squark sector are small, and each squark flavour turns out to be aligned with the associated quark flavour. Following the SLHA conventions and the standard MSSM literature, theT f matrices are decomposed asT where the overall strength of the trilinear scalar interactions for a (s)fermion species f is embedded in the three 3 × 3 A f (diagonal and real) matrices in the flavour space. As a consequence and for not too extreme values of the coupling strengths A f , typical MSSM scenarios only exhibit a flavour-conserving mixing of the third generation sfermions, any other mixing being subdominant and negligible. Such a mixing is modelled through the stop (θt), sbottom (θb) and stau (θτ ) mixing angles, and the stop (t 1 ,t 2 ), sbottom (b 1 ,b 2 ) and stau (τ 1 ,τ 2 ) mass-ordered physical states are related to the corresponding gauge eigenstates as Although the tree-level form of all mixing matrices introduced so far can be easily calculated from the Lagrangians of eq. (2.7) and eq. (2.8), loop correction effects are important. One-loop and known two-loop contributions are hence in general included in all available MSSM spectrum generators [66][67][68][69].

Generalities
Focusing in the following on NLO calculations in α S , we rotate the Lagrangians of eqs. (2.7) and (2.8) to the mass basis and omit from the discussion any term that is irrelevant with respect to the strong interaction, In this expression, L (QCD) SM denotes the QCD part of the SM Lagrangian involving quarks and gluons, the sum overq k refers to a sum over all twelve squark mass-eigenstates (ũ L ,c L ,t 1 ,ũ R ,c R ,t 2 ,d L , s L ,b 1 ,d R ,s R ,b 2 ) of masses mq k and the sums over q refer to sums over all six quark flavours. In the former sums, the mixing matrices associated with the first and second generation squarks are taken as 2 × 2 identity matrices, so that the first and second generation mass-eigenstatesq 1 andq 2 are subsequently identified with the left-handed and right-handed squarksq L andq R , respectively. In addition, mg stands for the gluino mass, T for the fundamental representation matrices of SU (3), P L,R for the left-handed and right-handed chirality projectors, g s = √ 4πα S for the strong coupling constant and the covariant derivatives are restricted to their QCD component.
Ultraviolet divergences appearing at the one-loop level are absorbed into the counterterms generated by the renormalisation of the above Lagrangian. Following the usual procedure, all bare bosonic fields Φ and fermionic fields Ψ are replaced by their renormalised counterparts, with the exception of third generation squarks for which matrix renormalisation is in order as they mix, q 1 δZq ,21 δZq ,22 Although the structure of the gluino-squark-quark interactions (the last line of eq. (2.20)) could induce the mixing of any squark flavour at the one-loop level, those effects are proportional to the corresponding quark masses. Considering n lf = 4 flavours of massless quarks, the first two generations are kept non-mixing so that gauge and mass eigenstates are equivalent. In addition, the bare parameters of the MSSM Lagrangian, generically denoted by y (for couplings) and m (for masses), are renormalised as y → y + δy and m → m + δm . (2.23) In this work, we calculate the various renormalisation constants appearing in the renormalisation procedure of the Lagrangian of eq. (2.20) in the on-shell (OS) scheme where the input parameters are physical observables such as the physical particle masses. There is however no unique definition of such a scheme in SUSY by virtue of existing interrelations between various mass and coupling parameters, which will be addressed in section 2.2.4. Fermion self-energy corrections Σ(p) are decomposed in terms of independent Lorentz structures, from which the OS fermionic wave-function renormalisation constants δZ L,R f and mass renormalisation constant δm f can be deduced. Imposing that the renormalised mass is the pole of the propagator and that the residue of the propagator pole equals one, we get where the prime denotes a derivative with respect to p 2 . Gauge-boson self-energy corrections D µν (p) are reduced, in the case of the gluon (that is the only relevant gauge boson as long as only α S corrections are concerned), to their purely transverse component, The wave-function renormalisation constant δZ g is obtained after imposing OS renormalisation conditions, and reads Finally, as above-mentioned, the scalar quark sector of the theory relies on matrix renormalisation, any off-diagonal element being vanishing in the case of the first and second generation non-mixing squarks. We rewrite the scalar self-energies Π ij (p) as where the indices i and j are either 1 or 2. This allows for the derivation of the diagonal (δZ ii ) and non-diagonal (δZ ij ) wave-function renormalisation constants, where m 2 i indicates the squared mass of the i th eigenstate, as well as of the mass renormalisation constant, We do not address the complex-mass renormalisation scheme in this paper, so that we kept implicit that only the real-part of the self-energies is considered in the above expressions. The case of complex-mass renormalisation conditions is especially delicate in the case of NLO computations within SUSY theories because their mass spectrum is arbitrary to a large extend thus making it necessary to implement the most general analytic continuation of the two-point functions [70].

The Standard Model sector
Beginning with the SM sector, the wave-function renormalisation constants δZ L,R q of the massless quarks q = u, d, s, c and δZ L,R Q of the massive bottom and top quarks Q = b, t, as well as the one of the gluon δZ g are given by where the B 0 and B 1 functions and their derivatives are the real part of the usual two-point Passarino-Veltman integrals [71] collected in appendix A. Moreover, n c = 3 and C F = (n 2 c −1)/(2n c ) stand respectively for the number of colours and for the quadratic Casimir invariant connected with the fundamental representation of SU (3), and the sum uponq k refers to a sum over all squark states. In addition, the bottom and top mass OS renormalisation constants δm Q (with Q = b, t) are given by Qi , (2.32)

Gluino renormalisation
Gluino renormalisation in the OS scheme is standard, and the corresponding wave-function renormalisation constant δZg ≡ δZ L g = δZ R g (as the gluino is a Majorana fermion) and mass renormalisation constant δmg read (2.33)

On-shell squark renormalisation
The naive on-shell scheme Using the standard OS formulas as presented in section 2.2.1 for deriving the wave-function renormalisation constants of the first and second generation squarks δZq and third generation squarks δZQ, we obtain Similarly, the corresponding mass renormalisation constants read This scheme however breaks weak interaction gauge invariance, as the physical squark masses are not allowed to be taken all independent. Left-handed up-type and down-type squarks of a given generation are indeed connected by SU (2) L so that they can consequently not be renormalised independently. Such a scheme is however useful and valid for many phenomenological applications relying on simplified models inspired by the MSSM in which only a few particles and a subset of all MSSM Lagrangian terms are supplemented to the SM, as for instance in the work of refs. [73,74] or for the numerical results presented in the following sections of this paper. In the latter case, the relations between the physical squark masses are ignored as the relevant terms are not present in the simplified model Lagrangian, so that all fields can be renormalised independently. This approach however breaks down as soon as one considers an entire generation of squarks and wants to retain SU (2) L gauge invariance as embedded in the MSSM.
In the rest of this subsection, we additionally present two of the most popular SUSY OS schemes, that will not be considered in our numerical simulation but that could easily be implemented in our framework as will be shown in section 2.3. Whilst the differences between all the possible schemes are formally of higher order, the corresponding higher-order contributions could be potentially large in some parts of the parameter space (for instance, when tan β is large). Moreover, the different schemes necessitate different sets of input parameters, which becomes relevant for comparing their respective predictions.
The 'm b on-shell' scheme As above-mentioned, gauge invariance under weak interactions implies that the masses of the lefthanded up-type and down-type squarks of a given generation are connected to a unique bare soft mass parameterm 2 Q . The tree-level squared mass matrices M 2 qu and M 2 q d associated with the up-type and down-type squarks of a given generation are indeed given, in the (q L ,q R ) basis, by where m Z , c 2β and s w stand for the mass of the Z-boson, cos 2β and sin θ w . We have furthermore indicated by a subscript n the relevant generation index, and m qu and m q d are the masses of the corresponding up-type and down-type quarks q u and q d . While for the first and second generations the off-diagonal terms vanish and those two matrices are diagonal, they must be further diagonalised for third generation squarks with the help of the two rotation matrices St and Sb, Consequently, one of the four mass parameters associated with the first or the second generation of squarks is a dependent parameter and cannot be renormalised by imposing naive OS renormalisation conditions. Similarly, care must be taken with the stop/sbottom sector where we have six quark and squark masses (m b , m t , mt 1 , mt 2 , mb 1 and mb 2 ), two mixing angles (θt and θb) as well as two soft trilinear interaction strengths (A t ≡ (A u ) 33 and A b ≡ (A d ) 33 ). All these parameters are related and thus receive one-loop α S corrections in a connected manner [75][76][77][78].
In the so-called 'm b on-shell' scheme, the renormalisation of the up-type and down-type squark sectors is performed simultaneously [79][80][81]. We consider the masses of the left-handed down and strange squarks as well as the one of the heaviest bottom squark as dependent parameters, with c w ≡ cos θ w , so that the corresponding counterterms are given by We have explicitly introduced in those expressions the dependence on the mixing angles whose renormalisation constants δθt and δθb are given by eq. (2.36).
As a result, the renormalised masses of the left-handed down and strange squarks and of the heaviest bottom squarks are shifted with respect to their pole masses m 2 d L ,pole , m 2 s L ,pole and m 2 b2,pole [76], where δmd L ,pole , δms L ,pole and δmb 2,pole stand for the naive OS renormalisation constants of eq. (2.35) and δm 2 d L , δm 2 s L and δm 2 b2 are the 'm b on-shell' counterterms of eq. (2.40). The tree-level masses are moreover given by eq. (2.39). These UV-finite shifts must in particular be accounted for when an entire MSSM spectrum is used, as typical MSSM spectrum generators solely output pole squark masses.
By virtue of eq. (2.37) and eq. (2.38), the strengths of the soft trilinear squark-Higgs interactions A t and A b also receive one-loop corrections in α S through their connection with the corresponding squark mixing angles, (2.42) The corresponding counterterms are given by as both the µ parameter and tan β do not receive α S corrections at one loop.
The 'A b /θb on-shell' scheme As a consequence of eq. (2.43), two of the three counterterms δm q , δA q and δθq are independent. There are thus various options for fixing the renormalisation conditions, that all lead to slight differences in the predictions. In the 'm b OS' scheme, the two δA q renormalisation constants are derived from the other counterterms. This is however known to yield potentially-unacceptably large threshold corrections to the bottom-quark pole mass due to the δA b counterterm when tan β is substantial [80,[82][83][84]. Whilst a fully DR renormalisation of the bottom sector (δm b , δA b and δθb) would avoid the problem, this is also known not to make manifest the decoupling of heavy particles. We therefore present here another commonly-used scheme in which the A b parameter is renormalised in the OS scheme via a kinematic condition on the coupling of the pseudoscalar Higgs boson A 0 to ab 1b2 pair. This approach relies on the proportionality of the A 0b 1b2 coupling to the product of the bottom Yukawa coupling (or the bottom mass) and the bottom trilinear coupling, so that shifts in one quantity can always be reabsorbed in the other one. In practice, we calculate the one-loop corrections to the above-mentioned vertex with appropriately chosen external momenta and include suitable wave-function corrections to avoid any infrared divergence, In the above expression, the F function originates from the one-loop corrections to the A 0b 1b2 vertex, where only the finite pieces of the loop integrals are retained (i.e. all pieces independent of 1/¯ in the conventions of appendix A) and µ R stands for the renormalisation/regularisation scale. As a consequence, the bottom mass counterterm is now a dependent parameter,

Renormalisation of the strong coupling
Our calculations require that the running of the strong coupling constant α S originates solely from the contributions of the gluons and n lf flavours of light quarks. We therefore renormalise the strong coupling by subtracting, at zero-momentum transfer, all massive particle contributions and MS contributions of all massless particles from the gluon self-energy [85][86][87]. They are then absorbed in the renormalisation constant of the strong coupling δα S with n lf = 4, where the sum in the last term includes all twelve squark species. The UV-divergent part of the renormalisation constant has been written explicitly in terms of the quantity where γ E is the Euler-Mascheroni constant and is related to the number of space-time dimensions The above renormalisation procedure however leads to a violation of SUSY as it introduces a mismatch between the strong coupling g s and the Yukawa interactionĝ s of a gluino with a squark and a quark. While these two couplings are equal at tree-level, as shown by the last term of the second line of eq. (2.7), the equality is destroyed by the difference in the number of fermionic gluino degrees of freedom and bosonic gluon degrees of freedom. This artificial breaking of SUSY is compensated by finite counterterms restoring SUSY invariance.
As we impose that the definition of the strong coupling g s is the SM one due to the decoupling theorem, only the quark-squark-gluino vertices and quartic squark interactions have to be shifted [88]. The SUSY restoring counterterm Lagrangian L (SCT) MSSM is then given, in the gauge eigenbasis, by where adjoint colour indices have been included and a sum over (s)quark flavours is understood for clarity.

Technical details on the model implementation in FeynRules
In order to calculate SUSY particle-production (total and differential) rates at colliders and to simulate MSSM signals by matching fixed-order results at the NLO accuracy with parton showers, we rely on the MG5 aMC framework [56]. Our methodology is based on the joint usage of the Feyn-Rules [58], NLOCT [59] and FeynArts [89] packages to automatically produce a UFO model [57] that can be used by MG5 aMC. However, there are substantial differences with respect to the procedure that has been followed for stop pair production [74], in the SUSY QCD case [73] and for slepton production [90], as a consequence of the non-trivial renormalisation procedure for the mixing angle and the trilinear scalar couplings detailed in section 2.2.4. After having implemented the model described in section 2.1 and its tree-level Lagrangian in terms of superfields, we make use of the superspace module of FeynRules [91] to re-express the MSSM Lagrangian in terms of the model physical degrees of freedom and four-component fermions. The renormalisation is then performed with the MoGRe package, that is introduced in appendix B and that is necessary for a flexible definition of the renormalisation scheme.
We firstly impose that all external parameters insensitive to QCD corrections are kept unrenormalised. We hence enforce vanishing renormalisation constants for all electroweak inputs (the Fermi constant G F , the inverse of the electromagnetic coupling at the Z pole 1/α and the Z-boson mass m Z ), the parameters of the Higgs sector (tan β, the α angle and the µ parameter), the slepton trilinear couplings (A e , A µ and A τ ), the soft masses associated with the electroweak particles and the electroweakino mixing matrices (U , V and N ) that are external parameters in the SLHA conventions [65]. Moreover, the first and second generation squark trilinear couplings (A u , A d , A c and A s ) are irrelevant as multiplied by a vanishing quark mass and will thus not be renormalised. These constraints are imposed by using the MoGRe`DefineUnrenormalizedParameter function introduced in appendix B.2.
Secondly, the quark mass dependence of the (remaining) trilinear squark-Higgs couplings (T u andT d ) as well as the one of the fermion Yukawa couplings must be made explicit to guarantee the correct functioning of NLOCT. This is achieved by making use of the RemovingInternalCst method introduced in appendix B.6.1. The same method is finally also used to replace all occurrences of the g s renormalisation constant in terms of the α S one.
Next, we indicate to the code that fields that are insensitive to the strong interaction at the oneloop level (the electroweak gauge and Higgs bosons) do not need to be renormalised. This is achieved by making use of the MoGRe`DefineUnrenormalizedField method detailed in appendix B.2. Whilst other purely electroweak fields such as electroweakinos or (s)leptons have in principle to be analogously tagged as unrenormalisable objects, they do not appear in any QCD vertex so that they will be automatically discarded by the code. We finally impose that all field wave-function renormalisation constants are real (via the MoGRe`RealFieldRenormalization method presented in appendix B.2).
In practice, the MoGRe package is initialised as which indicates to the code that the stop and sbottom fields mix and that the renormalisation of the four-scalar interactions can be ignored. While strictly speaking, four-scalar interactions cannot be ignored, restrictions omitting them are useful phenomenologically as these vertices rarely appear at tree-level. The constraints above-mentioned are then implemented as follows, The exact details of the renormalisation scheme must then be specified, as shown in appendix B.4. Focusing on the naive OS scheme, the stop and sbottom mixing matrices are renormalised on the basis of eq. (2.36), which is implemented as We subsequently make use of NLOCT to generate the ultraviolet counterterms and R 2 Feynman rules necessary to obtain, from the UFO interface of FeynRules, an NLO UFO model for the MSSM. Since we induce an artificial breaking of supersymmetry by the mismatch of the two gluino and (D − 2) gluon degrees of freedom, one needs to add to the model Lagrangian a set of finite counterterms allowing to restore supersymmetry when dealing with one-loop calculations. Enforcing the definition of g s to be the SM one, quark-squark-gluino and four-scalar interactions are the only interactions that need to be shifted. Those shifts are given by the following counterterm Lagrangian [88], where L shift is written in the gauge eigenbasis. Those counterterms are appropriately included in the MSSM UFO in two steps [74]. We first evaluate the associated Feynman rules with FeynRules and then provide the resulting set of rules to the UFO interface by means of the UVLoopCounterterms option of the WriteUFO method. In contrast to the previous approaches, the resulting model can be used, within MG5 aMC, beyond the simplified model context.
Prospino2 [93,101,105] and MadGolem [142][143][144] give predictions that are NLO-QCD accurate, while NLL-fast [145], NNLL-fast [146] and Resummino [147] also resum threshold logarithms. All codes, with the exception of MadGolem and Resummino, assume mass-degenerate squark spectra. Theoretical predictions at the differential level, that include both NLO and PS effects, are more recent and, so far, tackled on a process-by-process basis. The production of pairs of squarks [73,74,148,149], gluinos [74], electroweakinos [150], and sleptons [90,151,152] have all been considered in the past few years. Some of these computations have been carried out with MG5 aMC. We recall here that MG5 aMC makes use of the FKS method [153,154] (automated in the module MadFKS [155,156]) for dealing with IR singularities. The computations of one-loop amplitudes are carried out by switching dynamically between two integral-reduction techniques, OPP [157] or Laurent-series expansion [158], and tensor-integral reduction [71,159,160]. These have been automated in the module MadLoop [56,161], which in turn exploits CutTools [162], Ninja [163,164], Iregi [165], or Collier [166], together with an in-house implementation of the OpenLoops optimisation [167]. Finally, in the case of matching with PS, the MC@NLO formalism [34] is employed, whereas NLO multi-jet mergings rely either on FxFx [50] or UNLOPS [53].
We point out that the original MG5 aMC paper [56] had the goal of including as many information on Quantum Field Theories as possible in the meta-code, so as to allow it to simulate both SM and BSM processes by using the inputs in the form of UFO models constructed by codes such as FeynRules or Sarah [168]. Recently, the program has been upgraded, and can for instance now handle mixed-coupling scenarios, in particular QCD+EW simultaneous corrections [70]. However, the most general BSM calculations beyond LO feature a number of non-trivial characteristics that are absent in the SM. While fermion-flow-violating interaction vertices and non-renormalisable operators (which were not available at the time of the first release [56]) can now be handled, colouredsextet particles and the renormalisation-group running of new couplings are not yet included in MG5 aMC. Thus, we stress again that the current work will be limited to considering NLO QCD corrections to SUSY theories, whereby only quarks, gluons, squark, and gluinos can run in the loops. Moreover, real-emission contributions only consider additional massless SM particles in the final state, given that massive particle contributions are finite (and can thus be computed independently) and often numerically subleading (see the analogous discussion of ref. [70] that addresses Heavy Boson Radiation (HBR) in the context of the computation of NLO electroweak corrections).
We conclude this section by listing the UFO models that can presently be used for BSM simulations. These include simplified models, in which the SM is extended by colour-triplet and octet scalar particles [73,169], both gluinos and squarks [74] or sleptons [90], as well as by vectorlike quarks [170,171], a heavy top-philic scalar [172] or a spin-2 particle [173]. In the latter spin-2 case, new physics have also been previously explored in a semi-automated framework (in the sense where the virtual matrix elements are provided externally) based on MG5 aMC [174][175][176]. Various BSM setups in which the Higgs sector differs from the SM one have been released, such as the two-Higgs-doublet model [177,178], the Georgi-Machacek model [179], the Higgs characterisation model [180][181][182][183][184], and the SM effective field theory including dimension-six operators [185]. Higherdimension operators either affecting the sector of the top quark [186][187][188][189][190][191], dijet production [192], or Z-boson production [193] can be added as well. Moreover, the model library also allows for NLO+PS calculations in BSM models involving TeV-scale neutrinos [194], a left-right symmetry [195], as well as extra neutral and charged gauge bosons [196]. Finally, dark matter simplified models in which the dark matter particle is produced in s-channels are also available [197][198][199][200][201][202][203][204].
As was stressed in section 1, part of the present paper is devoted to creating a UFO model of the MSSM, which is still missing in an unrestricted framework (i.e. when going beyond the simplified-model approach). The lifting of such a restriction has also to do with the treatment of resonant contributions, also addressed here through the STR procedures.

Validation
Fixed-order NLO-QCD predictions for the total rates of specific two-to-two processes in the MSSM are currently available from three different standalone tools, namely Prospino2 [93,101,105], Resummino [147], and MadGolem [142][143][144]. A partial comparison of the results obtained with MG5 aMC and MadGolem has already been performed for coloured-scalar production, and agreement at the level of the numerical errors has been found [73]. Furthermore, the analytic expressions of all the R 2 counterterms of the MSSM model have been cross-checked against the results of ref. [205]. In this section, we employ MG5 aMC to compute total rates for several processes and specific choices of the MSSM parameters, and compare our predictions against those obtained with Prospino2 (that covers the production of any pair of strongly-or electroweakly-interacting superpartners in the case of a degenerate squark mass spectrum) and Resummino (that supports arbitrary SUSY mass spectra for the production of two electroweak superpartners). In the rest of this paper, all fermions are unambiguously four-component Dirac and Majorana ones, so that we replace the notation Ψ X by X. In particular, quarks, gluino and electroweakinos will be denoted by q,g, and χ, respectively.

Setup of the comparison
Although we shall focus on superparticle-pair production here, it should be clear that, in keeping with a general automation philosophy, MG5 aMC is not restricted to simulating processes with two-body final states. Furthermore, MG5 aMC lifts two other key limitations of current NLO QCD codes: firstly, the inability to tackle QCD-mediated production processes (Resummino); and secondly, the inability to support without approximation arbitrary (non-degenerate) squark mass spectra (Prospino2). In order to highlight these differences, we have opted to present a comparison of the cross sections, at the 13 TeV LHC, relevant to the following processes (where antisquarks are denoted with a star): in both regimes of degenerate and non-degenerate squark masses. We point out that we have explicitly checked thatχ +  [206,207], as the focus of our work is on the computation of cross sections and observables.
The model parameters corresponding to the considered SUSY benchmark point are specified in all three codes via a similar SLHA file [65], the contents of which are summarised in table 1. The complicated nature of the Prospino2 inputs prompted our use of two different sets of parton distribution functions (PDFs) for LO and NLO predictions, for which the appropriate value of α S (m 2 Z ) (equal to 0.08991 and 0.08314, respectively) had to be hard-coded: by default, Prospino2 uses a fixed value of α S (m 2 Z ) independent of the PDF set. We have used the LO and NLO central sets of the CTEQ6 PDFs [208], and additionally turned off the running of α S . Moreover, whilst Prospino2 keeps the exact dependence on the masses of the produced sparticles both at the LO and the NLO, the masses that appear in all of the internal squark propagators are set equal to some averaged value when working at the NLO. The K-factor, defined as the ratio of such an "averaged" NLO computation over the corresponding LO one, is then used to multiply the exact LO result to get the final (and hence, approximate) NLO prediction. In order to assess the quality of such an approximation, we have scanned the cross sections obtained with all three codes in the two different mass setups defined in table 1.
In the case of a spectrum with degenerate SUSY masses, the mass of the produced SUSY particle (M prod ) is set equal to 1.5 TeV, while the common mass of all other SUSY particles (M others ) is scanned over, in the range [100,1400] GeV. This insures that all of the masses that appear in internal propagators are equal to each other (i.e. the internal squarks are degenerate), with the possible exception of the propagators that involve the produced particles. For instance,

Parameter value
Parameter value Parameter value Table 1: SM and SUSY parameters of the benchmark point used for the comparisons performed in section 3. Dimensionful quantities are given in GeV. M prod denotes the mass of the produced particles in the processes of eq. (3.1), while M others denote the masses of all other SUSY particles in a degenerate mass setup (this quantity is scanned over, hence its range value). In the non-degenerate mass setup, the left-handed up squark mass is set equal to a fixed value of 1.4 TeV.
t 1t 1 production has diagrams with internalt 1 propagators: the corresponding mass is then kept equal to 1.5 TeV. In the case of a spectrum featuring non-degenerate SUSY masses, we use exactly the same configuration as for the degenerate case, except that this time the left-handed up squark mass (mũ L ) is set equal 1.4 TeV. This particular choice for breaking the degeneracy pattern allows for an increase of the sensitivity of the inclusive cross section to the mass splitting mũ L -M others . We stress that it is crucial to set the mass of the produced particles M prod to a value larger than all of the other masses, in order to insure the absence of any resonant real-emission contributions (see section 5). Such contributions would in fact complicate the comparison among the three codes, which adopt different strategies for handling them, leading in turn to potential non-negligible differences in their predictions.

Degenerate SUSY masses
We report in figure 1 the outcome of the comparison among the LO and NLO predictions of MG5 aMC, Prospino2, and Resummino, for a degenerate SUSY mass setup and the processes of eq. (3.1). As expected, we find no dependence on M others for inclusive stop (top right panel) and slepton (bottom right panel) pair production cross section at the LO, by virtue of the absence of internal SUSY particles of different flavours in the corresponding four-point tree-level amplitudes. This contrasts with the production of a pair of gluinos (top left panel) and charginos (bottom left panel), which both feature production modes that involve t-channel exchanges of squarks with different flavours.
The enhancement with respect to the LO results due to higher-order corrections is very significant for QCD-mediated processes (O(100%)), and considerably milder for the electroweak processes (O(10%)). One striking feature of the M others dependence of the NLO cross sections for botht 1t 1 andgg production lies in the characteristic kink appearing at M others 1330 GeV, which originates from the "resonant" anomalous thresholds [209] 200 400 600 800 1000 1200 1400 production at 13 TeV LHC 200 400 600 800 1000 1200 1400 In the degenerate mass regime, one expects to find complete agreement among the three codes, which is what figure 1 basically shows at both the LO and the NLO, except in the case of chargino pair-production at the NLO, where large differences can be seen between any two predictions. We point out that, for the other processes, the agreement between MG5 aMC and Prospino2 is at the level of 0.5% at the worst (the latter cross sections being larger than the former ones). This may originate from a small mismatch in the input parameters, whose settings are especially intricate in Prospino2 as many of them are directly hard-coded in different parts of its source code (the value of α S (m 2 Z ) is a prime example of this fact). At the NLO, the Prospino2 results for chargino pair-production also differ by a rather flat offset of 3% relative to the MG5 aMC results. This might again be due to a mismatch in the input parameters, which however we could not track down. Conversely, for chargino production the shape of the dependence on M others of Resummino predictions is significantly different from the MG5 aMC one, and a preliminary investigation of the Resummino code has revealed issues in its SUSY-induced renormalisation of a specific vertex; this will be addressed in an upcoming release of the program by its authors. Given the discrepancies found for this process, we have carried out a fully independent analytic computation of the cutconstructible parts of the corresponding virtual matrix elements, and found perfect point-wise

Non-degenerate SUSY masses
Even in the absence of the issues we have outlined in section 3.2, in the non-degenerate squark regime we do not expect to find complete agreement between MG5 aMC and Prospino2. Therefore, at least for the processes for which we found have agreement in the degenerate scenario, we can assess the quality of the mass-averaging procedure implemented in Prospino2 to derive the NLO cross sections when the spectrum is non-degenerate. Figure 3 shows that this approximation can lead to differences with respect to the exact results (as computed by MG5 aMC) of several percent. The exception is stop-pair production, where the effect of the SUSY masses being non-degenerate is small, as they appear only in the virtual amplitudes (and not at the level of real-emission diagrams).
As far as Resummino is concerned, this program has been designed to deal with the dependence on arbitrary SUSY mass spectra in an exact manner. In spite of this, we do not find agreement when comparing its predictions with those of MG5 aMC for chargino-pair production, whilst the agreement in the slepton-pair production case has only been found after a couple of bug fixes in Resummino (that have been implemented in version 2.0.2; unfortunately, these do not address the issue with charginos, which is still under investigation.

Summary of the comparisons
The comparisons presented in this section and the sometimes significant disagreement found among MG5 aMC, Prospino2, and Resummino results, underscore the need for a more comprehensive and robust implementation of NLO QCD corrections for SUSY processes, that can be reliably used by collider experiments. We believe that this is what is achieved currently by MG5 aMC, thanks to its highly automated approach, and its history of orthogonal cross-checks from applications and validation of the same framework to other models and simulations.

Total rates for supersymmetric benchmark processes in simplified scenarios
In this section we calculate the total cross sections at the NLO in QCD for several supersymmetric processes in the context of simplified models, that are typically employed for the interpretation of SUSY searches at the LHC. In these scenarios, one assumes that only the final-state SUSY particles are relatively light, while all of the other superpartners are decoupled by their very large masses. This setup allows us to avoid to deal with intermediate resonances in the real-emission contributions, which will be extensively discussed in section 5.

Setup of the calculation
We consider the processes pp →XỸ (where we denote byX andỸ two SUSY particles, which may also be identical) at the √ S = 13 TeV LHC, at its high-energy upgrade with √ S = 27 TeV, and at a potential future proton-proton collider, identified as the FCC-hh, with √ S = 100 TeV. IfX andỸ are of different species, we enforce their masses to be equal, mX = mỸ . All of the other SUSY particles which do not appear in the final states are decoupled by setting their masses equal to 15 TeV (30 TeV, 110 TeV) when √ S = 13 TeV (27 TeV, 100 TeV, respectively), with the exception of the two stop states whose masses are fixed to mt 1 = 16 TeV (32 TeV, 120 TeV) and mt 2 = 17 TeV (34 TeV, 130 TeV). We refer to appendix C for details on the complexity of a numerical implementation of the decoupling of heavy SUSY particles. The widths of particles are taken to be equal to zero, and the central values of the factorisation and renormalisation scales are set as follows: The theoretical uncertainty stemming from the missing higher-order corrections is estimated by a nine-point independent scale variation, (µ F , µ R ) = (ξ F , ξ R )M , with ξ F/R = 1/2, 1, 2. We use the NNPDF30 nlo as 0118 set of parton densities [210] as provided by the LHAPDF6 interface [211].

Production of a pair of SUSY particles of the same species
In this section, we consider six pair-production processes in whichX andỸ are of the same species. NLO total rates are shown in the left panel of figure 4, for √ S = 13 TeV proton-proton collisions resulting in the production of a pair of gluinos (gg, black circles), a pair of light stops (t 1t 1 , red diamonds), a pair of left-handed up squarks (ũ Lũ L , yellow triangles), a pair of left-handed and righthanded selectrons (ẽ + Lẽ − L andẽ + Rẽ − R , blue circles and brown triangles, respectively), and a pair of (opposite-sign) charginos (χ + 1χ − 1 , green triangles). In order to improve the visibility of the different curves, we have included a rescaling factor equal to 0.1 in the case ofũ Lũ L production. Our results include the bands associated with the theoretical uncertainties obtained from the independent variations of the renormalisation and factorisation scales, as well as from the PDF uncertainties; these are added in quadrature.
The NLO cross sections are found to span about 10 orders of magnitude when the SUSY mass M varies from 100 GeV to 3 TeV, with the strong production of squark and gluino pairs larger than the electroweak production of slepton or electroweakino pairs by orders of magnitude, for any giveñ M value. We also show the associated K-factors in the right panel of the figure, where each K-factor is defined as the ratio of the NLO total rate over the corresponding LO one evaluated at the central scale and with the central set of PDF. The depicted uncertainty therefore reflects the standard NLO cross section uncertainty, as extracted relative to the central LO predictions. The K-factors exhibit different behaviours for the different processes. Firstly, they are larger in the cases of strong squark and gluino production (K ∼ 1.5) than for Drell-Yan-like slepton and chargino production (K ∼ 1.2), as is expected from the strong/electroweak nature of such processes. Secondly, the Kfactor associated with thegg process shows a significant dependence on the SUSY massM , which can be traced back to the virtual amplitudes associated with the quark-antiquark contribution to the cross section and the large gluino colour charge [93]. Whilst subdominant for small SUSY masses, the quark-antiquark contribution becomes significant whenM increases (since the relative weight of the corresponding parton luminosity increases with respect to the gg one), and it therefore impacts the cross section to a more significant level. Conversely, theM dependence of the two Kfactors associated with the production of a pair of squarks is more moderate, and almost absent in the case of the electroweak processes. The right panel of figure 4 illustrates the benefits of higher-order calculations, as it shows that theoretical systematics are smaller at the NLO (filled areas) than at the LO (hashed areas). However, predictions relevant to largeM values are affected, both at the LO and the NLO, by significant uncertainties. This is because, in this region, the latter are dominated by PDF errors. In fact, by increasingM , the average Bjorken x's that enter the partonic cross sections also grow, and at large x's the PDFs are poorly constrained. Fortunately, one expects that a stronger constraining power of the searches for SUSY signals will go hand in hand with better-quality data for SM processes at large scales, which in turn will help reduce the PDF uncertainties.
Similar results as in figure 4 can be found when the centre-of-mass energy is set equal to √ S = 27 TeV (see figure 5) and 100 TeV (see figure 6). As is expected, the main difference with respect to the 13 TeV case is the increase of the cross sections at any givenM , due to the larger available centre-of-mass energies. By scaling up the SUSY massM to match the collider energy, the behaviours of the K-factors are essentially identical for the three collider scenarios.

Production of a pair of SUSY particles of different species
We now consider the production of two SUSY particles of different species, while still setting their masses equal to a common valueM . In figure 7, we present the dependence of the NLO cross sections onM for nine different SUSY pair-production processes, in proton-proton collisions at √ S = 13 TeV. We focus on two strong processes in which a gluino is produced in association with a left-handed up squark (gũ L , red diamonds) or antisquark (gũ L , green triangles), as well as four semistrong processes corresponding to the production of a gluino and a neutralino (gχ 0 1 , black circles), a gluino and a chargino (gχ + 1 , brown triangles), a left-handed up squark and a neutralino (ũ L χ 0 1 , yellow triangles) and a left-handed up squark and a chargino (ũ L χ − 1 , magenta squares). Our results finally also include predictions for three electroweakino pair-production processes in which the lightest chargino is produced in association with the lightest neutralino (χ + 1χ 0 1 , turquoise squares) or with the next-to-lightest neutralino (χ + 1χ 0 2 andχ − 1χ 0 2 , purple diamonds and red pentagons). For all of our predictions, the lightest neutralino is taken to be bino-like, whilst the next-to-lightest neutralino and the lightest chargino are both taken to be wino-like. Analogously to what has been done previously, we present the corresponding K-factors on the right panel of figure 7, and we include theoretical errors estimated by summing in quadrature the uncertainties stemming from scale variations and the PDF errors.
Gluino-squark cross sections (i.e. with agũ L or agũ R final state) are identical to each other, these processes being driven by strong interactions that are blind with respect to the (s)quark chirality. Owing to the larger up-quark density in the proton with respect to the antiup-quark one, the corresponding conjugate processes are suppressed by factors that range from a few units (for small SUSY masses) to almost two orders of magnitude (for large SUSY masses). This is illustrated in figure 7 forgũ L production (the cross section forgũ R is identical to the latter one, and is not shown). The K-factors associated with all these strong processes depend significantly on the SUSY massM , and vary from about 1.5-1.6 for smallM 's to about 1.9 for multi-TeVM . As was already observed in section 4.2, it is the gluino with its large associated colour charge that drives this dependence of the NLO K-factors. Theoretical systematics follow the same pattern as those relevant to same-species production, namely scale uncertainties are reduced at the NLO, whilst the total uncertainty increases for largeM because of the PDF behaviour. This is especially significant in the case of the antiquark density, so that we accordingly find a larger uncertainty for gũ L production than forgũ L production.
The other processes under consideration are of either purely electroweak or semi strong/electroweak nature. The rates are consequently reduced by several orders of magnitude, and the Kfactors turn out to be smaller than for the purely strong processes. As is shown in the right panel of figure 7, QCD corrections are in general larger for the semi-strong processes featuring a gluino or a squark in the final state (K ∼ 1.5), as is expected from their sensitivity to strong interactions that is present already at the tree level. Conversely, the purely-electroweak electroweakino pair production processes exhibit smaller K-factors of about 1.2, a typical value for the Drell-Yan-like electroweakino pair-production that occurs when all squarks are decoupled. The mass dependences of the K-factors are moreover modest, with the exception of the peculiar behaviour exhibited by thegχ + 1 ,χ + 1χ 0 1 , andχ + 1χ 0 2 processes forM ∼ 2.5 TeV. The dominant contribution to these three processes originates from a ud initial state. However, the NNPDF densities used here are mostly un- known at large x (x > 0.1) and are therefore associated with a large uncertainty. Furthermore, the NNPDF methodology (which relies on neural networks to perform the PDF fit) yields the odd shape of the cross sections and K-factors in this regime. Correspondingly, the PDF uncertainties related to these processes grow out of control forM > 1.5 TeV, and the shape of the central K-factor, in particular close toM ∼ 2.5 TeV, stems from the Born and real-emission contributions being affected differently by the corresponding partonic luminosities.
In order to better understand these peculiar features of thegχ + 1 ,χ + 1χ 0 1 andχ + 1χ 0 2 cross sections, we show in figure 8 the predictions obtained when the matrix elements are convoluted either with the NNPDF 3.0 NLO PDF set, or with the CT14nlo Hessian PDF set [212]. While the cross sections evaluated with CT14 PDFs show a more reasonable behaviour at largeM , this comes at the cost of introducing a theoretically-dominated bias on the predictions, as such a PDF set relies entirely, in the large-x region where there is no data point to constrain the fit, on the extrapolation of its parameterisation at the initial scale.
The results obtained by increasing the centre-of-mass energy to √ S = 27 TeV and 100 TeV are presented in figures 9 and 10, respectively. As far as the relative comparisons of these predictions with those relevant to √ S = 13 TeV is concerned, the same observations made at the end of section 4.2 for the case of same-species pair production apply here.
We conclude this section by pointing out that tables reporting the numerical values that correspond to the cross sections shown in figures 4-10 are provided as ancillary files on the electronic archive.

Perturbative computations in the presence of resonances
In this section, we discuss in general the problems posed to perturbative computations by the presence of narrow resonances, and outline the strategies (called Simplified Treatments of Resonances or STR for short) which one may employ to overcome such problems. STR include all of the Diagram Subtraction (DS) and Diagram Removal (DR) procedures defined so far in the literature, and must be seen as a systematic generalisation of the so-called on-shell subtractions (OSS). We also document here the implementation of the STR in MG5 aMC. Illustrative examples of their applications are given in section 6.

General features
producedures In any theory with a sufficiently rich particle spectrum, there is the possibility that the cross section for the production of a given asymptotic state δ is ill-defined in perturbation theory beyond the LO. Here, δ is such that its four momentum can, at least in principle, be reconstructed through measurements performed with a realistic detector, either directly or indirectly through its decay products 2 . This situation can occur in the following case. Let be an LO contribution to the production of δ; a and b denote the incoming partons that initiate the hard scattering, and X a set of final-state particles. The cross section for the process of eq. (5.1) may be inclusive or exclusive in X. At the NLO, real-emission corrections will receive contributions from processes that one can write as follows: where the nature of γ depends on the underlying theory whose perturbative expansion is considered. For example, in QCD γ can be a massless quark or a gluon, while in QED it can be a photon. Let us now suppose that a particle β (which must not resonantly contributes to the process of eq. (5.1)) exists, such that the two-body decay channel is kinematically allowed and that a + b −→ β + X , is a well-defined hard process, that we call the underlying resonant process. Equations (5.3) and (5.4) imply that, among the Feynman diagrams contributing to the process of eq. (5.2), there will be β-resonant ones, namely those that feature a propagator associated with β. In turn, this allows one to write the contributions of such diagrams to the differential cross section associated with eq. (5.2) as follows: dσ ab→δγX where the first term on the r.h.s. is the cross section for the production of an on-shell β. In eq. (5.5), we have denoted by m δγ and m β the invariant masses of the δγ pair and of the particle β, respectively. Even if one Dyson-resums the propagator in eq. (5.5), thus introducing a regularising Γ β factor that prevents the propagator from diverging at m δγ = m β , it may still happen that dσ ab→βX dσ ab→δX . In this case, the NLO contribution due to eq. (5.2) will be numerically (much) larger than its LO counterpart of eq. (5.1), thus 'spoiling' the perturbative expansion of the cross section for δ production.
Situations of this kind are annoying because potentially relevant for phenomenology, especially in SUSY theories where they are ubiquitous. Examples stemming from QCD corrections include in the SM (tW − associated production, whose underlying resonant process is tt production) and (δ, γ, β, X) = (q, q,g,q) , (5.8) (δ, γ, β, X) = (χ, q,q,q) (5.9) in SUSY (squark-pair and squark-neutralino production, whose underlying resonant processes are squark-gluino and squark-pair production, respectively 3 ). A further case is that of a simplified dark matter model, achieved e.g. by extending the SM with a dark matter particle χ and a mediator Y , so that (δ, γ, β, X) = (χ, q, Y, χ) (5.10) is the analogue of the processes of eq. (5.7)-(5.9). By adopting a commonly-used expression, which is strictly speaking incorrect but conveys the basic physics idea, one says that δX production interferes with βX production beyond the LO. The numerical dominance of the latter over the former implies that the corresponding cross section is, to a good approximation, a meaningful physical quantity (for example, we are used to talk about measurements of the tt cross section, which we compare with their perturbatively-computed counterparts). Conversely, the answer to the question of whether it is possible (and, if so, whether it is sensible and/or convenient) to deal with a perturbative, non-β-resonant, δX cross section depends on the context in which one is working. One can introduce two conceptually different classes of applications: 1. Definition of the non-β-resonant δX cross section as a measurable quantity, for a direct comparison with experimental results.
2. Use of non-β-resonant δX production in conjunction with βX production, as perturbative tools that render technically easier the computation of the δX cross section that includes both resonant and non-resonant contributions.
As a rule of thumb, applications of class 1 are mostly of interest to SM physics, while those of class 2 are relevant to both the SM and to SUSY (and, in general, to theories which are not confirmed experimentally, and whose signals need to be searched for). By using again the SM example of eq. (5.7), its class 1 applications entail the definition of the tW − cross section (see e.g. refs. [213,214] for recent ATLAS and CMS measurements of this quantity at the LHC), for which the underlying tt resonant process is considered as a background. On the other hand, in a typical class 2 application one would exploit the possibility of computing both the tW − and tt cross sections at the NLO (i.e. up to O(αα 2 S ) and O(α 3 S ), respectively) for a phenomenologically accurate description of W + W − b(b) production which is much less demanding, from a computational viewpoint, than the calculation of the W + W − bb cross section at the NLO.
The key point of the previous example is that the tW − cross sections that enter the two applications are not necessarily defined in the same way. In general, let be the amplitude associated with the process of eq. (5.2). The two quantities on the r.h.s. of eq. (5.11) denote the non-β-resonant and the β-resonant contributions, respectively. The matrix element will thus be proportional to: For both class 1 and class 2 applications, the contribution of the last term on the r.h.s. of eq. (5.12) must be minimised. In the context of NLO+PS simulations, this problem has been solved in ref. [215] by introducing two different types of procedures. In diagram removal (DR), one simply drops this contribution, whereas in diagram subtraction (DS) the non-β-resonant δX cross section will feature the linear combination: The pre-factor in the second term of eq. (5.13) is arbitrary to a large extent, but must obey the condition lim 14) The symbol P denotes a kinematic projection that maps a generic δγX configuration onto one that has m δγ = m β . Crucially, this map is fully local in the phase space, so that the difference in eq. (5.13) vanishes identically when m δγ → m β (also thanks to eq. (5.14)); such a locality condition is essential for the use of DS in event generators. Owing to the definition of the β-resonant amplitude, one has by neglecting production spin correlations. Thus, eq. (5.15) renders it intuitively clear that the difference in eq. (5.13) is constructed so as to avoid the double counting of the LO βX cross section in class 2 approaches. In practice, spin correlations cannot be neglected, and therefore A ab→βX is never used as such in DS procedures 4 ; this is just as well, since it helps guarantee the local cancellation between the two terms in eq. (5.13). We finally stress that there is ample freedom in the choices of the function f and projector P that enter the definition of a DS cross section through eq. (5.13). We shall exploit this fact in the following, by considering several different implementations.
As far as the second term on the r.h.s. of eq. (5.12) is concerned, DR procedures do not include it in the definition of the non-β-resonant δX cross section, while DS procedures do. Thus, for class 1 applications it is essential that both DR and DS results be obtained, and that their difference be less than the theoretical systematics. If that is not the case, non-β-resonant δX observables are simply not physically meaningful, and both DR and DS predictions must be discarded. Conversely, for class 2 applications, in which the emphasis is on obtaining the best approximation for the full δX cross section, one is interested in using DS approaches. DR results might also be kept, provided they are statistically compatible with the DS ones. We also point out that in the NLO+PS simulations of ref. [184] a third scenario has been considered (dubbed DR2 there, and originally introduced in ref. [116] within a fixed-order calculation), in which one keeps the second term on the r.h.s. of eq. (5.12), but does not perform the subtraction of eq. (5.13). As far as its usage in applications of class 1 and 2 is concerned, this approach is quite analogous to a DS one. However, given that no subtraction is carried out, we call it DR+I (for diagram removal plus interference).
In summary, the non-β-resonant δγX cross section can be defined as follows: (5.18) in the DR, DR+I, and DS approaches, respectively. Collectively, the procedures implied by eqs. (5.16)-(5.18) will be called STR (that stands for Simplified Treatments of Resonances). STR strategies have been pursued for a long while in the context of fixed-order calculations and for inclusive observables in both the SM and BSM theories (in particular in SUSY, where they are typically called OS subtractions) -see e.g. refs. [93,116,[217][218][219][220][221][222]. None of these earlier procedures is apt to be applied to exclusive event generation, and thus we believe one should refrain from using the DR or DS tags in association with them. As far as DR and DS procedures are concerned, either identical to or featuring variants of those originally proposed in ref. [215], results can be found in refs. [142,143,148,149,184,[223][224][225].
We now turn to giving some details about the implementation of eqs. (5.16)-(5.18) in MG5 aMC. In keeping with the general strategy that underpins the code, everything is fully automated and process-as well as model-independent 5 . We remind the reader that MG5 aMC is self-consistent, and thus that, in particular, it generates internally the Feynman diagrams and writes the corresponding amplitudes. This implies that the code stores the information on the topological structure of each diagram, and therefore knows where to find the resonances (which is essential in order to construct A ( ¡ β) ab→δγX and A (β) ab→δγX ). Furthermore, it can control the kinematical inputs and parameter settings in a diagram-by-diagram manner if needed. The immediate consequence of the previous observation is that the construction of the DR and DR+I cross sections of eqs. (5.16) and (5.17), respectively, is achieved in a straightforward (and unique) manner.
The case of the DS cross section is more involved, owing to the freedom in the definitions of the function f and of the projector P, although after having chosen f and P, eq. (5.18) uniquely determines the corresponding DS procedure. Unfortunately, it is impossible to parametrise the arbitrariness in the choices of f and P, and thus one must limit oneself to considering a finite number of physically-motivated options. We describe those implemented in MG5 aMC below, and point out that previous results in the DS approach [142,143,148,149,215,[223][224][225] have been obtained with a given (f, P) pair (with the exception of ref. [184], where two different forms of f have been compared).

Diagram-subtraction procedures
We start by pointing out that the discussion that follows is relevant to the last term on the r.h.s. of eq. (5.18), henceforth called the DS subtraction term. The other three terms in the definition of the DS cross section correspond to a straightforward tree-level calculation, and are thus not of concern here. We denote the kinematic of the process of eq. (5.2) as follows: where we have assumed that the set X is composed of n − 1 particles with momenta k i . It is convenient to introduce the following auxiliary momenta: Although the resonance β does not appear in the final state, the definition of its momentum in eq. (5.21) is physically meaningful, since we are solely dealing with β-resonant diagrams. In the centre-of-mass frame of the incoming hadrons: (1, 0, 0, 1) and with S being the squared hadronic centre-of-mass energy. Its parton-level counterpart reads thus: The action of the projection P on the partonic kinematic configuration is denoted as the following transformation: We also introduce, for consistency with eq. (5.25), the momentum of the projected resonance The DS strategies that we consider generally require the partonic incoming momenta to be changed. This can formally be seen as also stemming from the action of P, and thus be represented as follows: While the specific form of eq. (5.28) will depend on P, in all of our implementations we shall always choosex a andx b so thatx This implies that the original and projected partonic centre-of-mass frames will travel at the same speed w.r.t. the hadronic centre-of-mass frame. Equation (5.28) has two further implications. Firstly, the flux factor of the DS subtraction term is equal to 1/(2s). Secondly, its parton-luminosity factor is given by with H 1 and H 2 being the incoming hadrons. As far as the function f is concerned, we shall limit ourselves to considering the following form: is a generalised Breit-Wigner function (the standard one being obtained by setting x = m β ) in which κ is a normalisation factor that does not play any role. The rationale beyond eq. (5.31) is that its denominator will cancel, to some extent, an analogous term implicit in the projected matrix element of the DS subtraction term, that is thus replaced by the numerator of eq. (5.31) which supposedly models some of the off-shell-β effects. The precise extent of such a cancellation depends on the interplay of several factors (such as the choice of x, of the operator P, or of the PDFs), and cannot therefore be predicted a priori. This is one of the reasons why in eq. (5.31) x is treated as a free parameter.
In view of their use in the definition of P, we also introduce the following quantities. For any four-momentum p, we denote the boost to its rest frame by: This understands that p 2 = m 2 > 0; we implicitly assume that the boost is performed along p, and that it is such that: If in the rest frame of p we impose a 1 → 2 four-momentum conservation, with: We can now present specific details of the definition of the DS subtraction term in the DS procedures we pursue. The reader must keep in mind that the kinematic configuration of eq. (5.19) and its associated quantities eqs. (5.20)-(5.22) are thought to be given. Without loss of generality, we work in the incoming-parton centre-of-mass frame, q = 0, and the options described below are associated in MG5 aMC with the function f given in eq. (5.31) and with either of the settings: Other choices of x would be straightforward to implement.

Option A
This option follows the strategy first introduced in ref. [215]. We define: The momenta not associated with the resonance β are left invariant: In view of eqs. (5.43) and (5.44), we then define: As far as the momenta of δ and γ are concerned, we proceed as follows. First, we boost them in the rest frame of k β : Then, by keeping the information on e δ but discarding all the rest, we definek

Option B
This option generalises (to an arbitrary number of final-state particles) the strategy of ref. [142]. We define the mass of the recoil system X: Next, we define the energies of the projected resonance and recoil system as follows: The corresponding three-momenta are defined by preserving the direction of the original threemomenta, rescaling their lengths so as to impose the mass shell conditions These guarantee that k β = − k rec , since k β = − k rec . We also definē Finally, k δ and k γ are projected using the same procedure as in eqs. (5.46)-(5.49).

Option C
This option, which is currently not implemented in MG5 aMC, follows closely what is done for the phase-space generation as is carried out in the module MadFKS [155], in the case relevant to a massive FKS sister. This, in turn, generalises the massless-parton treatment of ref. [228]. The projected partonic centre-of-mass energy is defined as follows: with ς ≥ 1 being a free parameter 7 , typically of O(1). One then regenerates the final-state kinematic configuration, using the same random numbers as for the original one, ands instead of s (the configuration thus obtained has only a temporary role, and is denoted below in the same way as the original one). Next, a boost B is defined along the resonance three-momentum k β , so that: The projected momenta are then given by: Once again, the k δ and k γ momenta are projected using the same procedure as in eqs. (5.46)-(5.49).

Using DR and DS in MG5 aMC
All of the STR procedures described above can be employed within MG5 aMC by downloading the MadSTR plugin 8 , and by copying the directory MadSTR thus obtained inside the PLUGIN directory, which is part of any (recent) MG5 aMC release. The current version of MadSTR is compatible with MG5 aMC version 2.6 and higher; compatibility with versions 3 and higher, that are capable of carrying out mixed-coupling perturbative computations [70], will be added in the future. The plugin can be activated by using the following command to start MG5 aMC: istr STR procedure  Table 2: Possible values for the istr parameter, with the corresponding STR procedure. In the case of option A (and of option B when the condition in eq. (5.51) is not fulfilled), the parameters str include flux and str include pdf are also relevant (see text for details).

mg5_aMC --mode=MadSTR
One can then generate NLO processes and write them to disk as usual (with the generate and output commands -see ref. [56] for more details). The plugin will take care of identifying any potentially resonant contributions, of generating the associated underlying resonant processes, and of taking care of the extra bookkeeping in addition to that normally performed by MG5 aMC in nonresonant cases. Depending on the mass spectrum, the contributions for which the STR is needed are identified at run-time. Three parameters found in run card.dat, namely istr, str include flux, and str include pdf, can be used to choose the desired STR procedure and its associated options. More specifically, the type of STR is controlled by istr, which must assigned an integer value according to the options given in table 2. The other two parameters, str include flux and str include pdf, are active only if the STR is of DS type, and control the settings of the flux and the parton luminosity factors, respectively, in the DS subtraction terms. In particular, DS procedures imply changes to the partonic centre-of-mass energy (see eqs. (5.45) and (5.53)). In turn, this seemingly implies that the flux and the luminosity factors should be changed accordingly, when the partonic centre-of-mass energy is changed. However, this is not mandatory, since it is actually part of the definition of the projection inherent to DS procedures and, as such, is liable to be chosen by the user. It is in order to give one the possibility of exploring the consequence of this choice that the parameters str include flux and str include pdf have been introduced. By setting them equal to True (False), the flux and PDF factors are (not) re-defined. Note that the settings of the two parameters are independent from each other, and that True are the default values. More details in the context of a specific example will be given in section 6. Finally, the value of Γ β (which acts as a regulator when m δγ → m β ) can be controlled by changing the width of the corresponding particle β in the param card.dat file, for all β's that are potentially resonant. The code will set the widths of all coloured particles 9 equal to zero everywhere except in the resonant real-emissions diagrams and in the corresponding DS subtraction terms, in which the values provided in the param card.dat will be employed.

A case study: jets plus missing energy at the NLO+PS accuracy
We are now in the position to perform phenomenology studies in the MSSM with a generic particle mass spectrum. As an illustrative example, we choose the benchmark point presented in table 3, which is not excluded by current experimental searches at the LHC, and that features non-trivial decay patterns. In contrast with section 2, the bottom squarks are taken to be non-mixing. 5500 GeV mχ ± ,˜ ∓ 5500 GeV    table 3. q denotes the down-type (up-type) quark associated with the same-generation up-type (down-type) quark q.
In the scenario of table 3, the total widths and the relevant decay modes of the gluino and the squarks are as reported in table 4, together with the associated branching ratios, these results being computed at the LO with MadWidth [226]. Thus, the decay widths of the squarks and gluino are sufficiently small relatively to their masses (except fort 1 , which is however never resonant in the processes considered in this section) so that the narrow-width approximation is sensible. Conversely, and according to the parametrisation of long-distance effects e.g. in Pythia8.2 and Herwig (see refs. [229][230][231]), the sparticle widths are sufficiently large to avoid hadronisation before decay -in other words, no R-hadrons will be present in our simulations.
We shall now consider, at the 13 TeV LHC, multijet plus missing transverse-energy final states:  Table 5: Contributions to the signatures of eq. (6.1) stemming from the processes of eqs. (6.2)-(6.4). For each of these, we report the Born-level signature (second column), the underlying resonant process (third column), and the decay of the would-be resonant sparticles (fourth column). a signature that is typical of SUSY searches at hadron colliders. We shall compute the contributions to eq. (6.1) due to the following underlying processes: The subsequent hadronic decays of the final-state sparticles are carried out by the parton shower programme according to the results of  Table 6: Total inclusive and fiducial cross sections (in fb) at the √ S = 13 TeV LHC. The leftmost errors stem from scale variations, the rightmost ones from PDF uncertainties. We have set str include pdf=True and str include flux=True (see section 5.3).
in figure 11. In order to ensure that squarks only decay into light-flavoured jets, we restrict our simulations by solely considering the first two generations of squarks. We also distinguish, in the rest of this section, the light (anti-)squarkũ ( ) R from all of the other heavier (anti-)squarks of the first two generations, that we denote byq h .
The situation is summarised in table 5. For each of the processes of eqs. (6.2)-(6.4), which we also distinguish according to whether a light squark is present in the primary (i.e. before decay) process, we report the Born-level signature (second column), the underlying resonant process, defined according to eq. (5.4) (third column), and the relevant decays of the primary sparticles, according to eq. (5.3) (fourth column). The experimental signature of eq. (6.1) is obtained by considering the parton-level process definitions of eqs. (6.2)-(6.4) whose Born contributions do not feature any resonance 10 . Their corresponding real-emission contributions therefore include at most one SUSY resonance which is then necessarily subject to an STR procedure. In this way, our simulation setup is guaranteed to remain within the scope of STR applicability as described in section 5.1.
All of our simulations are NLO+PS accurate, whereby NLO matrix elements are matched with the Pythia8.2 [232] parton showers according to the MC@NLO formalism [34], automated in MG5 aMC. The resulting hadron-level events are clustered by making use of the anti-k T algorithm [233] with jet radius R = 0.4, as implemented in FastJet [234].
For our phenomenological analysis, we implement an event selection similar to the one of the CMS SUSY search of ref. [235]. Firstly, jets are required to have transverse momentum larger than 30 GeV, and pseudorapidity |η| < 2.4. Events are kept if they feature at least N jet ≥ 2 jets. Secondly, the total hadronic activity H T , defined as the scalar sum of the transverse momenta of all reconstructed jets, must be larger than 300 GeV. Thirdly, the missing transverse hadronic energy H miss the negative of the vector sum of the transverse momenta of all reconstructed jets with a pseudorapidity |η| < 5, must be larger than 300 GeV. Finally, the two leading jets and − → H miss T are imposed to be well separated in azimuth, ∆φ(H miss T , j 1,2 ) > 0.5. When a third and a fourth jet are within the acceptance defined before, we additionally impose ∆φ(H miss T , j 3,4 ) > 0.4. Total cross sections with (σ fiducial ) and without (σ inclusive ) the above cuts at the 13 TeV LHC are presented in table 6, for each of the six STR procedures listed in table 2. Gluino-pair production leads to the smallest cross sections, as a result of the large gluino mass and the correspondingly small gluon PDF at large Bjorken x's. The NLO inclusive (fiducial) cross section turns out to be equal to about 0.33 fb (0.23 fb), with a large K-factor of about 1.8 both for the inclusive and fiducial cases. This large K-factor originates from the large colour charge associated with the gluino, and the purely strong nature of the Born process, as was already discussed in section 4. Scale uncertainties are reduced by a factor of two with respect to the LO ones and, as a consequence of the typically large Bjorken x values associated with the production of a pair of 2 TeV gluinos, the theoretical error is dominated by the PDF uncertainties. Both scale and PDF errors are essentially independent of the STR procedure adopted, which is the reason why we only report them for the istr=2 case. As far as the STR dependence itself is concerned, it is about 3%, and thus much smaller than the other theoretical uncertainties. As is well known, this is an observable-dependent statement, and we shall show later that at the differential level things are more involved. Finally, we remark that the Monte Carlo integration errors are equal to about 0.2%, and have therefore been ignored.
Because of the smallerũ R mass, the cross sections of gluino-squark associated production and of squark pairs are 25 to 60 times larger than the gluino-pair ones. This behaviour is driven by that of the parton luminosities: valence quarks contribute significantly, since one is in an x-region where their PDFs are large. Correspondingly, the PDF uncertainties are much smaller than for gluino-pair production, as deep-inelastic scattering, fixed-target experiment data, LHC Drell-Yan, forward W -boson, and Z-boson data allow one to strongly constrain the fit of the valence quark densities at large x's. This implies that for these processes, at variance with the case of gluino pairs, the PDF errors are smaller than the scale ones, in spite of the fact that the latter are a factor of three smaller at the NLO than at the LO (as opposed to a factor of two for gluino-pair production). As far as the STR-option dependence is concerned, it is below 1% and thus, as in the case of the gluinos, largely subdominant with respect to the other uncertainties.
In figures 12, 13, and 14 we present six differential distributions for gluino-pair, gluino-squark, and squark-pair production, respectively. Such observables are directly relevant to the CMS SUSY search of ref. [235], and thus all results are obtained by applying the fiducial-volume cuts defined before. Each panel of each figure has the same layout, namely: the main frame presents the differential distributions in pb/GeV or pb, at the LO (dashed black histogram) and at the NLO (solid coloured histograms; there are six of these, corresponding to the six STR options we have considered). The statistical errors are shown as error bars, while the scale uncertainties are displayed as dark-blue bands; finally, the linear sums of the scale and PDF errors are represented as light-blue bands. The lower inset presents the ratios of the six NLO predictions over the LO one, the latter evaluated with the central scale choice and central PDF (i.e., these are the standard K factors).
Mirroring what has been found for total rates, PDF uncertainties are quite substantial for gluino-pair production, as one is always in a kinematic regime where large Bjorken x's are relevant. For the other two processes, such uncertainties are reduced, except in those phase-space corners where again one is sensitive to large Bjorken x's values, e.g. for large H T or H miss T . It is therefore obvious that, despite the progress made in the computations of the short-distance cross sections, precise simulations at large scales can only be achieved by constraining much more strongly the PDFs at x → 1. As far as the K-factors are concerned, they are found to be relatively flat for most observables, and close to 2 forgg production, to 1.5 forgq production, and to 1.4 forqq production, respectively. Distinctly non-flat K-factors are obtained in particular in the cases of the H T and jet-multiplicity distributions. This helps stress the fact that, in general, the re-scaling of LO+PS predictions by an overall constant factor is a dangerous procedure that may lead to unreliable results, which underscores one of the main motivations of the present paper.
It is interesting to observe that the dependence on the STR choice is generally very mild also at the differential level, which gives one confidence on the description of multijet SUSY-induced signatures through the production-and-decay picture of eqs. (6.2)-(6.4). A notable exception is the large-H T region (for all of the three processes, although it is particularly prominent in the case of gluino-pair production), where the STR dependence becomes extremely large, and appears to be pathological. We shall argue that, in fact, this behaviour allows the method to self-diagnose that it is being applied in regions where its founding assumptions are dubious at best. Let us start with a technical explanation. The STR procedures that differ the most from the median are those associated with istr=3,4 -these are DS procedures, option A (see table 2). In these cases, through a reshuffling operation the momenta of the incoming partons are changed, thereby implying that the corresponding Bjorken x's are changed too. This affects the value of both the flux factor and the PDFs. In the notation of section 5.2: , (6.5) which shows the effect on the flux and parton-luminosity factors (first and second terms on the r.h.s. of eq. (6.5), respectively) at the level of the subtraction cross section (last term on the r.h.s. of eq. (5.18)). If one generates a far off-shell kinematic configuration, m δγ m β , then x a x b x axb . Therefore, since eventually all PDFs decrease with increasing Bjorken x's, both terms on the righthand side of eq. (6.5) are large, which implies a strong suppression of the physical cross section of eq. (5.18) in this kinematic region (since the subtraction term is large); so strong, in fact, that it may become negative. Owing to the structure of the function f (m 2 ) given in eq. (5.31), this feature is more pronounced for istr=4 than for istr=3, which clearly shows in the figures.
In fact, although the mechanism we have just described is responsible for the behaviour of the cross sections in the large-H T region, its effects are particularly dramatic in gluino-pair production owing to x a,b 1 there (in other words, if the gluino had a smaller mass, the subtraction cross section would be significantly smaller). For such Bjorken x's, the central values of the PDFs are very suppressed and poorly constrained, and thus affected by very large residual uncertainties: indeed, we see from figures 12-14 that as soon as the STR-choice dependence becomes very large, so does the PDF uncertainty, the effect being larger when parton luminosities involving the gluon density are relevant. We have verified that, by "removing" the luminosity factor from eq. (6.5) by setting str include pdf equal to False (see section 5.3), the results obtained with istr=3,4 are much closer to the others, thus showing that the large STR-choice dependence is driven by the PDFs.
We point out that sizable differences due to STR options should be expected also for m δγ m β . However, while the Breit-Wigner distribution is symmetric around m β , the flux and PDF ratios are not. This implies that the explanation given above may not apply in this kinematic regime. However, regardless of which mechanism is responsible for an enhanced STR-choice dependence, the take-home message is the following: such a dependence is the signal that describing the process of interest by means of a cross section dominated by resonant production of (s)particles which subsequently decay is simply not adequate, and a full (unfortunately, more complicated) computation is required. Ultimately, the decision of where to stop trusting a simplified computation rests with the user, and it depends on many factors, in particular on whether one is interested in a class. 1 or class. 2 approach (see the itemised list below eq. (5.10) in section 5    Figure 13: Same as in figure 12 but for gluino-squark production.  Figure 14: Same as in figure 12 but for squark-squark production.

Conclusions
With the steady increase of the statistics accumulated by the LHC experiments, and the absence of positive results in the searches for new physics, it becomes necessary to improve the accuracy of the simulations of BSM signals, thus matching that of their SM backgrounds. By far and large, this currently means matrix elements computed at the NLO in QCD matched to parton showers (NLO+PS). Such a necessity has spurred some recent theoretical activity, whereby authors have addressed the needs of specific search strategies. The goal of this paper has been that of rendering such an improvement systematic (i.e. achievable for arbitrary processes and for a vast class of renormalisable theories). This is feasible thanks to the powerful and flexible environment constituted by the automated program Mad-Graph5 aMC@NLO and the physics models it can use for simulations. In particular, for the sake of the present work, two major limitations have been lifted. Firstly, at the level of the construction of models, we have overhauled the way in which FeynRules and NLOCT deal with on-shell renormalisation schemes, so that a much larger flexibility is achieved, that helps deal with the difficult cases especially relevant to supersymmetric theories. Secondly, at the level of the Mad-Graph5 aMC@NLO code proper, we have automated a variety of solutions to the problem posed by the presence of partonic channels that appear beyond the leading order in perturbative computations, and that feature narrow resonances; this problem is particularly acute in theories with a complex mass spectrum. Such solutions, that we have dubbed Simplified Treatment of Resonances (STR), generalise the so-called on-shell subtractions, and encompass the Diagram Removal (DR) and Diagram Subtraction (DS) strategies introduced in the last few years in the context of NLO+PS simulations. Technically, these two pieces of work have been implemented, respectively, in a plugin for FeynRules, called MoGRe, and in a plugin for MG5 aMC, called MadSTR.
As a proof-of-concept, we have generated an NLO UFO model for the MSSM with a widely-used renormalisation scheme, and we have studied processes that feature intermediate resonances at the LHC, at the NLO+PS accuracy and with a realistic set of final-state cuts. We have performed thorough self-consistency checks of our implementation, and compared some loop matrix elements generated by the code with those resulting from analytical computations. It is important to bear in mind that the core structure of the MadGraph5 aMC@NLO code, which has been only minimally affected by the present work, has by now been very extensively validated in countless simulations. As a further a posteriori validation, we have compared total-rate leading-and next-to-leading order results for benchmark 2 → 2 processes with those of the public codes Prospino2 and Resummino, with the restrictions that these two programs enforce. Such comparisons, whose details can be found in section 3, in some cases show disagreements (of different origins) among the various predictions, and further underscore the advantages of a general, process-independent, and automated implementation.
In conclusion, with the present work we have achieved, for the first time, the complete automation of NLO+PS simulations for supersymmetric-particle production at hadron colliders in the framework of the MSSM with a generic particle spectrum, and set up the tools for dealing with similarly involved theories by means of a user-driven framework. We point out, however, that we have not yet implemented and tested the general solution, introduced in ref. [70], to the problem posed by unstable resonances in the context of the complex-mass scheme [236,237]. Such a solution requires further developments in FeynRules and NLOCT, so that NLO UFO models contain the necessary routines for dynamically and automatically selecting, according to the particle spectrum, the appropriate Riemann sheets for the calculation of the UV counterterms. This is left to future work.
All of the computer programs relevant to this paper are publicly available -on top of Mad-Graph5 aMC@NLO, the MSSM model, together with the MoGRe plugin and ready-to-be-used Mathematica notebooks, can be found at http://feynrules.irmp.ucl.ac.be/wiki/MSSMatNLO, while the MadSTR plugin can be downloaded from https://code.launchpad.net/~maddevelopers/mg5amcnlo/MadSTRPlugin. This paper is accompanied by ancillary files, stored on the electronic archive, that collect NLO QCD results for total rates of pair-production supersymmetric processes, which we have refrained from including here for reasons of space.

A Conventions for one-point and two-point functions
In the analytical formulas presented in this paper, all A and B loop-integrals have been normalised as where we recall that D = 4 − 2 is the number of spacetime dimensions and µ R is the regularisation scale (taken equal to the renormalisation scale). The B µ vectorial integral has been further reduced to a scalar integral using Lorentz covariance, and the B 1 integral is connected to several B 0 integrals as the p 2 → 0 limiting case involving a derivative of the B 0 function with respect to the p 2 variable instead of the squared bracket. Explicitly, one gets and with the ultraviolet-divergent parts of the integrals being written in terms of the number of spacetime dimensions and the Euler-Mascheroni constant γ E , 1 = 1 − γ E + log 4π. Several special limits for the B 0 function and its B 0 derivative are useful, We stress that since this paper does not consider the complex-mass scheme, the renormalisation counterterms are defined using only the real part of the two-point functions. where one assumes that the current directory is the one containing the plugin.

B.1 The main method MoGRe$Renormalize and its options
The main function of the MoGRe plugin is called MoGRe$Renormalize and takes a Lagrangian as an input, as for instance in the following example where LMSSM stands for the MSSM Lagrangian. The user is allowed to specify three options that modify the behaviour of the method and that respectively address the treatment of the four-scalar vertices, loop-induced field mixing and the nature of the interaction in which the loop-corrections are evaluated. These options can be set following a standard Mathematica syntax, as for instance through the command This first indicates, through the Exclude4Scalars option, that all model four-scalar interactions have to be ignored in the renormalisation procedure, so that the corresponding counterterms are not evaluated. The default choice for this option (set to True in our example) is False. While strictly speaking ignoring the renormalisation of the four-scalar vertices is incorrect, these vertices rarely appear at tree-level so that the associated counterterms are often not necessary. Avoiding their calculation and their inclusion in the final UFO model therefore allows for an increase of the efficiency of the computations, both at the NLOCT and MG5 aMC level.
Secondly, we provide information on the different sets of fields that mix at the one-loop level through the FlavorMixing option of the MoGRe$Renormalize method. If set to True (the default choice), all fields carrying the same quantum numbers and lying in the same spin and colour representations are assumed to mix. In contrast, all loop-level mixings are forbidden if this option is set to False. The user has also the possibility, like in our example, to provide a list with the different sets of fields that mix at the one-loop level. Any physical field can be used in such a list, and the code further checks whether the input is compatible with the representation of the involved fields under the model gauge groups. In our example, we have forbidden any loop-level mixing, except the one of the two stop-eigenstates (denoted st1 and st2 in the FeynRules model) and the one of the two sbottom-eigenstates (denoted sb1 and sb2 in the FeynRules model), as given by eq. (2.22).
Finally, the last option indicates which type of interaction should be renormalised, among all the interactions declared in the MR$InteractionOrderHierarchy option of the FeynRules model implementation. In the MSSM implementation, two types of interactions are available, namely the QCD and QED ones. In the above example, that matches the physics goals of this paper (NLO QCD corrections for the MSSM), we have solely selected the QCD interaction type QCD.
Before describing how the bare Lagrangian is technically renormalised, we detail in the next subsections various methods that can be used to simplify the model and modify the way in which MoGRe$Renormalize works.

B.2 Simplifying the procedure
In general, all external parameters have to be renormalised, which yields a heavy renormalisation procedure for complex models like the MSSM. However, some parameters may not need to be renormalised, like the electromagnetic coupling constant that does not receive any correction at one loop in QCD. Whilst this type of information can be useful to speed up the renormalisation procedure, the programme cannot guess it at this stage where no calculation has been done yet. For this reason, the user is allowed to declare the parameters that should not be renormalised with the MoGRe`DefineUnrenormalizedParameters method that takes, as arguments the corresponding symbols as implemented in the FeynRules model. The arguments can also be provided as a unique list. Going back to the considered example where only QCD corrections matter, the command leads to declaring all parameters connected to the electroweak sector (namely the electroweak inputs, the Higgs sector parameters, the electroweak gaugino and scalar soft masses as well as the chargino and neutralino mixing matrices) to be insensitive to QCD corrections at one loop. The renormalisation of all relevant internal parameters is accordingly and automatically simplified through their functional dependence on the above parameters.
Similarly, the code assumes by default that all fields get renormalised, although this may not be the case in practice. For instance, the weak boson two-point functions are insensitive to QCD corrections at one loop. This type of information can be passed to the code by means of the MoGRe`DeclareUnrenormalizedFields method that takes as arguments all fields that should not be renormalised. The arguments can be provided either sequentially or under the form of a list. For instance in the MSSM implementation worked out in this paper, we could use MoGRe`DeclareUnrenormalizedFields[ seL, seR, smuL, smuR, stau1, stau2, sne, snm, snt, A, W, Z ]; although the sleptons are in principle not necessary as they do not appear in any QCD vertex (but we keep them here for illustrative purposes). This prevents all charged sleptons, sneutrinos and electroweak bosons from being renormalised. All remaining fields will be renormalised, the associated wave-function renormalisation constants being taken complex by default. Reality conditions can be enforced through the usage of the MoGRe`RealFieldRenormalisation method, that takes as an argument the symbols associated with the concerned fields. All relevant symbols can be given again either under the form of a list or of a sequence. In the case where no argument is provided, all wave-function renormalisation constants are taken real, as with the following example,

MoGRe`RealFieldRenormalisation[ ]
For mixing fields for which matrix renormalisation is in order, the method acts on all the elements of the renormalisation matrix.

B.3 Restrictions
The FeynRules model implementation may contain (external or internal) parameters that are numerically vanishing when default values are accounted for. While it is in general safer to keep these parameters all along the renormalisation procedure, so that they are renormalised and get associated with potentially non-vanishing renormalisation constants, this is often not necessary and only leads to heavier subsequent calculations. The MoGRe package by default takes care of the removal of these zero parameters both from the tree-level Lagrangian and from the rules dedicated to the replacement of the bare quantities by the renormalised ones. In this way, those parameters will not be renormalised as not present in the bare Lagrangian anymore, and the code will make sure that they do not re-appear through the renormalisation of other quantities. For instance, the CKM matrix is often taken diagonal, so that MoGRe by default removes all its off-diagonal elements. This behaviour can be turned off by issuing the command EnforceZeros = False;

B.4 Specifying the renormalisation scheme
As detailed in section 2.2.4, there is no unique way to define an OS renormalisation scheme in the MSSM and in SUSY in general. This consisted in the main reason that has led to the development of the MoGRe package. Scheme-specific renormalisation conditions can be added by the user by making use of a dedicated method named AddRenormalizationCondition. The latter takes two arguments, a renormalisation constant (associated with either a parameter or a field) and a function of different parameters, fields and other renormalisation constants. As a result, the first renormalisation constant will be considered equal to the function given as the second argument. The way in which field renormalisation constants should be input follows the FeynRules syntax, wave-function renormalisation constants being provided as for non-fermionic, left-handed fermionic and right-handed fermionic fields respectively. In the diagonal case, the two field symbols fld1 and fld2 are equal. In a non-diagonal field mixing case, they can be different. For a parameter prm, the corresponding syntax reads In a second step, still prior to any computation, simplifications are performed and the Lagrangian is put under an internal format allowing for a more efficient run. More precisely, the Lagrangian is truncated from its constant terms, and all parameters that are vanishing are removed except if the EnforceZeros flag has been set to False (see section B.3). The Lagrangian is then expanded so that all (unphysical) gauge eigenstates are replaced by physical mass eigenstates. In the case where four-scalar interactions are requested not to be renormalised (see section B.1), they are removed from the Lagrangian. They will however be reintroduced at the very end of the renormalisation procedure, as those interactions could appear into given loop diagrams. Moreover, MoGRe requires that all tree-level bilinear terms are canonically normalised and that all kinetic and mass mixings have been appropriately rotated away by the user. In practice, the code ignores all bilinear terms provided by the user and reintroduces the canonical ones directly on the basis of the model field content.
As a last initialisation step, MoGRe makes use of the FeynRules model information on the external and internal parameters that must be exchanged during the renormalisation procedure (the FeynRules FR$LoopSwitches option [58]). In the FeynRules syntax, FR$LoopSwitches consists in a list of 2-tuples of parameters, FR$LoopSwitches = { { prm1 , prm2 } , .... } in which prm1 is external and prm2 is internal. Prior to the renormalisation of the model, prm2 is made external and prm1 internal, the dependence of this last parameters on the other parameters being derived by the code. For instance, in many publicly available FeynRules models, the Wboson mass m W is derived from the other electroweak inputs. OS renormalisation however requires m W to be external. The FR$LoopSwitches is then used to trade it, for example, with the Fermi constant G F . We refer to the FeynRules manual [58] for more information.
Some parameters can moreover be doubly defined in FeynRules models, like the Yukawa couplings and the fermion masses in the SM that are taken as different input parameters, although being numerically equal. The idea behind this trick consists in allowing for massless fermions and non-zero Yukawa couplings at tree-level. However, this makes the theory ill-defined when renormalisation is at stake, as all these parameters must be enforced to be equal when counterterms are evaluated. This can be achieved by means of the FR$RmDblExt FeynRules parameter, that consists in a list of Mathematica replacement rules mapping one of the doubly-defined parameters to the other. For instance, replaces every single Yukawa coupling (normalised to be exactly equal to the associated mass parameter) by the corresponding fermion mass. Such a replacement is also enforced in MoGRe in the case where the FR$RmDblExt parameter exists. Concerning the model introduced in this work, all parameters are uniquely defined so that this is irrelevant. Along these lines, we recall that we rely on NLOCT [59] for the analytical computation of the various counterterms. This implies that care must be taken with any coupling depending on particle masses, like the Yukawa couplings or the trilinear scalar interaction strengths of eq. (2.17). By default, NLOCT renormalises them in the MS scheme regardless of any finite pieces that are relevant in the OS scheme. A correct treatment therefore requires to replace them by their analytical expression (that involves the fermion masses), prior to the call to NLOCT. This can be achieved straightforwardly with the RemovingInternalCst method of MoGRe that re-expresses a given parameter in terms of the others. Concretely, this method removes a given parameter and the associated renormalisation constant from the model and replaces them by the corresponding analytical expressions. Whilst the replacement rule of the parameter itself is taken from the model file, the one of the renormalisation constant is derived on the fly.
In the context of the MSSM implementation presented in this paper, we have implemented so that all Yukawa and trilinear scalar interactions have been replaced according to their dependence on the external quark masses. We have also made use of this method to define the renormalisation of the strong coupling through α S and not g s .

B.6.2 Field renormalisation
In order to get the list of fields that should be renormalised, the programme starts by extracting from the Lagrangian all relevant interaction terms on the basis of the information passed as the value of the CouplingOrders option of the MoGRe$Renormalize method (see section B.1). The corresponding field content is subsequently extracted and the relation between the bare and renormalised quantities are derived following eq. (2.21). Matrix renormalisation is by default considered for what concern fields lying under the same representation of the gauge and Poincaré groups. An expansion over all flavour indices is then performed and the field mixing restrictions passed as the FlavorMixing option of the MoGRe$Renormalize method (see section B.1) are finally enforced.
The wave-function renormalisation constants associated with the left-handed and right-handed chiralities of a Majorana fermion being equal, the code simplifies the resulting expression by mapping the right-handed one onto the left-handed one, This has the advantage to prevent NLOCT from calculating twice the same quantity and to yield more compact expressions for the counterterms.

B.6.3 Parameter renormalisation
The renormalisation of the model parameters is accounted for in three steps. External parameters (including internal parameters that have been made external with FR$LoopSwitches, see section B.6.1), particle masses and internal parameters (including external parameters that have been made internal with FR$LoopSwitches) are handled separately. All parameters and masses are directly renormalised according to eq. (2.23), after removing all doubly-declared parameters defined through the FR$RmDblExt variable (see section B.6.1). The code additionally takes care of deriving all relations connecting the renormalisation constants of the internal parameters to those of other (external or internal) parameters. Those relations are truncated at the one-loop level. For instance, the renormalisation constant associated with the third generation trilinear coupling strength of eq. (2.17), δ(T u ) 33 , can be written in terms of the top-quark mass and (A u ) 33 renormalisation constants, The obtained set of relations do not include any dependence on the renormalisation constants that would be associated with unrenormalised parameters (see section B.2) and vanishing parameters have been removed (see section B.3).
As a last step, the code applies all renormalisation conditions that have been provided by the user via the AddRenormalizationCondition method. The relations between bare parameters and the corresponding renormalised ones are modified so that the dependent renormalisation constants are replaced by their functional form. This also concerns the rules relating the renormalisation constant of an internal parameter to other parameters and their renormalisation constants.

B.6.4 Renormalisation of the Lagrangian
All the previously derived parameter and field redefinitions are finally applied to the Lagrangian. A truncation at the one-loop level is performed by the code, so that each Lagrangian term is at most linear in the renormalisation constants. The Lagrangian is then formatted so that NLOCT can be called to derive the UV and R 2 counterterms, following the techniques detailed in ref. [59]. As a consequence, while MoGRe lifts some of the limitations inherent to a joint use of FeynRules and NLOCT, it naturally inherits all limitations that are strictly bounded to NLOCT. For instance, couplings independent of the particle masses are renormalised in the MS scheme and α S has to be renormalised either in the MS scheme or as described in section 2.2.5. This will be addressed in future work.

C Decoupling of heavy SUSY particles in MG5 aMC
Arbitrary SUSY mass spectra often feature SUSY particles with masses much larger than the collision energy scale, and which are therefore expected to decouple in the corresponding crosssection computations.
In the tree-level matrix elements, this decoupling property applies individually to each treelevel Feynman diagram featuring a propagator of the decoupling heavy SUSY particle. In the loop matrix elements, however, this decoupling is realised in a more complicated way involving cancellations among several loop diagrams featuring the decoupling particle(s) running in the loop. These cancellations become more severe as one approaches the decoupling limit, and numerical inaccuracies in the loop computations will eventually spoil them, yielding incorrect predictions.
The solution to this problem simply amounts to completely removing the problematic heavy modes from the process definition, effectively enforcing the decoupling property of the resulting matrix elements. It is however rather impractical having to settle for different process definitions depending on the masses of the heaviest particles in the spectra, and it is therefore desirable to determine quantitatively when this explicit removal procedure really becomes mandatory from a numerical standpoint. To this end, we consider the one-loop matrix element for the pair production of gluinos from initial-state gluons (i.e. gg →gg), for the reason that this particular process features the worst numerical behaviour in the decoupling limit amongst all 2 → 2 one-loop SUSY matrix elements. We fix the kinematic configuration to √ s = 4mg = 2 TeV, and report in figure 15 the numerical evaluation of this matrix element as a function of the mass of a decoupling down squark. We observe that the numerical evaluation starts to depart away from the decoupling limit for downsquark masses around 500 TeV in double precision and 5000 TeV in quadruple precision. However, these numerical instabilities will start to significantly impact the accuracy of the final result only when reaching even larger masses of about 10 3 TeV and 10 4 TeV respectively. In view of these results, we conclude that a conservative recommendation is to manually enforce the decoupling of certain heavy particles by removing them explicitly from the process definition when their masses are larger than about a thousand times the characteristic energy probed by the observable. We stress however that in any case we tested that the associated numerical instabilities are correctly detected by MadLoop and will therefore adequately be reported to the integrator, which will in turn set the corresponding weight to zero if the accuracy is too poor. If such exceptional configurations g g → g g, complete ME Decoupling limit Ninja DP Ninja QP COLLIER The other lines denote a numerical evaluation using various one-loop reduction algorithms and implementations of scalar master integrals: Collier (black) and Ninja in double (DP, green) and in quadruple (QP, red) precision.
occur too frequently, a clear warning is issued, hence preventing the user from inadvertently using incorrect results.