The Inclusion of Two-Loop SUSYQCD Corrections to Gluino and Squark Pole Masses in the Minimal and Next-to-Minimal Supersymmetric Standard Model: SOFTSUSY3.7

We describe an extension of the SOFTSUSY spectrum calculator to include two-loop supersymmetric QCD (SUSYQCD) corrections of order $\mathcal{O}(\alpha_s^2)$ to gluino and squark pole masses, either in the minimal supersymmetric standard model (MSSM) or the next-to-minimal supersymmetric standard model (NMSSM). This document provides an overview of the program and acts as a manual for the new version of SOFTSUSY, which includes the increase in accuracy in squark and gluino pole mass predictions.

message.The higher order corrections included are for the MSSM (R−parity conserving or violating) or the real R−parity conserving NMSSM only.CPC Classification: 11.1 and 11.6.Does the new version supersede the previous version?: Yes.Reasons for the new version: It is desirable to improve the accuracy of the squark and gluinos mass predictions, since they strongly affect supersymmetric particle production cross-sections at colliders.Summary of revisions: The calculation of the squark and gluino pole masses is extended to be of next-to-next-to leading order in SUSYQCD, i.e. including terms up to O(g 4 s /(16π 2 ) 2 ).

Introduction
Near the beginning of LHC Run II at 13 TeV centre of mass collision energy, hopes are high for the discovery of new physics.A much-studied and long-awaited framework for new physics, namely weak-scale superymmetry, is being searched for in many different channels.In the several prominent explicit models of supersymmetry breaking meditaion, the best chance of observing the production and subsequent decay of supersymmetric particles is via squark and/or gluino production.Squarks and gluinos tend to have the largest production cross-sections among the sparticles and Higgs bosons of the MSSM, because they may be produced by tree-level strong interactions, as opposed to smaller electroweak cross sections relevant for the other superparticles.In R−parity conserving supersymmetry (popular because of its apparently viable dark matter candidate), squark and/or gluino production may result in a signal of an excess of highly energetic jets in conjunction with large missing transverse momentum with respect to Standard Model predictions.This classic LHC signature was searched for at LHC Run I at 7 and 8 TeV centre of mass energy, but no significant excess above Standard Model backgrounds was observed.Jets plus missing transverse momentum channels then provided the strongest constraints in the constrained minimal supersymmetric standard model, for instance (CMSSM), along with many other models of supersymmetry breaking mediation.
The higher collision energy at Run II allows for new MSSM parameter space to be explored.In the event of a discovery of the production of supersymmetric particles, one will want to first interpret their signals correctly, and then make inferences about the superymmetry breaking parameters.In order to do this with a greater accuracy, we must use higher orders in perturbation theory.Of particular interest is the connection between the supersymmetric masses of the squarks and gluinos and the experimental observables (functions of jet and missing momenta).The experimental observables may be used to infer pole (or kinematic) masses of the supersymmetric particles, which may then be connected via perturbation theory to the more fundamental supersymmetry breaking parameters in the Lagrangian.These could then be used to test the underlying supersymmetry breaking mediation mechanism [1,2].Conversely, if no significant signal for supersymmetry is to be found at the LHC, we shall want to interpret the excluded parameter space in terms of the fundamental supersymmetry breaking parameters.Again, this connection is sensitive (for identical reasons to the discovery case) to the order in perturbation theory which is used.In the most easily accessible channels at the LHC (squark/gluino production), it is useful therefore to use higher orders in the QCD gauge coupling, since this is the largest relevant expansion parameter.
State-of-the art publicly available NMSSM or MSSM spectrum generators such as ISAJET [3], FlexibleSUSY [4], NMSPEC [5], SUSPECT [6], SARAH [7], SPHENO [8], SUSEFLAV [9] or previous versions of SOFTSUSY [10], do not have the complete O(α 2 s /(16π 2 )) two-loop corrections to gluino and squark masses included, despite their being calculated and presented in the literature [11,12,13].Here, we describe their inclusion into the MSSM spectrum calculation in the popular SOFTSUSY program, making them publicly available for the first time.As we emphasised above, we expect them to be useful in increasing the accuracy of inference from data of supersymmetry breaking in the squark and gluino sectors.
The paper proceeds as follows: in the next section, the higher order terms that are included are briefly reviewed.We then provide an example of their effect on a line through CMSSM space.After a summary, the appendices contain technical information on how to compile and run SOFTSUSY including the higher order terms.

