Squark production in R-symmetric SUSY with Dirac gluinos: NLO corrections

R-symmetry leads to a distinct realisation of SUSY with a significantly modified coloured sector featuring a Dirac gluino and a scalar colour octet (sgluon). We present the impact of R-symmetry on squark production at the 13 TeV LHC. We study the total cross sections and their NLO corrections from all strongly interacting states, their dependence on the Dirac gluino mass and sgluon mass as well as their systematics for selected benchmark points. We find that tree-level cross sections in the R-symmetric model are reduced compared to the MSSM but the NLO K-factors are generally larger in the order of ten to twenty per cent. In the course of this work we derive the required DREG $\to$ DRED transition counterterms and necessary on-shell renormalisation constants. The real corrections are treated using FKS subtraction, with results cross checked against an independent calculation employing the two cut phase space slicing method.


Introduction
The Minimal Supersymmetric Standard Model (MSSM) is one of the most studied extensions of the SM. Often, an unbroken R-parity is assumed which implies that supersymmetric particles can only be produced in pairs. In this paper we consider a distinct realisation of supersymmetry (SUSY) based on an unbroken, continuous R-symmetry. The basic feature of R-symmetry is that particles and superpartners have different R-charges where the differences are unambiguously prescribed by the SUSY algebra. For definiteness we focus on the Minimal R-symmetric Supersymmetric Standard Model (MRSSM) [1] but our discussion will apply also to more general R-symmetric models.
From the phenomenological point of view, the MRSSM is an appealing model, with immediate restrictions following from R-symmetry. The model contains Dirac instead of Majorana gauginos and adjoint scalar superpartners for all gauge fields. This also leads to significant changes in the Higgs sector due to the presence of additional singlet and triplet scalars. The µ-and A-terms of the MSSM are forbidden; mixing between left and right handed squarks or sleptons is forbidden, and large contributions to flavour and CP violating observables are suppressed [1,2].
In a recent series of papers, the electroweak sector has been investigated and it has been shown that the MRSSM can accommodate the experimentally measured W and Higgs boson mass as well as the electroweak precision variables and is compatible with direct detection searches for dark matter and LHC searches of electroweak particles [3,4]. For related studies, see [5][6][7].
In this paper we focus on the strongly interacting sector of R-symmetric SUSY. It is characterised by Dirac gluinos, scalar gluons and left and right handed squarks which are mass eigenstates and have opposite R-charges. The phenomenology of this SUSY QCD sector of the MRSSM has already been studied at the tree level [8,9], where it was shown that conservation of R-charge is responsible for a rather drastic suppression of the inclusive squark production cross section leading to lower bounds on squark masses. This suppression is even more amplified when the gluino mass is increased. 1 Corrections to this observable at the next-to-leading order (NLO) have only been approximately estimated in ref. [11] by using global MSSM instead of MRSSM K-factors. It is the purpose of this paper to describe an exact NLO SUSY-QCD calculation in the MRSSM at the example of squark-squark (qq) and squark-antisquark (qq † ) production at the LHC. We will expose and explain differences to the analogous calculation in the MSSM.
The paper is structured as follows. The next section describes the strongly interacting sector of the MRSSM. In section 3 we evaluate the leading order cross sections for the production of colour charged MRSSM particles. Section 4 describes the evaluation and renormalisation of the virtual amplitudes. Our results have been obtained with two independent codes which use different regularisation schemes [12]. We provide a list of counterterms including the transition counterterm from dimensional regularisation to dimensional reduction. Section 5 is devoted to the description of real corrections. We present two alternative ways of dealing with infrared singularities used in this work: the two cut phase space slicing method and the FKS subtraction. We also discuss the treatment of the on-shell resonances. In section 6 we proceed with a phenomenological analysis. We give a detailed comparison of K-factors in the MSSM and the MRSSM, explaining the physical origins of their differences. Then we give an overview of the NLO corrected cross sections in different regions of parameter space, including an uncertainty discussion. Section 7 contains our conclusions, and the appendix collects a list of Feynman rules, implementation details, and numerical results for verification purposes.

Details of the model
The field content of the MRSSM is enlarged in comparison to the one of the MSSM. The necessity of this arises from R-symmetry. As can be seen from table 1 the gluinog L , as superpartner of the gluon, has R-charge +1. In order to account for a non-vanishing Dirac gluino mass it needs to be partnered with a new fieldg R with R-charge −1. The Dirac nature of gluinos manifests itself in an N = 2 supersymmetric gauge sector including scalar gluons O, called sgluons, which transform in the same representation as the gluon under gauge transformations. The corresponding Lagrangian of the strongly interacting part of superfield boson fermion left-handed (s)quark the MRSSM, including one massless quark of arbitrary flavour reads Note that the vector superfield of the gluon V s in the first line of eq. (2.1) transforms for each term in the appropriate representation of SU (3) C , i.e. the fundamental, the antifundamental and the adjoint one. The soft breaking Lagrangian accounts for squark, gaugino and sgluon masses. These mass terms arise from a hidden sector spurion. For gauginos the D-type spurion is given by W α = θ α D and mediates supersymmetry breaking at the mediation scale M : dθ 2 W α M W α s O. After integrating out the spurion one obtains [3,13] L soft = − m 2 where D a is the usual auxiliary field in the SU (3) sector. The complex sgluon state has to be decomposed into two real fields: 3) The mass of the CP-even scalar O s receives an additional contribution from the gluino Dirac mass whereas the mass of the pseudoscalar O p is solely given by the soft breaking parameter: The  3 Squark and gluino production at the leading order As a starting point for the NLO SUSY-QCD calculation, we recall the tree-level production of squarks and gluinos in the MRSSM 2 and point out differences to the familiar MSSM. A detailed study including comparisons of Dirac, Majorana and hybrid gluinos can be found in ref. [16]. As in the MSSM there are six partonic channels contributing to squark and gluino production. Three of which for squark-antisquark and squark-squark production and three for gluino-antigluino and squark-(anti)gluino production q i q i →gg, gg →gg, The inclusion of charge conjugated processes is understood if they exist. The indices denote quark flavours. The corresponding Feynman diagrams are shown in figure 1. In the following Tree-level diagrams for squark and gluino production in the MRSSM. For simplicity, only one quark flavour is shown. In the third and last line, also a charge conjugated process exists.
we give analytic formulae for the partonic cross sections for several channels in the MRSSM and explain if/how they differ to their respective analogue in the MSSM. We sum over the n f − 1 first (s)quark flavours (where n f ≡ 6). 3 The gauge coupling g s from the qqg-vertex is identical to its supersymmetric analogueĝ s from the qqg-vertex. Taking all squarks as mass degenerate, we obtain the leading order partonic cross sectionŝ whereŝ is the squared partonic centre-of-mass energy and the following abbreviations of ref. [17] are used In comparison to the MSSM there are two main differences. Firstly, an overall restriction stemming from an unbroken R-symmetry is that the final state particles' R-charges must sum up to zero. Hence in the MRSSM, only diagrams for qq →q Lq † L , qq →q Rq † R and qq →q LqR can be drawn as shown in figure 1. Processes with a squark-antisquark (squarksquark) pair of different (same) "chiralities" in the final state are forbidden by R-symmetry.
For the allowed channels in the MRSSM the chirality projectors lead to the replacement for the gluino propagator in the first and third line of figure 1. On the level of the cross section, this manifests in the substitution 3.11) and the vanishing of the term proportional to δ ij forqq production in ref. [17], compared to the MSSM expressions. Expandingσ B (q i q j →qq) appropriately for large values of mg, shows that the leading term in the MSSM is proportional to m −2 g , whereas in the MRSSM it is m −4 g , as expected from eq. (3.10) 4 . A second feature of the MRSSM is that gluino and antigluino are no longer indistinguishable particles, which manifests in twice as much degrees of freedom available for the gluino in the final state. In contrast to the above mentioned feature, this characteristic induces an increase of some cross sections. It strikes very clearly in the cross section of gg →gg, which is doubled in comparison to the MSSM. For the process qq →gg, both features appear: The first line of eq. (3.6) is twice as much as its MSSM analogue in ref. [17] but the following two are not. This is because those originate from t-and u-channel diagrams where only one instead of two squark "chiralities" occur. In addition, the MRSSM result misses the t-and u-channel interference. In the last channel, i.e. qg →qg, both features are present and cancel exactly. On the one hand, R-charge allows only the production of "right-handed" squarks in association with the gluino. On the other hand, there is a distinct antigluino which can be produced with a "left-handed" squark. As forqq production, the charge conjugated process exists.

