Soft gluon resummation for associated gluino-gaugino production at the LHC

We perform a threshold resummation calculation for the associated production of gluinos and gauginos at the LHC to the next-to-leading logarithmic accuracy. Analytical results are presented for the process-dependent soft anomalous dimension and the hard function. The resummed results are matched to a full next-to-leading order calculation, for which we have generalised the previously known results to the case of supersymmetric scenarios featuring non-universal squark masses. Numerically, the next-to-leading logarithmic contributions increase the total next-to-leading order cross section by 7 to 20% for central scale choices and gluino masses of 3 to 6 TeV, respectively, and reduce its scale dependence typically from up to $\pm12$% to below $\pm3$%.

1 Introduction Supersymmetry (SUSY) is an attractive extension of the Standard Model (SM) of particle physics. As the maximal space-time symmetry, it relates bosons to fermions and predicts the existence of spin partners of the SM particles that lead to a stabilisation of the Higgs mass, to the unification of the gauge couplings at high energies, and to a viable dark matter candidate. In many scenarios of the Minimal Supersymmetric Standard Model (MSSM), the dark matter candidate is the lightest neutralino, a mixed fermionic state composed of the superpartners of the photon, the Z-boson and the CP -even neutral Higgs bosons. The search for SUSY particles is consequently an important research focus at the LHC. Squarks and gluinos would be most copiously produced there through the strong interaction, and their masses could therefore already be constrained by ATLAS and CMS in Run I and the first year of Run II to lie beyond 1 TeV [1,2]. In contrast, pairs of sleptons and gauginos would be produced electroweakly. Their mass limits are therefore still considerably lower, i.e. in the range of a few hundreds of GeV [3,4]. In most cases, the LHC analyses are based on simplified scenarios with cascade decays of the squarks and gluinos to jets and leptonic decays of the sleptons and gauginos, all accompanied by missing transverse energy from the escaping lightest neutralino.
The experimental analyses rely on the availability of precise theoretical predictions for the production cross sections. For many years, the state of the art were next-to-leading order (NLO) calculations, which typically lead to an increase over the leading-order (LO) cross section [5][6][7] and a reduction of the theoretical uncertainty due to a stabilisation of the renormalisation and factorisation scale dependence [8][9][10][11][12][13][14][15][16][17][18][19]. Since SUSY particles are often produced close to threshold, large logarithms can spoil the convergence of the perturbative series. The current state state of the art is therefore to include also the resummation of leading logarithms (LL) and next-to-leading logarithms (NLL), that has been developed originally for SM processes [20][21][22][23][24], in the context of slepton production [25][26][27][28][29], gaugino production [30][31][32][33][34][35], as well as squark and gluino production [36][37][38][39][40][41][42]. In some cases, also next-to-next-to-leading logarithms (NNLL) and beyond have been resummed [39,[41][42][43][44][45][46][47][48][49], partly using soft-collinear effective theory (SCET). SCET and perturbative QCD results have been compared analytically and numerically, e.g., in Refs. [50,51]. Since the full NNLO results and two-loop matching coefficients for gluino-gaugino associated production are unknown, we for consistency present results only at the NLO+NLL level, even though the two-loop soft anomalous dimensions are in principle known for arbitrary processes in QCD, but are in practice not trivial to extract for our specific process. We therefore have to leave a consistent NNLO+NNLL calculation for future work. Nevertheless, an estimate of approximate NNLO+NNLL vs. NLO+NLL effects can be obtained from a recent calculation for stop pair production [49], where they typically (for a stop mass of 2 TeV and an LHC energy of 13 TeV) amount to an increase of the total cross section by 5% and a further stablisation with respect to scale variations. Note that NLO calculations have also been combined with parton showers for a variety of SUSY [52][53][54][55][56] and GUT processes [57][58][59][60], and these calculations usually agree well with resummation calculations within the theoretical uncertainties.
In this paper, we present a threshold resummation calculation for the associated production of gluinos and gauginos at the NLL+NLO accuracy. This is one of two channels (the other one being the associated production of squarks and gauginos, left for future work) for which NLO calculations have been computed previously [11][12][13][14], but where a resummation calculation has so far not been performed. Its production cross section is of intermediate strength, as it involves both strong and weak couplings. It can become phenomenologically relevant in particular in the case that gluino pair production is beyond the current LHC reach due to an exceedingly large gluino mass. This could very well be realised in Nature, as in Grand Unified Theories (GUTs) one expects the gaugino mass parameters M i to unify, similarly to the corresponding gauge couplings α i , so that after renormalisation group running M 3 ≃ 6M 1 at the weak scale. The gluino with mass mg = M 3 is then typically much heavier than the electroweak gauginos with masses of the order of the bino (M 1 ), wino (M 2 ≃ 2M 1 ) or higgsino (µ) mass parameters.
The outline of this paper is as follows: In Sec. 2, we describe briefly our analytical calculations at LO and NLO and the refactorisation of the cross section, then in more Figure 1. Tree-level t-(left) and u-channel (right) Feynman diagrams for the associated production of a gluino and a gaugino at hadron colliders. The dashed lines represent squark exchanges.
detail the necessary calculation of the hard matching coefficient as well as the matching procedure and inverse Mellin transform. Our numerical results are presented in Sec. 3 for a typical benchmark scenario that satisfies Higgs mass, flavour-changing neutral current, muon magnetic moment and LHC constraints. Our conclusions are given in Sec. 4. The coupling conventions employed in this paper are listed in App. A, and the calculation of the soft anomalous dimension is presented in App. B.

Soft gluon resummation
We start the presentation of our work with a description of our analytical results. Our LO and NLO calculations are presented in Sec. 2.1 and verified to agree with those obtained previously in the limit of degenerate squark masses. Sec. 2.2 describes briefly the refactorisation and resummation formalism up to the NLL accuracy. In Sec. 2.3 and App. B, we present in more detail the required calculations of the process-dependent hard matching coefficient and the soft anomalous dimension. Both are analytically checked against each other in Sec. 2.4, where also the matching to the NLO calculation and the inverse Mellin transform are performed.

Production of gluinos and gauginos at leading and next-to-leading order
At hadron colliders, the associated production of a gluino and a gaugino with four-momenta p 1 and p 2 and masses mg and mχ proceeds at leading order (LO) through the annihilation of a quark and an antiquark, both taken as massless, with four-momenta p a and p b , Here, we distinguish possibly different quark flavours with a prime and label the gaugino mass eigenstate by the index j (j = 1, ..., 4 for neutralinos and j = 1, 2 for charginos).
The contributing processes are shown in Fig. 1. They are mediated by the exchange of a virtual squark (q) in the t-(left) and u-channel χ are the usual partonic Mandelstam variables. The corresponding squared matrix elements are given by where the label c of the squark masses mq and couplings L ( ′ ) , R ( ′ ) , L ( ′ ) and R ( ′ ) refers to those appearing in the complex conjugate diagrams. Our conventions for the notations of couplings can be found in App. A. The SU(3) colour factors are C A = 3 and C F = 4 3. The relation between the electromagnetic and weak couplings e = g sin θ W depends on the weak mixing angle θ W , and g s (µ r ) is the (renormalisation-scale dependent) strong coupling constant. The total spin-and colour-averaged squared amplitude is then where the sum is performed over all squarks in the propagators and where we have inserted a relative minus sign between the t-and u-channel for the crossing of a fermion line. We have hence generalised the results of the literature [5][6][7][8][9][10][11] by allowing for arbitrary squark mixings and mass eigenstates in the squark couplings and propagators, as we have already done in our calculations for slepton [25][26][27][28][29] and gaugino pair production [30][31][32][33][34][35]. This will in particular allow us to extend existing studies of SUSY flavour violation [61][62][63][64][65][66] to new processes. Integration over the two-particle phase space dPS (2) = dt (8πs) leads to the total partonic cross section and convolution with the factorisation-scale dependent parton distribution functions (PDFs) to the total hadronic cross section [67]. Here, τ = M 2 S denotes the ratio of the squared invariant mass of the produced SUSY particle pair and the hadronic centre-of-mass energy and is related to the partonic momentum fractions x a and x b by z = τ (x a x b ) with z = 1 at LO. For non-mixing squark exchanges, the NLO corrections are well-known [11][12][13][14][15]. They involve one-loop self-energy, vertex correction and box diagrams interfering with those at tree level as well as squared real gluon and quark emission diagrams, which can involve intermediate on-shell squarks. We have recalculated the full NLO cross section for general squark mixings and mass eigenstates using dimensional regularisation of ultraviolet and infrared divergencies, MS renormalisation for couplings and wave functions substituted with a finite shift in the quark-squark-gluino Yukawa coupling to restore SUSY [68], onshell renormalisation for all squark and gluino masses, and the dipole subtraction method to deal with infrared and collinear divergences [69,70]. This method exploits the fact that the NLO cross section σ (1) can be split into two parts σ {2} and σ {3} with two-and threeparticle kinematics, respectively, and a collinear counterterm σ C , that removes initial-state collinear singularities, (2.8) The two-and three-particle cross sections are individually regularised by subtracting from the real corrections σ R an auxiliary cross section σ A . The latter captures all infraredsingular behaviour, can be analytically integrated over the singular phase space regions, and is added back to the virtual corrections σ A . In the limit of mass-degenerate squarks, we achieve full numerical agreement with our previous calculation [12][13][14], that employed the phase-space slicing method, and the public code Prospino [8].

Refactorisation
Beyond LO, z = M 2 s describes the energy fraction going into the hard scattering process. The energy fraction carried by an additional emitted gluon (or quark) is then 1 − z and describes the distance from partonic threshold. As is well known, large logarithms α s 2π n ln m (1 − z) with m ≤ 2n − 1 remain at higher orders in α s = g 2 s (4π) even after the cancellation of soft and collinear divergences among the virtual and real emission corrections [71,72]. They arise due to restrictions on the phase space boundary of the latter and spoil the convergence of the perturbative series when z → 1, i.e. for soft emitted gluons. They must therefore be resummed to all orders in α s for reliable predictions in the threshold region.
Resummation of these logarithmic corrections to all orders and exponentiation can be achieved when the kinematical and dynamical parts of the cross section are fully factorised. By applying a Mellin transform to the quantities F = σ AB , σ ab , f a A and f b B with y = τ , z, x a and x b in Eq. (2.7), the hadronic cross section can be written as a simple product Large logarithms in z → 1 then turn into large logarithms of the Mellin variable N .
Applying the eikonal Feynman rules, using renormalisation group equation properties and the evolution equations, the partonic cross section can then be refactorised and written in the resummed form [20][21][22] σ (res.) where we have identified the renormalisation and factorisation scales µ = µ r = µ f for brevity. The hard function which is non-singular when z → 1 or N → ∞ can be expanded perturbatively in a s = α s (2π) and is discussed in more detail in Sec. 2.3. Like the soft wide-angle function ∆ ab→ij,I to be discussed in App. B, it is sensitive to the colour structure of the underlying hard process from which the gluon has been emitted, and can be decomposed into irreducible colour representations I. Together with the functions ∆ a,b describing soft collinear radiation, the soft wide-angle function ∆ ab→ij,I can be exponentiated in a closed form [73] ∆ with λ = a s β 0 L, L = lnN andN = N e γ E , which contains all the enhanced logarithmic terms. The terms G (1) ab and G (2) ab→ij represent the leading-logarithmic (LL) and next-toleading logarithmic (NLL) approximations [20][21][22][23][24] where the one-loop coefficient D (1) ab→ij,I depends on the soft anomalous dimension and is process dependent. While it vanishes for gaugino pair production with coloured particles in the initial state only [32], it must be included for associated gluino-gaugino production, where soft gluons can also be emitted from the final-state gluino. The associated collinear divergence is, however, screened due to the massive emitter. The coefficients in the above equations read where C a = C F,A for quarks and gluons andΓ ab→ij,II is the modified diagonal soft anomalous dimension of the colour representation I.

Hard matching coefficient
The resummation of logarithmically enhanced contributions at threshold can be improved by including in the hard function H ab→ij, but also the N -independent contributions of the NLO cross section which beyond NLO are multiplied by threshold logarithms. The hard matching coefficient C ab→ij (M 2 ) described in Sec. 2.1, keeping only the N -independent terms. In the case of associated gluino-gaugino production, this is again simplified by the fact that there is only one colour basis tensor, so that the index I can be dropped.
In Eq. (2.8), the three-particle contributions to σ (1) can safely be neglected, since all infrared divergences are canceled after subtracting from dσ R the auxiliary cross section dσ A and the finite terms are phase-space suppressed near threshold [41,42]. The virtual corrections dσ V and the integrated dipoles ∫ 1 dσ A are straightforward to transform into Mellin space in the invariant-mass threshold approach, since they are proportional to δ(1 − z) and thus constant in N . The corresponding analytical results can be found in Refs. [13] and [69,70], respectively. The collinear remainder σ C is usually split into contributions from two insertion operators P and K [70] The former is directly related to the regularised Altarelli-Parisi splitting distributions at O(α s ), while the latter depends on the factorisation scheme and on the regular parts of the Altarelli-Parisi splitting distributions. For the initial quark a, we obtain after transforming to Mellin space and similarly for the incoming antiquark b with t → u. Non-diagonal operators give only 1 N -suppressed contributions, and there are no initial-state gluons at LO. In the limit of C A → 0, one recovers the well-known results for Drell-Yan like processes [32].
For the hard matching coefficient, only the N -independent parts of the above results are needed. The logarithmic terms terms can be used to check the resummed cross section when re-expanded to NLO (see below). Since the 1 N terms have been systematically neglected, the collinear improvement of the resummation formalism suggested in Refs. [74,75] has not been performed in contrast to our calculations for uncoloured slepton [27] and gaugino pair production [32], but similarly to the calculations for coloured squarks and gluinos [40][41][42]. This is in particular due to the fact that analytic results for the subleading terms in the Mellin transform of the collinear remainder can not be obtained. For the Drell-Yan process, the constant terms in the hard function H (1) ab→ij (M 2 , µ 2 ) are sometimes also exponentiated based on the argument that they factorise the complete Born cross section, include finite remainders of the infrared singularities in the virtual corrections and are thus related to the corresponding singularities in the real corrections giving rise to the large logarithms [76]. However, similarly to gaugino pair production [32], gluino-gaugino associated production does not proceed through a single s-channel diagram (see Fig. 1) and the virtual corrections factorise only at the level of amplitudes, and not at the level of the full cross section [12][13][14]. Resumming these finite terms is therefore not justified in this case.

Matching and inverse Mellin transform
While near threshold the resummed cross section is a valid approximation, far from it the normal perturbative calculation should be used. A reliable prediction in all kinematic regions is then obtained through a consistent matching of the two results with (2.29) Here, the resummed cross section σ where H (0) and H (1) are the first and second order parts of the hard matching coefficient , the leading logarithms α s (2π) ln 2N have the coefficient 4C F , which agrees with the leading logarithmic contribution to the hard matching coefficient arising exclusively from the K-operators in the collinear remainder in Eq (2.27). The coefficient 4C F also governs the scale-and more precisely the ln(µ 2 f s)dependent part of the next-to-leading logarithms α s (2π) lnN , which agrees with the corresponding parts of the quark and antiquark P -operator expectation values in the collinear remainder in Eq. (2.26). In contrast, the C A -terms depending on µ f cancel there, while the remaining NLL terms are From the K-operators in Eq. (2.27), we get in addition the NLL contributions Having computed the resummed and the perturbatively expanded results in Mellin space, we must multiply them with the N -moments of the PDFs according to Eq. (2.11) and perform an inverse Mellin transform in order to obtain the hadronic cross section as a function of τ = M 2 S. Special attention must be paid to the singularities in the resummed exponents G (1,2) ab , which are situated at λ = 1 2 and are related to the Landau pole of the perturbative coupling a s . To avoid this pole as well as those in the Mellin moments of the PDFs related to the small-x (Regge) singularity f a A (x, µ 2 0 ) ∝ x α (1 − x) β with α < 0, we choose an integration contour C N according to the principal value procedure proposed in Ref. [77] and the minimal prescription proposed in Ref. [78]. We define two branches 21 773 1300 315 1892 2288  425  1758 2754 552  714 1553 -2200   Table 1. Higgs sector and soft SUSY breaking parameters in our pMSSM-13 benchmark model. All values except the one for tan β are in GeV.
where the constant C is chosen such that the singularities of the N -moments of the PDFs lie to the left and the Landau pole to the right of the integration contour. Formally the angle φ can be chosen in the range [π 2, π[, but the integral converges faster if φ > π 2.
The Mellin moments of the PDFs are obtained by fitting to the parameterisations tabulated in x-space the functional form used by the MSTW collaboration [79] f which has the advantage that it can be transformed analytically with the result Here, y = A 2 + 1 and B ′ (x, y) = B(x, y) Γ(y) = Γ(x) Γ(x + y). We have verified that we obtain good fits not only for the MMHT2014NLO118 [79], but also for the CT14NLO fits [80] up to large values of x and for all typical factorisation scales, even though the latter are obtained with an ansatz including an exponential function. The fit to the NNPDF 30 nlo as 0118 PDFs [81] is slightly less stable in the large-x region. They will therefore only be used for estimates of the PDF uncertainty. We compute in this case first an (approximately PDF-independent) K-factor, i.e. the ratio of NLL+NLO over NLO cross sections, using stable (e.g. CT14NLO) PDFs and then multiply with it the NLO calculation convoluted with NNPDFs directly in x space.

Numerical results
We now turn to our numerical results. For our calculations, we have used the Particle Data Group (PDG) values for the Standard Model parameters, in particular for the value of the strong coupling constant at the Z-pole α s (M Z ) [82]. In our LO and NLO/NLL+NLO calculations, it is evaluated in the one-and two-loop approximations, respectively, with five active quark flavours. All light quarks including the bottom quark are taken as massless. The top quark is decoupled. Its (pole) mass enters only in the gluino self-energy and has little numerical influence on the production cross sections.

Benchmark scenario
Our results are given for a specific phenomenological MSSM (pMSSM) benchmark scenario with 13 free parameters. These parameters are listed together with the corresponding fitted numerical values in Tab. 1. Our scenario is inspired by the benchmark point II of Ref. [66],  that has been obtained with a Markov Chain Monte Carlo (MCMC) scan using also PDG values for the Standard Model parameters. This 19-parameter scan has been performed with a focus on non-minimal flavour violation (NMFV). It therefore included seven flavourviolation parameters and checked in particular the most stringent flavour-changing neutral current (FCNC) constraints from rare B-and K-decays. Since we are not interested in NMFV, we set these parameters all to zero. This reduces the SUSY contributions to the rare meson decays. To compensate for the reduced mass splitting in the top squark sector and obtain a Higgs-boson mass compatible with the measured value, we have instead changed the sign and increased the absolute value of the trilinear coupling A f . We then still obtain a neutralino lightest SUSY particle and in addition a light top squark, which leads to a viable dark matter candidate and allows in general for sufficient stop coannihilation to reproduce the observed dark matter relic density [83][84][85][86][87]. While we continue to impose the GUT relation between the bino and wino mass parameters, M 1 ≃ M 2 2, we allow the gluino mass parameter M 3 to vary independently, which brings us to 19 − 7 + 1 = 13 free parameters. For our default scenario, we still impose the GUT relation for M 3 ≃ 6M 1 . The physical SUSY mass spectrum is obtained with SPheno 3.37 [88] and shown schematically in Fig. 2. Apart from the FCNC, the Higgs-boson and neutralino dark matter data, also the observed value of the anomalous magnetic moment of the muon, which is unaffected by the gluino mass, is reproduced in this scenario. In addition, it satisfies the increasingly stringent constraints that are imposed on the masses of the SUSY particles from direct search results at the LHC. For example, with the 2015 data from Run II, ATLAS and CMS exclude gluino masses up to 1400 and 1280 GeV, assuming masses of the lightest neutralino of up to 600 and 800 GeV, respectively [1,2]. Mass-degenerate light charginos and second-lightest neutralinos produced electroweakly have been excluded at Run I up to 465 and 720 GeV in the case of massless lightest neutralinos [3,4]. These exclusion limits are however not valid within the general pMSSM, as they have been obtained assuming direct production cross sections and simplified decay scenarios.
As we have already mentioned in the previous section, our calculations allow for arbitrary squark mixings and mass eigenstates in the appearing couplings and propagators. The mixing of squark interaction eigenstates is numerically relevant only for third-generation (and in particular top) squarks. Since the top and bottom quark PDFs are small, the mixing influences predominantly the gluino self-energy diagram with little (below the percent level) numerical effect. In contrast, the masses of first-and second generation squarks in our scenario span a large range from about 600 to 2300 GeV. Averaging over these masses leads to cross sections that are almost a factor of two larger. This is already true at LO, since the dominant effect comes from squark propagators already present there. If the LO calculations are performed without averaging and corrected by a mass-averaged K-factor, as it is, e.g., done in Prospino [8], differences of about 5% remain. These differences originate in particular from intermediate squarks at NLO, which in the general case can sometimes be on-shell, while they cannot be on-shell in the mass-averaged case.

Invariant mass distribution
The associated production of a gluino and the lightest neutralino will be difficult to observe at the LHC, as the latter escapes directly undetected. It is therefore more promising to study the associated production of a gluino with the second-lightest neutralino (or the lightest chargino of often equal mass), since it will decay into an additional Z (or W ) boson, whose leptonic decay products will then lead to an identifiable signal and better background suppression. As our default PDFs, we use CT14NLO at NLO and NLL+NLO and CT14LL together with the one-loop approximation for α s at LO (see above) [80].
In the upper panel of Fig. 3, we show the invariant-mass distribution given by Eq. (2.33), for the production of a gluino with a mass of 1892 GeV and a second-lightest neutralino with a mass of 630 GeV, where both masses have been chosen such that they lie beyond current LHC limits even in simplified scenarios. The cross sections peak at about 3.2 TeV and then fall off towards higher invariant masses M . Due to additional radiation, the maximum is shifted from LO (blue) to NLO (green) and NLL+NLO (red) towards slightly smaller values of M . At the same time, the scale uncertainties (shaded bands) are significantly reduced from LO to NLO and then again to NLL+NLO. The second reduction is more clearly visible in the lower panel of Fig. 3, where it amounts to a change from ±12% to ±3% at high invariant masses. The scale errors have been obtained here in the usual way from individual variations of the factorisation and renormalisation scales by a factor of two about the average mass of the two produced final-state particles, (mg + mχ0 2 ) 2, excluding relative factors of four. In the high-mass region, the corrections from threshold resummation at NLL increase the central NLO cross section by up to 10% (black line) and more closer to the threshold.

Scale uncertainty of the total cross section
After integrating over the invariant mass M in Eq. (2.33), we obtain the total production cross section for a gluino and a second-lightest neutralino. For a process that depends already at LO on the strong coupling constant, one expects the significant (approximately logarithmic) scale dependence to be already reduced at NLO. This is clearly visible in Fig. 4, where the NLO result (green dashed curve) shows the characteristic maximum at approximately half the central renormalisation and factorisation scale (upper panel). The scale dependence is further reduced at NLL+NLO (red full curve), as it was already the case for the invariant-mass distribution (see above). When the NLL+NLO result is re-expanded to NLO (blue dotted curve), it becomes a good approximation to the full NLO result in particular for large scale choices, when the logarithmic terms dominate the cross section. This is also true when the renormalisation (central panel) and factorisation scales (lower panel) are varied individually and not together. The lower two panels also demonstrate nicely the interplay of the renormalisation and factorisation scale behaviour in the NLO and NLL+NLO cross sections, that together produce the stabilised behaviour in the upper panel with a large plateau in particular at NLL+NLO. We have verified that we obtain similar results also for larger gluino masses of, e.g., 3 TeV.

Gluino mass dependence of the total cross section
Since the gluino mass is unknown, it is interesting to compute the total cross section for associated gluino-neutralino production as a function of the gluino or gaugino mass. The gluino mass dependence is shown in the upper panel of Fig. 5 and in tabular form in Tab. 2 in App. C. As expected, the cross section falls steeply with the gluino mass from 3 to 0.01 fb in the range mg ∈ [500, 3000] GeV. With LHC luminosities of currently a few fb −1 and in the near future a few 100 fb −1 at √ S = 13 TeV, these cross sections will soon be observable. In the high-luminosity phase of the LHC, expected to collect up to 3000 fb −1 , even larger gluino masses can be reached that may kinematically no longer be accessible in the strong production of gluino pairs.
As the lower panel of Fig. 5 shows, the NLO scale uncertainty on the total cross section of ±10% at 500 GeV decreases only slightly towards higher gluino masses. This is in sharp contrast to the NLL+NLO prediction, that has already a smaller scale error of ±7% at 500 GeV and that becomes much more reliable with an error of only a few percent at large gluino masses. At the same time, the NLO cross section is increased at NLL+NLO by 7 (black line) to 20% for gluino masses of 3 to 6 TeV, respectively.

Gaugino mass dependence of the total cross section
Similarly to the gluino mass dependence, the total cross section for gluino-gaugino associated production decreases with the gaugino mass. Since the mass of the second-lightest neutralino (and the almost identical one of the lightest chargino) is a dependent physical mass obtained after diagonalisation of the neutralino (or chargino) mass matrix with the mixing matrix N , we vary instead the bino mass parameter M 1 , which fixes immediately also the wino mass parameter M 2 through the GUT relation M 2 ≃ 2 M 1 . As one can see in the upper panel of Fig. 6, the mass eigenvalue of the second-lightest neutralinoχ 0 2 increases linearly with M 2 up to M 2 = µ = 773 GeV, where a typical avoided crossing occurs. At higher values of M 2 , it is the mass eigenvalue ofχ 0 4 that depends linearly on M 2 , while the mass ofχ 0 2 stays constant. Accordingly, the decomposition ofχ 0 2 changes from winotype (large mixing matrix element N 22 , red curve) to higgsino-type (large mixing matrix elements N 23 and N 24 , yellow and blue curves). Both features will of course influence the production cross section of the process pp →gχ 0 2 . The dependence of the cross section on the wino mass parameter is shown in Fig. 7 and in tabular form in Tab. 3 in App. C. As expected, it falls with M 2 as long as the physical mass ofχ 0 2 changes. For M 2 > µ, the neutralino becomes higgsino-like and couples mostly via quark Yukawa couplings. Since the heavy quark PDFs in the proton are small, the cross section starts to fall even faster than before despite the fact that the gaugino mass remains constant and the available phase space no longer changes. Interestingly, the cross section increases again somewhat for M 2 > 1150 GeV, which can be explained with a slightly increasing bino component (see Fig. 6). The scale uncertainty (shaded bands) is drastically reduced from LO (blue) to NLO (green), in particular at low values of M 2 . However, the NLO scale dependence increases with M 2 up to M 2 = µ = 773 GeV due to the rising contribution of the large logarithms. This can be seen more clearly in the lower panel of Fig. 7. Beyond this mass, the scale dependence remains constant as expected. A similar trend is observed in the NLL+NLO scale dependence (red), albeit at an again much lower level. To be specific, it rises only from ±1 to ±3% compared to a scale dependence at NLO that rises from ±3 to ±6%. The K-factor (black line) increases also with M 2 from 1.02 to 1.07.

Parton density uncertainty of the total cross section
While the scale uncertainty is expected to be reduced due to the resummation of large logarithms, the PDF uncertainty is normally not improved by this procedure. It is usually estimated by propagating the experimental uncertainties on the fitted data points through to the PDFs (e.g. linearly via a Hessian method) and leads to the production of orthogonal eigenvector PDF sets corresponding to a 90% confidence level. This method is employed by the CTEQ and MMHT collaborations and has been found to produce comparable results to the more intricate Lagrange multiplier method [79,80]. The uncertainty on the cross section is then obtained from where n = 28 and 25 is the number of eigenvector directions in the CT14 and MMHT2014 fits, respectively. Since the MMHT PDF sets are only available for 68% and not 90% confidence level, the corresponding error must be multiplied by the standard factor of 1.645 for compatibility. The NNPDF collaboration uses instead a Monte Carlo method, where the PDF uncertainty is obtained by sampling the available replicas [81]. Since PDF uncertainties are usually not produced for LO fits and since the results are very similar at NLO and NLL+NLO, we will only study them at the level of NLL+NLO cross sections. For NNPDF, we estimate the PDF uncertainty using the K-factor method described in Sec. 2.4. The results are shown in Fig. 8. As one can see, the estimates by the three groups overlap to a large extent. While the three individual PDF uncertainties are still relatively small and comparable to the scale uncertainty with about ±5% at small values of the gluino mass of about 500 GeV (cf. Fig. 5), they increase rather than decrease with mg and reach a level of ±20% at 3 TeV. This is of course due the fact that the PDFs are much less constrained at large than at intermediate values of the parton momentum fraction x. Compared to the central prediction with CT14, those with MMHT2014 and NNPDF30 lie systematically higher by a few percent. The increase of the uncertainty towards larger x is less pronounced in MMHT and more pronounced in NNPDF. It will therefore be interesting to study the impact of threshold-improved PDFs in future work, as it was done for squark and gluino production [89,90]. It is important to note that the scale uncertaintites computed in the previous sections and the PDF uncertainty computed in this section are independent and therefore usually added in quadrature for a reliable estimate of the total theoretical uncertainty.

Conclusion
We have presented in this paper a threshold resummation calculation at the NLL+NLO accuracy for the associated production of gluinos and gauginos at the LHC. This process is of intermediate strength compared to the strong production of gluino (and squark) pairs and the electroweak production of gaugino (and slepton) pairs. It can in particular become relevant should the gluinos prove to be too heavy to be pair-produced at the LHC. This situation would not be unexpected, if one takes the GUT relation between the gaugino masses seriously, which predicts M 1 = M 2 2 = M 3 6 after renormalisation group running at the weak scale. Lighter gluinos, e.g. with a possible cosmological impact through their coannihilation with gauginos, typically require non-minimal assumptions such as non-universal gaugino masses or vector-like supermultiplets [91][92][93][94]. Conversely, associated gluino-gaugino production could also become phenomenologically important in the less likely case that the gauginos lie beyond the kinematic reach of the LHC, so that e.g. the gravitino becomes the lightest SUSY particle and its associated production with (relatively light) gluinos an interesting search channel [7].
Our investigations required the (re-)calculation of the full NLO corrections, which we generalised to the case of non-degenerate squark masses, and of the process-dependent soft anomalous dimension and hard matching coefficients, which we could show to be consistent with each other. For a typical benchmark scenario, obtained in a recent MCMC fit of the Higgs boson, FCNC, muon magnetic moment and LHC data, we presented numerical predictions for the invariant-mass distribution and the total cross section as a function of the gluino or gaugino mass. The resummation of the NLL contributions increased the NLO cross sections at large invariant mass by up to 10% and stabilised them dramatically with respect to the scale dependence. As expected, the PDF uncertainty was, however, not reduced. It will therefore be interesting to study the impact of threshold-improved PDFs in the future.
Numerical predictions for other SUSY scenarios are available from the authors upon request. The calculation will also be included in the next release of the public code Resummino [35]. Its application to the associated production of gluinos and gravitinos would not only require the inclusion of the additional s-channel gluon exchange, but also a full NLL+NLO calculation for the gluon-initiated diagrams. The NLL+NLO calculation for the associated production of squarks and gauginos is in progress and will be presented elsewhere.

A Coupling conventions
The conventions for the couplings appearing in our calculations are defined in Fig. 9. The electromagnetic and (renormalisation scale dependent) strong coupling constants are denoted by e and g s (µ r ), respectively. T a βα and f abc are SU(3) colour matrices in the fundamental representation and structure constants, γ µ and P L,R Dirac matrices and chirality projection operators. The latter are associated with generic MSSM coupling constants L ( ′ ) , R ( ′ ) , L ( ′ ) and R ( ′ ) which involve squark and gaugino mixing matrices and can be found together with the quartic squark couplings X and Y in Refs. [95,96].
−g s f abc γ µ Figure 9. Interaction vertices appearing in the associated production of gluinos and gauginos at LO and NLO. All momenta are ingoing, and arrows describe charge/fermion flow.

B Soft anomalous dimension
If the calculation of the (modified) soft anomalous dimension is performed in the axial gauge with a gauge vector n µ , it is given bȳ where one sums over the two incoming particles and where n 2 = −n 2 − i [21,40]. The dimensionless vector v k is given by the momentum of the incoming massless particle k rescaled by 2 s. Here, the soft anomalous dimension has been modified (subtracted) for the soft functions of the two incoming Wilson lines annihilating into a colour-singlet, i.e the Drell-Yan process, effectively isolating the gauge dependence of a single line and making the soft functions separately gauge invariant. Soft anomalous dimensions are computed from the renormalisation constants Z IJ of Wilson-line operator products by taking the residues of their ultraviolet poles in = 4 − D in D dimensions, Here, we only need the one-loop corrections, depicted in Fig. 10. If we denote by k and l the eikonal lines, between which the gluon is spanned, by C kl IJ the colour mixing factors and by ω kl the kinematic parts of the one-loop corrections, we obtain for the correction to the colour basis tensor c I Γ ab→ij, ω b1 ω a1 Figure 10. One-loop diagrams contributing to the soft anomalous dimension for the associated production of massive colour-octet gluinos and colour-singlet gauginos with momenta p 1 and p 2 from massless colour-triplet quarks and antiquarks with momenta p a and p b . The self-energy contributions of the latter vanish.
We compute these diagrams in an irreducible s-channel colour basis with tensors c J , which one obtains in general after decomposing the reducible initial-state or final-state multi-particle product representations and which has the advantage of rendering the anomalous dimension matrices diagonal at threshold [21,40]. Since we have only one coloured final-state particle, the gluino with adjoint colour index i, i ′ , there is also only one colour basis tensor c J = Tr(T i T i ′ ) = T F δ ii ′ with T F = 1 2, similarly to the case of prompt photon production with an associated gluon jet [97], and we can drop the associated indices I, J. At LO, it leads to the colour factor Tr(T i T i ) = T F δ ii = C A C F computed in Sec. 2.1. At one loop and after removing the LO colour factors, we obtain for the four diagrams in Fig. 10 where j is the colour index of the exchanged gluon. The kinematic part of the one-loop corrections can be written in a general form as where q is the loop momentum, ∆ k,l and δ k,l are signs associated with the eikonal Feynman rules, and P stands for the principle value [21,40] The integrals can be solved [21] (also for two coloured final-state particles with unequal masses [40]) with the results [98] Here, we have combined the signs in S kl = ∆ k ∆ l δ k δ l , so that S ab = 1, S a1 = 1, S b1 = −1, and S 11 = −1. The double poles in in the first three integrals involving at least one massless particle have canceled among themselves. As one can easily see, the scalar products are The function L k = [L k (+n) + L k (−n)] 2 depends in a rather complicated way on the gauge vector n [21,98]. However, all gauge-dependent terms disappear after the inclusion of the self-energies of the two incoming Wilson lines. Combining colour factors, signs, soft integrals and simplifying the result leads tō Due to the LSZ reduction formula, only half of the self-energy contribution ω 11 has been taken into account. All terms proportional to C F have vanished, so that only terms proportional to C A remain. In the massless limit and before subtracting the initial-state self-energies, one recovers the well-known result for associated gluon-photon production, i.e. the C A -term in Eq. (2.26) of Ref. [97]. Our modified soft anomalous dimension can also be compared to the one obtained for associated top-quark and W -boson production in Eq. (3.8) of Ref. [99] after adjustments of the colour factors. The result in Eq. (3.1) of Ref. [100] is slightly different, since Feynman gauge and not axial gauge was used there. The final result for the soft wide-angle emission function in associated gluino-gaugino production is therefore At the production threshold, where the final-state particle velocities vanish and we find D qq→gχ = −C A , (B. 18) in accordance with Ref. [40].

C Lists of total cross sections
In Tabs. 2 and 3 we list the total cross sections for the associated production of a secondlightest neutralino and a gluino at the LHC with a centre-of-mass energy of 13 TeV in tabular form. In Tab. 2, these cross sections are presented as a function of the gluino mass, while in Tab. 3 they are presented as a function of the wino mass parameter M 2 . These tables thus correspond to Figs. 5 and 7. They also include listings of the respective scale and PDF uncertainties.   Table 3. Total cross sections for p p →χ 0 2g at √ S = 13 TeV as a function of the wino mass parameter M 2 using CT14NLO PDFs.