Higher Order Terms
Two-loop contributions to fermion pole masses in gauge theories were calculated in Ref. [11] and these results are specialised to compute the gluino pole mass.Depending upon the ratio m q/m g, a two-loop correction of up to several percent was found.On the other hand, Ref. [12] calculated the two-loop contributions to scalar masses from gauge theory and these results are specialised to the two-loop squark masses, where corrections up to about one percent were noted.In both cases, the results were obtained in the DR ′ renormalization scheme [14], consistent with the renormalization group equations used in softly broken supersymmetric models.The version of SOFTSUSY described here now contains these computations [11,12].A library for computing two-loop self-energy integrals, TSIL [15], is included within the SOFTSUSY distribution; there is no need to download it separately.
In addition, the gluino result has been improved by re-expanding the gluino self-energy function squared mass arguments around the gluino, squark, and top quark pole masses, as described in Ref. [13].This requires iteration to determine the gluino pole mass, and hence is slower than simply evaluating the self-energy functions with the lagrangian mass parameters, but the accuracy of the calculation is improved significantly [13].Typically, 4 to 6 iterations are required to reach the default tolerance of less than 1 part in 10 5 difference between successive iterations for the gluino pole mass.This iteration is not the bottleneck in the calculation, as we find that the CPU time spent on the squark masses is about twice that spent on the gluino mass, even with 6 iterations for the latter.We did not implement this improvement for the squarks since the effect is less significant, and the necessary iterations would be much slower.

Kinky Masses
In implementing these 2-loop pole masses, we encountered an interesting issue that does not seem to have been noted before, as far as we know.Consider the self-energy and pole mass of a particle Z that has a three-point coupling to particles X and Y. Let the tree-level squared masses of these particles by z, x, and y respectively.If z happens to be close to the threshold value ( √ x + √ y) 2 , then the computed 2-loop complex pole mass of Z will have a singularity proportional to where ϵ is infinitesimal and positive, and The reason for this can be understood from the sequence of Feynman diagrams shown in Figure 1.After reduction to basis integrals, the result of Figure 1(b) contains terms proportional to V(x, y, u, v) and B(x, y ′ ), in the notation of refs.[16,15].(Here the prime represents a derivative with respect to the corresponding argument, and the external momentum invariant is s = z.)Then, for example, one can evaluate: We have found that in general such singularities do not cancel within the fixed-order 2-loop pole mass calculation.This includes e.g. the pole mass of the gluino, where we checked that there are singularities in the pole mass as one varies the top mass and one of the stop masses very close to the 2-body decay threshold.Similarly, there are singularities in the 2-loop pole masses of the top squarks, in each case when the top mass and the gluino happen to be very close to the 2-body decay thresholds.We have also checked that this behavior occurs in a simple toy model involving only three massive scalar particles.This singular behavior might be quite surprising, because the pole mass is supposed to be an observable, and therefore ought to be free of divergences of any kind.The resolution is that, similar to problems with infrared divergences, the singularity is an artifact of truncating perturbation theory.The 3-loop diagram of Figure 1(c) will diverge like 1/(−δ) 3/2 , and similar diagrams of loop order L will diverge like 1/(−δ) L−3/2 .These contributions can presumably be resummed to give a result that is well-behaved as δ → 0, although proving that is beyond the scope of the present paper.
The above singularity behavior is not tied to the presence of massless gauge bosons.However, if massless gauge bosons are present, then there is another, less severe, type of singular threshold behavior, due to the Feynman diagram shown in Figure 2. In the notation of refs.[16,15], this arises due to the near-threshold dependence of the corresponding self-energy basis integral: Again, the singularity in the computed 2-loop pole mass is an artifact of the truncation of perturbation theory, and would presumably be removed by a resummation including all diagrams with 2, 3, 4, . . .massless gauge bosons exchanged between particles X and Y.
Note that it requires some bad luck to encounter any of these threshold singularity problems in practice, as there is no good reason why the tree-level masses should be tuned to the high precision necessary to make the problem numerically significant.Nevertheless, it could lead to a small1 "kink" in the computed pole mass if one performs a scan over models by varying the input masses.To avoid this possibility, in our code we test for small δ, and then replace the offending 2-loop self-energy basis integrals B(x, y ′ ) and V(x, y, u, v) and M(x, x, y, y, 0) by values interpolated from slightly larger and smaller values of the external momentum invariant.Explicitly, if δ < t (we choose t = 0.04 in the current version of SOFTSUSY, with the value controlled by the quantity THRESH TOL in src/supermodel/supermodel/sumo params.h), the integral in question I(s) is replaced by This provides a pole mass that varies continuously in scans of input masses near thresholds, and the difference between our value and the value that would be obtained by a proper resummation should be small.