Figure 2:
The production cross sections for squarks and gluinos at the LHC with √ S = 13 TeV in the MSSM (left) and the MRSSM (right). It is summed over five flavours, all possible squark "chiralities" and for squark-squark and squark-gluino production also the charge conjugated processes are taken into account. The top row contains results in function of the gluino mass (with squark masses fixed to 1.5 TeV), the bottom one in function of squark masses with the gluino mass fixed to 2 TeV. The PDF set used is MMHT2014LO [18]. As renormalisation and factorisation scale µ R = µ F = m 1 +m 2 2 has been chosen, where m i are the final state's particle masses.
In conclusion, we obtain both suppressing as well as amplifying effects in the MRSSM. In contrast to the Dirac gluino, the presence of sgluons has no effect on the discussed treelevel processes. Figure 2 shows the convolution of the partonic cross sections with parton distribution functions (PDFs) in the MSSM and the MRSSM for the 13 TeV LHC. For details regarding the parameters used, see the caption of figure 2. We sum over all "chiralities" and five flavours as well as charge conjugated processes, when distinct. Due to the prohibition of the above mentioned "chirality" states in the MRSSM, we obtain a suppression in thẽ qq † andqq production (by a factor of 1.3 and 13.2 at BMP2, respectively), which increases with the gluino mass. Furthermore, an amplification of gluino production, due to its Dirac nature, is visible. In combination with experimental results, the absence of detected gluinos translates into a larger exclusion limit of the gluino mass than in SUSY-QCD [8]. We will thus focus on the region of parameter space with a rather large gluino mass when discussing our NLO results.

Virtual corrections
The present and the subsequent sections describe the calculation of the NLO SUSY-QCD corrections to squark production processes in the MRSSM. 5 The inclusion of gluino final states we postpone for future work, since they are not as important in the motivated scenario where squarks are significantly lighter than gluinos. Thus, on hadron level we consider squark-squark production pp →qq and squark-antisquark production pp →qq † . On the partonic level we have to consider NLO corrections to qq →qq for the first process and NLO corrections to qq →qq † and gg →qq † for the second process.
In the present section we focus on the virtual corrections. They involve ultraviolet and soft and collinear infrared divergences. These divergences are removed after renormalisation and combination with real corrections.
In intermediate steps, the divergences need to be regularised. However, computations of SUSY-QCD corrections suffer from the fact that no regularisation is at the same time directly compatible with SUSY and the standard definition of PDFs, see also ref. [12] and ref. [19] for a recent review. Hence we have done the calculation in two complementary ways.
• The first calculation uses dimensional regularisation (version HV in the notation of ref. [12]) and Passarino-Veltman reduction of one-loop integrals. In the HV scheme all four-dimensional quantities are treated in D = 4 − 2 dimensions, except that particles outside loops are kept unregularised. In this calculation SUSY is broken in intermediate steps and SUSY-restoring counterterms must be added.
• The second calculation uses dimensional reduction/the four-dimensional helicity scheme (version FDH in the notation of ref. [12]), helicity methods and integrand reduction techniques. In the FDH scheme, only space-time and momenta are treated in D dimensions, while gluons are kept quasi-four-dimensional. Like in HV, only particles inside loops are regularised. In this calculation a finite shift is needed in the strong coupling renormalisation, and transition rules have to be applied to convert infrared divergent amplitudes back to the HV scheme.
The two calculations are in full agreement. In the following we provide the details, first on the general renormalisation scheme, then on the two calculational procedures. We need to consider the one-loop amplitudes for the three partonic processes qq →qq † , gg →qq † , and qq →qq. Sample Feynman diagrams for the calculation of renormalisation constants are shown in figure 3. Diagrams contributing to the production of (anti)squarks are shown in figure 4. Diagrams of type a) are exactly the same in the MRSSM and in the conventional MSSM and give equal results. Diagrams of type b) can be drawn in the MRSSM but unlike in the MSSM do not contribute. The result is proportional to the Majorana mass and leads therefore to R-charge violation in the MRSSM. Diagrams of type Figure 4: Examples of Feynman diagrams relevant for the NLO corrections toqq andqq † production. Categories as in figure 3.

One-loop diagrams and renormalisation
c) exist in the MRSSM but are absent in the MSSM because of the appearance of sgluons, additional strongly interacting particles.
The ultraviolet renormalisation requires the introduction of coupling, mass, and field renormalisation. We define all SUSY masses in the on-shell renormalisation scheme to have physical masses as input. The quark masses are set to zero, except for the top quark mass. For simplicity we take all squarks as mass degenerate. We define the field renormalisation in the on-shell scheme, which leads to a correctly normalised S-matrix but also to infrared divergent renormalisation constants. The strong coupling is renormalised in a decoupling scheme as described below.
With these definitions, standard methods lead to the following results for the quark/squark field and squark mass renormalisation constants, computed with one-loop SUSY-QCD cor-rections: We use standard Passarino-Veltman integrals and their derivatives, see e.g. ref. [20] for the definition. The renormalisation constants are computed in both the HV and FDH scheme; we denote terms only appearing in HV by the subscript HV. We do not have to distinguish left-and right-handed squarks since we do not take electroweak corrections into account.
Though not strictly necessary for our calculations, we also quote the result for the gluino field renormalisation. Here the left-and right-handed parts renormalise differently, since they are part of different superfields, i.e. the right-handed part of the gluino does not couple to (s)quarks. This is reflected by a gluino self-energy which is not left-right symmetric and has the following basic structure: where the dots stand for contributions not stemming from (s)quarks. The results read The corresponding gluino mass renormalisation constant is Note that the contribution stemming from the quark-squark loop is halved when compared to the MSSM result, due to the non-coupling right-handed gluino. As a consequence, the left-handed Dirac gluino allows only for a "left-handed" squark or a "right-handed" antisquark in the loop. Finally, the field renormalisation constant of the gluon is given by The renormalisation of the strong coupling has to be treated in a special way in order to make use of the experimental determination of α s in the SM 5-flavour scheme and for compatibility with available PDF sets. The renormalisation has to be matched to the SM MS 5-flavour scheme, and contributions from heavy particles to the renormalisation of g s have to be subtracted at zero momentum. In practice, we have separated the loop-diagrams used in the determination of δg s into contributions from light and heavy particles. From loop-diagrams involving solely light particles we have kept only the part corresponding to the MS-scheme, whereas for diagrams involving heavy particles we have adopted zeromomentum-subtraction. This leads to with the typical UV-divergent constant ∆ defined as in ref. [20]. As a result, the renormalisation constant δg s contains additional µ-dependent terms, where µ is the MS renormalisation scale. These have the effect of decoupling heavy particles from the running of g s , which is then given by the SM 5-flavour β-function. As a side remark, we mention that if δg s was defined in pure MS-scheme, the β-function would vanish at one-loop-level.

Method 1: HV and Passarino-Veltman reduction
In our first method we have performed the calculation in HV regularisation, i.e. usual dimensional regularisation for internal particles, while particles outside of loops are kept unregularised.
It is well known that dimensional regularisation breaks SUSY due to a mismatch between degrees of freedom of the gluon and the gluino in D = 4 − 2 = 4 dimensions. Since dimensional reduction [21,22] is known to preserve SUSY at the one-loop level (for reviews of checks see e.g. [23][24][25]), the required SUSY-restoring counterterms can simply be obtained from comparing renormalisation constants in dimensional regularisation and dimensional reduction. Appropriate transition counterterms are known for physical parameters in generic SUSY models at one-loop level [26], for the full MSSM at one-loop level [27], and for SUSY-QCD at two-loop level [28]. For our calculation, only one such transition counterterm is needed: the one for the squark-quark-gluino vertex. Denoting the renormalisation constant for this vertex by δĝ s , we find that it has to satisfy where δg s is given by eq. (4.10). The result for δg restore s is the same as in the MSSM, since SUSY-breaking is only associated with the gluon, which has the same couplings in the MSSM and MRSSM.
The implementation of this calculation has been done with the help of several Mathematica packages. The generation and processing of amplitudes has been performed by FeynArts [29] and FormCalc [30,31]. The model file for the MRSSM containing the treelevel vertices was generated by SARAH [32][33][34][35]; the one-loop counterterms Feynman rules were included by hand into the model file. The output has been passed on to a C++ program which performs the evaluation of loop integrals using LoopTools [36] and does the integration using the CUBA library [37].

Method 2: FDH and integrand reduction approach
In our second calculation we employ the FDH regularisation scheme. We follow the notation of ref. [12], so FDH is the same as standard dimensional reduction, but "regular" particles outside of loops are kept unregularised. The FDH scheme is advantageous not only because it preserves SUSY at one-loop level, but also because it allows the use of powerful and efficient helicity methods for the evaluation of loop amplitudes.
In the past, dimensional reduction was often not applied to SUSY-QCD calculations because of an issue with QCD factorisation discussed in refs. [38,39]. In the meantime it has been understood that factorisation behaves differently depending on whether FDH as defined above, or whether DRED, where also particles outside loops are regularised, is employed [12,40]. 6 In case of FDH, factorisation holds as expected, however the infrared anomalous dimensions are different from the HV scheme and the transition rules found in [41,42] apply. In case of DRED, factorisation is complicated by the appearance of external -scalars with separate couplings and anomalous dimensions. The understanding of the infrared behavior of all these schemes and their transition rules has been extended to the multi-loop level in refs. [43][44][45][46][47]. The upshot of these results is that FDH can be used to regularise ultraviolet and infrared divergences for SUSY-QCD processes. However in order to combine the results with real corrections evaluated in HV and to convolute with usual PDFs, we need to convert the amplitudes from FDH to HV. The appropriate transition rules for the squared amplitudes are, in the notation of ref. [12] |M = C A /6 and the sum running over all external partons.
An alternative to Passarino-Veltman reduction for NLO calculations with Monte Carlo methods is the usage of helicity methods and an integrand reduction approach, see e.g. ref. [48] for a review. For these methods, several computer codes already exist. They contain a general implementation of reduction methods and require the input of model dependent information. We summarise these parts first and then describe the implementation used for our calculation.
Needed for our calculation are the relevant Feynman rules of the MRSSM and the counterterms. The necessary renormalisation constants of masses and wave functions are given in eqs. (4.1) to (4.9). The renormalisation of the coupling constant g s given in eq. (4.10) includes the finite term containing the expression 1 − 1| HV which marks the wellknown transition between DR and MS scheme needed for FDH. Including this transition rule allows us to combine matrix elements calculated in FDH with PDFs given in the MS scheme. The remainder of the calculation is done automatically as described in the following.
For our second calculational method, we use GoSam 2 [49,50] 7 to calculate the virtual corrections using helicity methods and integrand reduction methods. The information concerning the MRSSM is passed to GoSam using the UFO interface [60]. For this, the additional strongly interacting particle content of the MRSSM was added to the SM implementation in FeynRules [61]. Numerous checks were performed to verify that the SARAH and FeynRules model files give the same tree-level results for the relevant processes.
The UV renormalisation of the MRSSM has been added by hand to GoSam. 8 The counterterms are implemented using OneLOop for the loop functions and added to the matrix.f90 template in the GoSam interface. The exact counterterm structure has to be fixed once after generating the considered process with GoSam.
Comparing the computational speed between both Methods both implementations lead to similar running times for the considered processes with the timings being within an order of magnitude of each other. The difference are in general not relevant when combining with the real corrections of section 5 which take up the majority of the computing time for the full calculation.

Real corrections
The singularities remaining in the virtual matrix elements, which are of infrared (IR) origin, cancel after combination with the real corrections. The singularities left in the real emission processes are then removed by mass factorisation [64,65].
A multitude of methods that implement the above cancellation on a numerical basis have been developed. In the following two subsections, we will discuss the two approaches used in this work: the two cut phase space slicing method (TCPSS) [66] and the Frixione-Kunszt-Signer (FKS) subtraction [67,68].