Illustration of Results
We now illustrate the effect of the O(α 2 s /(16π 2 )) corrections to gluino and squark pole masses, taking the constrained minimal supersymmetric standard model (CMSSM) pattern of MSSM supersymmetry breaking Lagrangian terms.In the CMSSM, at the gauge unification scale ∼ O (10 16 ) GeV the gaugino masses are set equal to M 1/2 , the scalar masses are set to a universal flavour diagonal mass m 0 and the soft supersymmetry breaking trilinear scalar s /(16π 2 )) terms whereas m 1-loop is the pole mass calculated at one-loop order.In the right panel, ∆(σ) = σ 2-loop /σ 1-loop − 1, where σ 2-loop is the 13 TeV LHC production cross-section calculated using the O(α 2 s /(16π 2 )) terms in pole masses whereas σ 1 loop is the 13 TeV LHC production cross-section calculated with one-loop pole masses.The label 'gg' refers to gluino-gluino production, 'sg' to squark-gluino plus antisquark-gluino, 'ss' to squark-squark plus antisquark-antisquark, 'sb' to squark-antisquark, and 'tb' to stop-antistop production.The cross-sections in each case are obtained at NLL+NLO using NLL-fastv3.1-13TeV[17,18,19,20,18,21,22,23,24]. .couplings are set equal to a massive parameter A 0 times the relevant Yukawa coupling.Here, for illustration, we set M 1/2 = 1 TeV, A 0 = −2 TeV, the ratio of the MSSM's two Higgs vacuum expectation values tan β = 10 and the sign of the Higgs superpotential parameter term µ to be positive.We then allow m 0 to vary, and plot the relative difference caused by the new higher order terms in Fig. 3a.For this illustration, we have not included two-loop corrections to gauge and Yukawa couplings and three-loop renormalisation group equations for the superpotential parameters [25], although with them, the results are qualitatively similar.The overall message from the figure is that differences of percent level order in the pole masses of gluinos and squarks arise from the higher order corrections.In Fig. 3a, the plot does not extend to larger values of m 0 /M 1/2 , because there is no phenomenologically acceptable electroweak symmetry breaking there (the superpotential µ term becomes imaginary at the minimum of the potential, indicating a saddle point).The rough size of the two-loop correction is consistent with that estimated in the previous literature [11,12,13].Note that there are wiggles in the g, t1 , and t2 curves, near m 0 /M 1/2 = 1.35, 2.9, and 1.9, respectively.These are the remnants of the instabilities mentioned in section 3.1 after smoothing by the interpolation procedure described there, with t = 0.04.These correspond to the thresholds for the decays g → t t1 and t1 → t g and t2 → t g, respectively.Although the remaining wiggles are visible on the plot, they are small in absolute terms, and are of order the uncertainties due to higher-order corrections.
In Fig. 3b, we show the concomitant effect on the various next-to-leading order sparticle production cross-sections at the 13 TeV LHC by using NLL-Fast3.1 [17,18,19,20,18,21,22,23,24].NLL-Fast3.1 calculates at the nextto-leading order (NLO) in supersymmetric QCD, with next-to-leading logarithm (NLL) re-summation.However, it is based on fixed interpolation tables with only 3 significant digits, resulting in a visible jaggedness of the points in the figure.Nonetheless, we see that the change in the gluino mass due to next-to-next-to leading order (NNLO) effects leads to a large 15-25% reduction in the production cross-section.Other sparticle production cross-section modes shown decrease by more than 5%, thus accounting for these NNLO effects is important in reducing the theoretical uncertainties.Some of the curves terminate when m 0 /M 1/2 > 2.3 because NLL-Fast3.1 considers the production cross-sections to be too small to be relevant, and so returns zeroes.
In Fig. 4, we illustrate the renormalization scale (Q) dependences of the computed gluino and squark pole masses, by varying the scale Q at which the masses are calculated by a factor of 2 around M S US Y = √ mt 1 mt 2 .We see the In the case of the gluino, the solid line uses the re-expansion of squared mass arguments about the gluino and squark pole masses as described in [13]; this is the default used by SOFTSUSY for the gluino when the higher-order pole mass calculations are enabled.
expected reduced scale dependence going from 1-loop to 2-loop order.In the case of the gluino, we also see that the 2-loop calculation using a re-expansion around both the squark and gluino pole masses, as described in [13], displays less renormalization scale dependence than the case where only the gluino mass is re-expanded around its pole mass in the gluino pole mass calculation, or the case in which no re-expansion is used.This is consistent with the suggestion in [13] that the perturbation expansion is made more convergent by re-expanding around pole masses rather than around the running masses.For this reason, the re-expansion method is used as the default by SOFTSUSY when the higher-order pole masses are enabled.In the cases of the squark masses, the improvement of the renormalization scale dependence of the computed pole mass is less significant going from 1-loop to 2-loop order.We also note that in each case, the 2-loop correction is much larger than the scale dependence in the 1-loop result.Thus, as usual, the renormalization scale dependence gives only a lower bound, and not a reliable estimate, of the remaining theoretical error.