Method 1: Two cut phase space slicing (TCPSS)
A review of TCPSS can be found in ref. [66]. The method has been implement and tested by one of us in the context of the (S)QCD particle production on the process of sgluon pair production in refs. [69,70]. In this method the three-body real emission phase space is decomposed with respect to the additional parton radiation into: soft (S), hard collinear (C) and the hard non-collinear (HC) regions by introducing two parameters δ s and δ c , taken to be numerically small. This can be schematically written as where in the last step the hard part H is split into collinear (HC) and non-collinear (HC) pieces. Within the soft and collinear approximations the divergences are then dimensionally regularised and extracted analytically. We will now present technical details needed to carry out this calculation.

Soft emissions
The soft phase space S is defined by the condition that the energy of an outgoing gluon E 5 , in the rest frame of colliding partons, fulfils where δ s 1. We label incoming parton momenta as p 1 and p 2 , such thatŝ = (p 1 + p 2 ) 2 . In the soft limit and 4 − 2 IR dimensions the 2 → 3 amplitude can be written as where p 5 is the gluon momentum and is the non-abelian eikonal current (whose sum extents over all partons except for the finalstate gluon), which is colour-connected to the 2 → 2 process amplitude M 2 through the colour operator T of the particle f . Only the first term in eq. (5.3) contributes to the singular part of |M 3 | 2 , as the interference term is regular in the limit δ s → 0.
Similarly, in the soft limit the three-body phase space factorises, with a phase space measure where θ 1/2 describe the direction of the gluon emission in the rest frame of colliding partons. The gluon phase space is integrated over the solid angle, and E 5 is as given in eq. (5.2). For convenience, the necessary expressions for the D-dimensional angular integrals have been gathered in appendix D.2 of ref. [69].
The expressions for the HV regularised real emission matrix elements have been generated using FeynArts and FormCalc with the model file described in section 4.2. Expanding in terms of IR we find the double and single poles, and a finite part. The double-pole terms agree with the well-known minimal structure, being proportional to the four-dimensional Born cross sections given in eqs. (3.3), (3.4) and (3.5), which can be found in ref. [12]: For numerical verification we show the cancellation of these terms and the virtual matrix elements for a single phase space point in appendix C. The single-pole coefficient is not cancelled completely between virtual and soft contribution. The remaining terms have the form where N is the number of colours and n f −1 is the number of massless quark flavours. These uncancelled terms come from the phase space region where the gluon is collinear with an incoming parton, but its energy is non-zero. They do not cancel out with the virtual contributions as they have different kinematics. As discussed in the next subsection, these are the terms that can be absorbed by a redefinition of the PDF at NLO.

Collinear emissions
In the collinear limit, defined by the condition [71] where δ c 1 and i = 1, 2, the real-emission cross section factorises at the level of the absolute squared matrix element. Contrary to the definition of the collinear region in ref. [66], eq. (5.11) decouples soft and collinear regions.
The double differential hadronic hard-collinear cross section is given by Note that there are two possible ways in which qq in the initial state can be obtained from the proton-proton system. The P ik (z, ) are the D-dimensional unregulated Altarelli-Parisi (AP) splitting kernels [72] The δ ik in the integration boundaries of eq. (5.12) ensures that the integral is taken up to z = 1 − δ s for kernels which are singular as z → 1 (P qq and P gg ). The remaining P gq kernel is obtained by the replacement P gq (z, ) = P qq (1 − z, ).
The Bjorken variable in f k/p is rescaled so that in the Born configurationσ B is taken atŝ = x 1 x 2 S. 9 Collinear singularities will cancel out with the renormalised PDFs. The first-order correction to i-th flavour PDF in the MS prescription is given by where P + ij (z) are the '+' regulated AP splitting kernels where the associated '+' prescription is defined as For quark radiation, the integral will be taken up to 1 as there is no soft singularity.
For partonic processes which have a soft singularity there is a mismatch in the z integration boundary between eq. (5.12) and eq. (5.16). To account for that, we write eq. (5.19) as 20) where the term in the last line is of (at least) O(δ s ) and can be neglected. Eq. (5.16) then reads where now the unregularised, four dimensional AP splitting kernels appear and the softcollinear factors A sc 1 for the splittings with a soft gluon (g) are given by The first term is just the Born partonic cross section convolved with the scale-dependent PDFs. The terms A sc cancel out with eqs. (5.9) and (5.10). Combining now eqs. (5.12) and (5.24) gives, together with the LO cross section, a final result for the hard-collinear In appendix C we verify the cancellation of single poles for all partonic processes considered in this work at one phase space point.

Method 2: FKS subtraction
As an alternative approach for the calculation of the real corrections we make use of the FKS subtraction scheme. Using this subtraction method, suitable subtraction terms for individual soft and collinear singularities in the squared matrix elements are constructed which allow for a convergent numerical integration over the phase space.
This scheme is implemented as automatised method MadFKS [73] in the Monte Carlo program MadGraph5_aMC@NLO [74,75] (MG5aMC@NLO). The application of the FKS scheme is only dependent on QCD specifics and has been used for many different previous calculations. The program is well tested by applications to the SM [76][77][78] and to a multitude of BSM models including models with SUSY [79], extra dimensions [80], leptoquarks [81] and dark matter candidates [82,83] as well as the Two-Higgs-Doublet Model [84,85] and the Georgi-Machacek model [86].
An interface between MG5aMC@NLO and GoSAM for SM calculations is already implemented [87] as specified with BLHA1 [88]. To allow for calculations in BSM models we extend this to the BLHA2 [89] conventions. The appropriate changes are summarised in appendix B. 10 The correctness of the implementation has been tested thoroughly by ensuring that all appearing divergences cancel between the virtual and real part for all subprocesses in various regions of phase space.

Comparison of TCPSS and FKS subtraction
In TCPSS, the dependence on the (unphysical) regulators δ s and δ c should vanish after adding the hard non-collinear emissions. We verified extensively that this is indeed the case. Examples for the cancellation of the cut parameter dependence are shown in figure 5 using MMHT2014nlo68cl PDFs [18] interfaced via LHAPDF6 [91]. In figure 5a we plot the hard non-collinear (blue) and (summed) virtual, soft and collinear parts (red) in function of the phase space slicing parameter δ s for the uu →ũ LũR channel for BMP2 and δ c = 10 −6 . In the bottom subplot, the sum is then compared with the FKS result ( with a shaded 10 For another recent application of GoSam with MadFKS see ref. [90].  Figure 5b shows the same for the uū →ũ Lũ † L process and BMP1 as a function of the δ c . The comparison with the FKS subtraction reveals that a relative accuracy of 10 −4 requires the cut parameters δ s = 10 −5 and δ c = 10 −6 , which is what we use in our numerical analyses. As the TCPSS method relies on a partial cancellation of two large contributions at the level of an integrated cross sections, it is inherently slower than local subtraction schemes like FKS. For the case at hand, this means an order of magnitude slow down for the same final precision. This directly translates into an order of magnitude speed difference between the standalone C++ code and the MG5aMC@NLO as for the standalone calculation most time is spend evaluating the 2 → 3 hard non-collinear part. is not exclusive to BSM models as a similar issue appears in the case of the real corrections to the SM tW production (see e.g. ref. [92]), which contain resonant contributions from the top pair production.

Treatment of resonances in real emission diagrams using diagram removal
A popular way of dealing with this problem is the diagram removal technique. In this approach one either completely discards all resonant diagrams at the level of the amplitude or their square at the level of the squared matrix element. In this work we employ the first of those solutions. 12 This version of the DR approach is also the default way of dealing with resonant divergences when using MadFKS in MG5aMC@NLO.
The removal of a subset of Feynman diagrams in general violates gauge invariance and care has to be taken to ensure that the effect is numerically small. For the case of the MSSM this has been studied in depth in refs. [94,95] (see also refs. [96,97] for other implementations). Also, it requires a careful choice of a gauge not to spoil the factorisation of collinear singularities.
To understand this last point, consider the gluon splitting g(p 1 ) →q * (p 1 − p 5 )q(p 5 ) connected to some bigger amplitude throughq. We use the * -symbol to emphasise thatq is in general not on-shell. The final state quark momentum p 5 can be parametrised through Sudakov decomposition as where n µ is a reference null vector, 1 − z is the fraction of gluon energy carried by the quark, and p 5,⊥ · p 1 = p 5,⊥ · n = 0. In the limit of p 5,⊥ → 0 the vectors p 1 and p 5 become spatially parallel and diagrams with this splitting develop a (p 1 − p 5 ) −2 = (1 − z)/p 2 5,⊥ singularity. Due to the helicity conservation, physical gluons cannot decay to a pair of on-shell quarks. Therefore, in the physical gauge the matrix element will always have additional power of p 5,⊥ in the numerator. As the one particle phase space for radiation of q is given by dΦ 1 ∼ dp 2 5,⊥ , the collinear singularities of the interference terms are integrable since those terms scale as dp 2 5,⊥ /p 5,⊥ . This is no longer true in an unphysical gauge, where longitudinally polarised gluons might appear. For the full amplitude with a ug initial state, 12 An alternative approach would be to employ the diagram subtraction method (sometimes also referred to as the Prospino scheme [93]). We have checked that switching to this choice changes the total cross sections at a percent level. We therefore postpone the detailed studies of this method to the publication documenting the RSymSQCD code. which is gauge invariant, longitudinal states decouple through a Ward identity as no triple gluon vertices appear. This is no longer true after the removal of resonant diagrams, though.
We therefore calculate the real emission matrix elements in the light cone gauge. A convenient choice of the gauge-fixing vector η is the momentum of the other incoming parton p 2 . This choice allows to avoid spurious divergences present in the polarisation sum as p 1 · p 2 =ŝ/2. 13 To study the numerical impact of this gauge choice we consider two families of gauge vectors (assuming the momenta p 1,2 oriented in the ±z-direction) For δ = 0 the η − choice is equivalent to the choice of η = p 2 , while the η + choice is singular as p 1 · η + = 0. In the limit of δ 1 results of both gauge choice converge as η + η − . Both of these features are clearly seen in figure 7, where we plot the cross sections for gq →ũ Lũ † L q with q = u,ū (7a) and gu →ũ LũRū (7b) in function of the δ parameter for BMP2. For gu →ũ LũRū the gauge dependence is enhanced since resonant diagrams were giving a substantial contribution to this amplitude.
We point out, that keeping the choice of the non-singular gauge η − , the DR subtracted processes give per mille level contribution to the total cross section. Hence, for all intents and purposes, the gauge dependence is not a problem and for studies in the following sections we use the customary choice of η = p 2 .

Results
When the virtual corrections of section 4 are combined with the real contributions from section 5 all IR divergences cancel as expected and we are able to calculate the NLO SUSY-QCD cross section forqq andqq † production in the MRSSM. In this section we first point out and explain the consequences of R-symmetry on NLO cross sections in supersymmetric QCD before exhibiting the results of our calculations forqq andqq † production. We conclude with remarks on the scale dependence and uncertainties of the computed observables. 14 A useful quantity which helps to understand these different aspects is the K-factor. We define a K-factor K as ratio of NLO over LO total cross section following ref. [17]: where σ NLO is calculate using NLO PDFs and σ LO using NLO or LO PDFs, depending on the argument. If no argument is given we use LO PDF sets as default case for σ LO . The K-factor depends on the process, masses of the fields in the model, and the centre-of-mass energy, as well as the choice of renormalisation and factorisation scale. In the following, we set √ S = 13 TeV and µ R = µ F = m 1 +m 2

2
, where m i are the final state's particle masses. For a consistent definition we identify the values of the Dirac gluino mass in the MRSSM with the Majorana gluino mass in the MSSM, the values of the squark masses in both models, as well as neglecting left-right squark mixing in the MSSM. In the MRSSM the K-factor has an additional dependency on the sgluon mass parameter. To fix all remaining SM parameters, we set the top quark mass to m t = 172 GeV and take the strong coupling constant α MS s (m Z ) from the used PDF set.

Effects of R-symmetry
The results of R-symmetry at NLO may be pinned down to two features: the presence of a Dirac instead of Majorana gluino and the existence of sgluons. The conservation of R-charge, already discussed at LO in section 3, also needs to be commented on when comparing K-factors to the MSSM. The effect of R-symmetry is however only present in the virtual corrections, i.e. the real corrections do not differ in the MSSM and MRSSM. The effects discussed in the following comprise the unique features of MRSSM neglected in the study of ref. [11] by only taking MSSM NLO K-factors into account.

Dirac nature of the gluino
In the presentation of renormalisation constants in section 4.1, we already saw ramifications of the Dirac gluino whose components couple differently. Now, we study the effect of replacing the Majorana with a Dirac gluino in a physical process.
To this end, consider diagram 1a) and 1b) (referring to the first diagram of category a) and b), respectively) of    Figure 9: Dependence of the K-factor on the sgluon mass forqq andqq † production in the MRSSM.
latter is only non-zero in the MSSM, since the gluino undergoes a chirality flip to its noncoupling right-handed component. Taking only objects with spinor indices into account these diagrams evaluate to M 1a) ∝v(p 2 )P R (/ t + mg)P L/ qP R (/ t + mg)P L u(p 1 ) =v(p 2 )/ t / q/ tP L u(p 1 ), Here q µ is the quark's loop momentum and mg is the Majorana gluino mass of the MSSM, which is zero in the MRSSM. Hence diagram 1a) is the same in both models, while diagram 1b) vanishes in the MRSSM. Alternatively, diagram 1b) can also be understood to be zero in the MRSSM as R-charge is not conserved at all vertices of the diagram as the R-charge of the gluino flows in the opposite direction than the one of the squark appearing in its self energy insertion. On similar grounds we are able to understand why diagram 2b) in figure 4 does not contribute in the MRSSM. Its matrix element has the following proportionality in the MSSM and is therefore zero in the MRSSM. Note that its analogue, depicted as diagram 2a) in figure 4 is proportional to the u-quark mass and therefore zero in both models. Diagram 3b) of figure 4 appears only forqq † production and is proportional to the Majorana mass of the MSSM squared for the same reasoning as discussed before. The number of diagrams affected by the difference between Dirac and Majorana mass is only a subset of all contributing virtual graphs. Therefore, the effect is most apparent when we use the dependency on the gluino mass and make it large compared to all other appearing scales. This can be seen when comparing the plots of figure 8. On the right plot, where the mass scales are not too different from each other, the MRSSM and MSSM line have a similar dependency on the gluino mass, where the magnitude of the K-factor is affected by the sgluon mass as described below. If, however, the squark mass is reduced, the behaviour changes drastically for large gluino masses. Then, the additional contributions in the MSSM lead to a K-factor rising with gluino mass, instead of falling like in the MRSSM.