Summary and Conclusions
Two-loop O(α 2 s /(16π 2 )) are now included in the public release of SOFTSUSY, and so are available for use.We have demonstrated that along a typical line in CMSSM space, they are responsible for around O(1%) relative changes in the squark and gluino masses, which will change their various production cross-sections by around 5 − 25%.The largest 2-loop SUSYQCD pole mass correction is for the gluino, and it increases as the squarks become heavier than the gluino, reaching 2% when m 0 /M 1/2 = 4.These sources of theoretical error can now be taken into account and consequently reduced by using the new version of SOFTSUSY.In turn, the connection between measurements and various Lagrangian parameters (in particular, soft supersymmetry breaking sparticle mass parameters) is made more accurate.Thus, fits of the MSSM or NMSSM to cosmological and collider data will have a smaller associated theoretical error (as would studies of the unification of sparticle masses were there to be a discovery and subsequent measurement of sparticles).Other increases in MSSM or NMSSM mass prediction accuracy await future work: Higgs mass predictions, gauginos and sleptons, for instance.

Figure 2 :
Figure 2: Feynman diagram for the self-energy of a particle that couples to two lighter particles, each of which couples to a massless gauge boson.Near threshold, this diagram has a ln(−δ) singularity, where s = ( √ x + √ y) 2 (1 + δ), with x and y the squared masses of the internal particles.

Figure 4 :
Figure4: Scale dependences of gluino and squark pole mass predictions for the CMSSM with m 0 = M 1/2 = 1 TeV, A 0 = 0, tan β = 10 and µ > 0, comparing 1-loop and 2-loop approximations.The upper left panel is for the gluino, upper right panel for ũL , lower left panel for t1 , and lower right panel for t2 .In the case of the gluino, the solid line uses the re-expansion of squared mass arguments about the gluino and squark pole masses as described in[13]; this is the default used by SOFTSUSY for the gluino when the higher-order pole mass calculations are enabled.