Sgluon non-decoupling effects
The sgluon necessarily appears in the MRSSM as the scalar superpartner of the additional component of the Dirac gluino. The sgluon affects the MRSSM prediction in several ways. Obviously, it enters the virtual corrections as new matter content, contributing to the β function of the strong coupling such that it becomes zero at one-loop level. Furthermore, the mass of the pseudo-scalar (scalar) component is determined (mostly) by the SUSY breaking mass parameter m 2 O , which is independent of the gluino mass as can be seen in eq. (2.4). It is therefore possible that the sgluon can be much heavier than the gluino and not be directly observable at the LHC. Still, effects of the sgluon could be experimentally accessible, even if superheavy, as a mass splitting between superpartners in a SUSY multiplet leads to non-decoupling, so called super-oblique, contributions [99].
Super-oblique effects lead to a physical difference between the gauge coupling g s and the gaugino couplingĝ s at loop level. It can be understood in the context of an effective field theory (EFT) where the sgluon is integrated out. In the EFT, the sgluon decouples from the gauge boson and gaugino self-energies (first two diagrams in category c) of figure 3) which changes the corresponding field renormalisation constants and therefore induces a change in the corresponding couplings g s andĝ s . This produces a one-loop difference of in the theory without sgluon. The production ofũ LũR at the LHC is mediated by the gluino alone at tree level. Therefore, the super-oblique contribution from the sgluons can be written as product of the total tree-level cross section and the difference given by eq. (6.5) such that The result forqq † production is more complicated as the gluino only appears in a subset of the diagrams at tree level. Hence, the non-decoupling correction to the cross section affects only a part of the contributions, not allowing for a simple factorisation as before.
The sgluon can influence the prediction of the MRSSM in another way, namely via the sgluon-squark-antisquark vertex induced by the gluino mass parameter, see 9 of the Feynman rules in appendix A. This leads to additional terms in the renormalisation constants, stemming from the third Feynman diagram in category c) of figure 3, as well as additional virtual corrections to the processes, see the last three diagrams of category c) in figure 4. We have checked that these contributions are small compared to the super-oblique corrections.
In figure 8, the effect of super-oblique corrections is visible when comparing lines with different sgluon masses. As expected, their spacing follows the logarithmic dependency described by eq. (6.6) for small gluino masses. When comparing the lines for m O = 3 TeV and m O = 10 TeV at large gluino masses it can be seen that the differences between the K-factors is reduced, as the assumption that sgluons can be integrated out is not valid. Then, also ratios of the sgluon masses to the Mandelstam variables become relevant. Figure 9 shows the effect at very large sgluon masses for bothqq as well asqq † production. For both processes we find the logarithmic enhancement for large sgluon masses, more prominently forqq production as discussed before. Super-oblique corrections lead to a difference in the K-factors of roughly twenty (five) per cent when the sgluon masses changes by two orders of magnitude forqq (qq † ) production.

R-charge forbidden processes
Forqq production we have seen that in the MRSSM only the production of left-and righthanded squarks together is allowed, while the production of squarks with the same "chirality" is forbidden by R-charge conservation. This is relevant when comparing the MRSSM Kfactor for e.g. pp →ũ LũR to the MSSM: should we compare it to the MSSM K-factor for pp →ũ LũR , or to the MSSM K-factor for total squark production, pp →ũ LũR ,ũ LũL ,ũ RũR ? The second K-factor is commonly available via tools like Prospino [93] or NNLLfast [100] 15 and used by the experimental collaborations for MSSM analyses. The first, however, allows a more direct comparison between the models. This point is depicted in figure 10. The left-hand plot compares the ratios of K-factors (MRSSM over MSSM) for the same process, i.e. pp →ũ LũR which is equivalent to the ratio of NLO cross sections. It therefore allows to view the effects of the sgluons and Dirac-gluinos at NLO. On the other hand, the right-hand plot contrasts the K-factors for full u-squark production, i.e. K(MSSM) is the usual K-factor includingũ LũL andũ RũR production. Thus, this plot does not solely illustrate NLO effects but also the conservation of R-charge, already present at LO.
On the left of figure 10 it can be seen that the sgluon and Dirac mass effects can lead to a difference of the K-factor of five to eight per cent between the MRSSM and the MSSM when considering only theũ LũR production. The ratio is smaller than one for mq mg and larger than one in the other parts of parameter space with the maximum at the ratio of mg/mq = 2. However, when the K-factor for the summed production in the MSSM is used the ratio deviates significantly from one. Using the usual, summed K-factor of the MSSM to estimate NLO MRSSM cross sections from LO ones would consequently lead to a systematic underestimate between 10% and 23%.  We have described the qualitative differences of the NLO corrections between the MRSSM and MSSM in the previous section, especially highlighting the role of the Dirac gluino mass and the appearance of the sgluon. In the following, we give an overview of the quantitative features of the corrections in the MRSSM analysing the variation of the K-factor. Figures 11 and 12 summarise the dependency of the K-factor on the different masses of the strong sector forqq andqq † production, respectively. The left (right) plot of each figure is given for a sgluon mass of 3 TeV (300 TeV) and shows the K-factor depending on the Dirac gluino and common squark masses.

Squark production at NLO
As discussed before, the change of the sgluon mass from 3 to 300 TeV leads to a global enhancement of the K-factor of around twenty per cent in the whole parameter plane for qq production which originates from the super-oblique corrections. The increase for theqq † production is reduced to about five per cent as only part of the contributions receive the relevant corrections.
Forqq production, the total NLO contributions can decrease the LO cross section by 10% and can increase them by more than 50% of the LO prediction assuming sgluons are close in mass to gluino and squarks. The relative size of the corrections falls with rising gluino mass while increasing squark masses lead to an enhancement. This feature is already present in the MSSM [17] and not influenced by the Dirac nature of the gluino or the presence of the sgluon. In the scenario of interest, with a heavy Dirac gluino and rather light squarks, the size of NLO corrections is reduced compared to the remainder of the parameter space.
The K-factors forqq † production are in general smaller than forqq production. They yield up to 30% corrections to the LO production cross section. The corrections are largest for small squark and gluino masses. For small squark masses, the gluino mass does not influence the K-factor significantly. In this region pure QCD corrections are dominant and only with a large sgluon mass do the effects described in section 6.1.2 become important. The K-factor is smallest for large squark masses and increases in this parameter region with the gluino mass.

Full results including uncertainties
We illustrate in figure 13 the summarising results of our study of the NLO SUSY-QCD corrections for the production of squarks in the MRSSM. The plots show LO and NLO cross sections forqq andqq † production and their dependence on the squark and gluino mass as well as the corresponding K-factors. The uncertainty bands are calculated by summing the  uncertainty of the scale dependence 16 and the PDF uncertainty in quadrature. A detailed study of both sources of uncertainty is described below.
As we have seen before and has been known from the MSSM, the plots highlight that NLO QCD corrections are in general large and usually positive. Comparing the LO and NLO predictions including their respective 68% CL uncertainty bands we see that the bands overlap for most of the parameter space and the NLO corrections show no unexpected behaviour.
Comparing the cross sections for the MRSSM calculated with NLO precision to the ones of the MSSM, the distinctions between both models are similar to the ones already discussed at tree level in section 3. These differences stem from the conservation of R-charge in the MRSSM so that only different (same) "chiralities" of squarks can be produced inqq (qq † ) production as has been discussed in detail in section 6.1.3.
Effects from the existence of sgluons becomes relevant only for sgluon masses above 16 To estimate the scale dependence, we varied both, µF and µR like described in section 6.3.1.  hundreds of TeVs. The Dirac nature of the gluino in the virtual contributions is of influence in the region of light squarks and large gluino mass but does not alter the behaviour of the cross section significantly compared to the tree-level differences.
The relative uncertainties of the cross sections are reduced from 50% or more to below 20% when going to NLO. The dominating uncertainty component comes from the scale variation. Therefore, a further reduction of the uncertainties would require going to next-to NLO and/or including effects from the resummation of threshold corrections as has been done for the MSSM predictions [102,103].
We show in figure 13 the K-factor defined in equation (6.1) using the LO PDF sets for the LO cross section. Alternatively, it is also possible to use the NLO PDF sets for both, NLO and LO cross section. The difference is studied based on table 3, where the LO cross sections and NLO K factors using either LO or NLO PDF sets are given for the benchmark points of table 2 and both squark production processes of this paper.
In general, using the NLO PDF sets for the LO cross section leads to an enhancement of the K factor (reduction of the LO cross section) arising from the difference in the strong coupling constant between the two sets. (Here, with α MS s (m Z ) = 0.135 for MMHT2014LO, α MS s (m Z ) = 0.120 for MMHT2014nlo68cl.) This can be seen when comparing the changes in the K-factors between BM2 and BM3. They differ by the choice of squark mass which is used as renormalisation and factorisation scale relevant for the extraction of the correct α MS s (µ R ). For BM2 with a larger squark mass the difference in K-factors is enhanced compared to BM3.
Additional effects due to difference in the PDF fits for quark and gluon are also present. This leads to an enhancement ofqq † compared to theqq production with rising squark mass.
For the results of this work we assume that all squark masses are at the same value. We also investigated the scenario, where the constraint of equal squark masses was dropped. We checked that the quantum corrections are well behaved in this case. The phenomenological effects of such a scenario, especially for final state squarks of different masses, will be studied in forthcoming work. In the following, we describe in detail the scale and PDF uncertainties as well as differential distributions of the NLO SUSY-QCD corrections.   Figure 14: Differential cross section ofqq andqq † production as function of the transversal momentum p T or the pseudo-rapidity η for BM1. The integration error is given.

(Fixed) Renormalisation and factorisation scale dependency
A major motivation to calculate perturbative corrections to the prediction of a cross section is the reduction of the theoretical uncertainty. One way of quantifying this is the variation of renormalisation and factorisation scale. Figure 15 shows the variation of scales forũ LũR andũ Lũ † L production. To achieve a qualitative understanding renormalisation and factorisation scale are set equal and varied together. The prediction of the cross section at LO changes by more than 50% when the scale is varied from the reference value at µ R = µ F = mq by a factor of two. This variation uncertainty is reduced to below 20% when taking NLO corrections into account. This difference in the scale dependency also leads to a strong dependence of the K-factor, as shown in the lower plot, which can become smaller than one for a certain choice.
The commonly accepted range for the variation of the scales is [µ R,F /2, 2 · µ R,F ] where, for the most conservative approach, both scales are varied independently and one takes the envelope of those nine values as final estimate. This is the procedure we use to estimate   Table 4: Cross section (in fb) ofũ LũR production for the benchmark points of table 2 using different PDF sets. The first set is the standard NLO PDF of this paper. The other columns compare different PDF sets at α s (m Z ) = 0.118. The factorisation and renormalisation scale are set to the common squark mass. All uncertainties are given at 68% CL. For this the uncertainties of the CT14 set have been rescaled by 1.642 as the default is given at 90% CL. [104] the scale uncertainty in figure 13.

PDF uncertainty
Another main source of uncertainty for the cross section prediction at a hadron collider comes from the use of PDFs as they can not be calculated from first principles but need to be extracted from data. Details on this and the application at LHC Run II can be found in ref. [104]. Here, we only aim to achieve a basic understanding of the relevant PDF uncertainties. For this we compare the cross section for our three BMPs defined in  and both processes calculated with our default set, MMHT2014nlo68cl, against one of the same group with a different fit value for α s (m Z ), MMHT2014nlo68clas118. Additionally, we compare to cross sections calculated using the CT14 [105] and NNPDF3.0 [106] PDF set with α s (m Z ) = 0.118. They are summarised in the tables 4 and 5 .
The difference of the result for the MMHT2014nlo68cl and MMHT2014nlo68clas118 sets stems only from the strong coupling constant, α s (m Z ) = 0.120 and α s (m Z ) = 0.118 respectively, and gives an estimate for the uncertainty coming from this input parameter. For the considered points the uncertainty is at most a few per cent. The different PDF sets for α s (m Z ) = 0.118 all agree within their uncertainty as expected from the comparison done in ref. [104]. The uncertainty of the individual sets is below five per cent forqq production. This is understandable as for this mainly the up quark PDF is relevant which is rather well determined, see figure 6 of ref. [104]. The PDF uncertainty forqq † production on the other hand is increased and can reach 10% to 15% depending on the PDF set. This originates from the appearance of initial gluons and antiquarks as in this case the PDFs are not as well known, see figure 5 of ref. [104].
As the study of PDF uncertainties is not in the focus of this work, the PDF uncertainty in figure 13 contains only the uncertainty of the MMHT2014nlo68cl set.

Dynamical scales and differential distributions
The main concern of this work was to perform the calculation of NLO corrections for the MRSSM and highlight the physical differences to the MSSM. We addressed it by focusing on global K-factors and using fixed renormalisation and factorisation scales. Here we give a brief comment on going beyond that. The MG5aMC@NLO framework allows in principle to use dynamical renormalisation and factorisation scales (as opposed to fixed ones studied in section 6.3.1) and to study differential distributions at the NLO. We have verified that these components work as expected. The global K-factors vary by less than five per cent when choosing between a fixed (like the common squark mass) and a dynamical (like the total transverse mass in an event) scale. As there is no physically preferred scheme we follow the works on the MSSM [94][95][96][97] and use fixed renormalisation and factorisation scales.
In figure 14, we show sample differential distributions in the transversal momentum p T and the pseudo-rapidity η for the benchmark point BM1. The top row shows it forqq while in the bottomqq † production is given. The K-factor can actually vary substantially over the different kinematic regions changing by 50 per cent or more. Still, for the regions containing the bulk of the contributions (500 GeV < p T < 1000 GeV and −2 < η < 2) the global K-factor is a good approximation.

Conclusions
The Minimal R-symmetric Supersymmetric Standard Model is a well-motivated, viable model. It provides a realisation of SUSY distinct from the familiar MSSM with a significantly modified coloured sector: The gluino carries R-charge and is of Dirac instead of Majorana nature, and there are scalar and pseudoscalar colour octets, the sgluons. Left-and right-handed squarks carry opposite R-charge. In a motivated region of parameter space, the squark masses are around the TeV-scale while the gluino and sgluons are somewhat heavier.
Here we have presented an analysis of squark production at the LHC in the MRSSM taking into account NLO corrections of the entire strongly interacting sector. Both the tree-level results and the NLO corrections show important features and differences to the MSSM: • In the MRSSM only opposite-chirality squarks can be pair-produced, i.e.q LqR production is allowed, whileq LqL andq RqR are forbidden by R-charge conservation. For squark-antisquark production, the converse statement is true. Owing to this, the overall squark production rate in the MRSSM is lower than in the MSSM.
• When comparing K-factors between the two models one has to be careful to distinguish in the MSSM between the K-factor for total squark production and the K-factor specific for e.g.q LqR production. The first is the one usually quoted; it differs from the MRSSM counterpart by 10 to 20%. The second is usually not quoted in the MSSM, but it is more directly comparable between the two models.
• Even the comparison betweenq LqR production in the MRSSM and MSSM reveals several differences between the two models. The MRSSM Dirac gluino enters the loop corrections differently from the MSSM Majorana gluino, and large sgluon masses lead to non-decoupling, logarithmic enhancements of the NLO corrections. To a lesser degree, these effects also lead to differences in theqq † production.
Our study also has technical aspects. During the last decade several efforts lead to fast and automatised ways to evaluate SM processes at NLO and beyond. In recent years, there has been a push to expand this also to BSM models. Usually, this added capability is tested in the context of well-known and understood models like the THDM or the MSSM. In this paper we were able to achieve full agreement between a calculation based on a subset of publicly available tools and an independent calculation for NLO QCD corrections of the MRSSM. The latter calculation uses the techniques, described as method 1 in this paper and results will be available via the soon to be published RSymSQCD code [107]. The MRSSM provides a valuable test of the implemented mechanisms as it contains unique features complementary to the previous test cases. It therefore shows how far the latest available machinery can be pushed. It also illustrates the possibility mentioned in ref. [12,19], to compute QCD corrections to SUSY processes in different renormalisation schemes.
The results obtained here are not specific for the MRSSM but can be applied more generally also to other models with R-symmetry and a difference in R-charge assignment. This only requires that the particle content under consideration is the same and all the SUSY QCD vertices given in the appendix are present.
It will be now of high interest to identify the allowed squark mass range in the MRSSM in the light of current LHC data. In line with simplified studies [9,11] we expect that significantly lighter squarks are allowed in the MRSSM, thanks to the fewer allowed production channels. On the other hand, our NLO results show that the assumption that the K-factors in the MSSM and MRSSM are equal is not correct. The MRSSM K-factors are higher and dependent on the hierarchies between squark, gluino and sgluon masses and should be taken into account to obtain valid limits.
In SUSY models, it is often not possible to define a consistent, conserved fermion number. 17 To compute Feynman diagrams in an unambiguous way, a fermion flow can be introduced [108,109], which does in general not agree with the flow of fermion number.
Even though there are no Majorana particles in the MRSSM, it is useful to adapt this procedure. The reason for this is the existence of fermion number violating processes like qq →qq, which results in a clash of arrows in the associated Feynman diagrams.
In the following we list Feynman rules for the MRSSM which are new or different to the ones in the MSSM. We labelled the lines with the corresponding quantum field operators of the respective Lagrangian term (thus this labeling does not coincide with the labels of external particles in section 3 and 4). The fermion flow on vertices is always directed from an unbarred to a barred spinor field and indicated by an extra arrow next to the diagrams. For calculations involving fermions one needs to multiply the Feynman rules in the opposite direction of the fermion flow.
The Feynman rules 4b and 5b are the complex conjugates of 4a and 5a , respectively. Applying a flipping rule to a vertex one has to reverse the curved arrow, i.e. the fermion flow and replace Ψ with Ψ C . In addition one has to add a minus sign for Feynman rule 1 . 17 For Dirac particles like quarks, the direction of fermion number flow is given by the arrow on the Dirac propagator.

B The MRSSM in MadGraph and GoSam
To allow the use of MadGraph5_aMC@NLO (MG5aMC@NLO) together with GoSam for the MRSSM, several changes are necessary in the programs. For MG5aMC@NLO: • Replacing the path to the model in the template gosam.rc of MG5aMC@NLO, • Modifying the write_lh_order method in the madgraph/iolibs/export_fks.py file to produce a LH file conforming to BLHA2 instead of BLHA1.
• Adapting the SubProcesses/BinothLHA_OLP.f also from BLHA1 to BLHA2 by using the OLP_SetParameter function to pass α s from MG5aMC@NLO to GoSam and the OLP_EvalSubProcess2 function to get the virtual matrix element from GoSam.
For GoSam: • Adjusting the naming of the strong coupling in the olp_module.f90 template for the aMC interface to the UFO model name so that it compiles.
• The default SM renormalisation of the virtual matrix element is switched off and replaced with a model and subprocess specific renormalisation in matrix.f90.

C Validation
The selected parameter and phase space point is the same for all MRSSM processes: Real parts contain soft-collinear mass factorisation counterterms. The agreement of the finite parts are not as precise as the one of the poles. This is due to a limited precision of our LoopTools installation. If however double precision is replaced with quadruple precision the agreement is improved.