Radiative corrections to stop-antistop annihilation into gluons and light quarks

We present the full one-loop SUSY-QCD corrections to stop-antistop annihilation into gluons and light quarks within the Minimal Supersymmetric Standard Model including Sommerfeld enhancement effects from the exchange of multiple gluons between the incoming particles. These corrections are important as stop (co)annihilation becomes the dominant contribution to the relic density for scenarios with a small mass difference between the neutralino and the stop which are less constrained by current LHC searches and consistent with the observation of a 125 GeV SM-like Higgs boson. We discuss important technical details of our one-loop, real emission, and resummation calculations where we pay particular attention to the cancellation of infrared divergences and the associated application of the dipole formalism for massive initial scalars. The corrections have been implemented in the dark matter precision tool DM@NLO which allows us to study numerically the impact of these corrections on the annihilation cross section. We find that for the chosen reference scenario the dominant correction comes from the Sommerfeld effect and that the pure NLO correction is below 3%. The inclusion of these radiative corrections is still large enough to decrease the relic density by more than 10% and shift the cosmologically preferred parameter region by a few GeV relative to the standard MicrOMEGAs result. Therefore, the inclusion of these corrections is mandatory if the experimental errors are taken as upper and lower bounds of the theory value.


INTRODUCTION
There is compelling evidence from astrophysical observations that there is a yet unknown type of matter called dark matter (DM) which does not interact electromagnetically but manifests itself through its gravitational effects on baryonic matter [1]. The most promising candidate for dark matter is a weakly interacting massive particle (WIMP) as it is consistent with structure formation due to its non-relativistic velocity and naturally leads via the freeze-out mechanism to the correct relic density of cold dark matter (CDM) Ω CDM h 2 = 0.120 ± 0.001 (1) as determined by the Planck satellite within the ΛCDM model [2]. The indicated uncertainty corresponds to the 1σ interval, and h stands for the present Hubble expansion rate H 0 in units of 100 km s −1 Mpc −1 .
As the Standard Model (SM) does not accommodate a suitable DM candidate there is the need for physics beyond the SM. A widely studied extension is the R-symmetric Minimal Supersymmetric Standard Model (MSSM) [3] as it contains not only an appropriate WIMP candidate in the form of the lightest neutralinoχ 0 1 , but also offers a solution to the hierarchy problem and allows * michael.klasen@uni-muenster.de † karol.kovarik@uni-muenster.de ‡ luca.wiggering@uni-muenster.de for the unification of gauge couplings at high energies. In order to make a theoretical prediction for the relic density of the neutralino under the assumption of the freeze-out scenario, one has to solve the Boltzmann equation for the DM number density n χ where n eq χ denotes the density in chemical equilibrium and H the Hubble rate [4,5]. Today's neutralino relic density is then given by where n 0 χ denotes the present value for the number density, m χ the DM mass and ρ c today's critical density. The number density equation in Eq. (2) is only an all order expression in the zero temperature limit since the phase space distribution functions of the SM particles are no longer exponentially suppressed by energy conservation for more than two particles in the initial or final state. This in principle forbids the usage of Maxwell-Boltzmann statistics and the neglect of Bose enhancement and Fermi blocking factors for 2 → 3 processes appearing at the oneloop level in the collision term. However, in Ref. [6] the additional thermal corrections where found to be suppressed by a factor T F /m χ 1 compared to zero temperature O(α s ) corrections with T F being the freeze-out temperature. The thermal corrections are therefore negligible at the current level of experimental precision of the dark matter relic density justifying the zero-temperature arXiv:2210.05260v2 [hep-ph] 10 Jan 2023 approach. The thermally averaged effective cross section involves a sum over all supersymmetric particles with odd R-parity where σ ij corresponds to the cross section for the annihilation of i and j into all possible SM particles. For the following analysis it is important to recall that the ratio n eq i /n eq χ is Boltzmann suppressed n eq i n eq with T being the temperature at time t. A direct consequence of Eq. (5) is that besides neutralino annihilation only those processes involving other particles from the odd sector in the initial-state with a small mass difference to the DM candidate can contribute significantly to σ eff v . Especially for large neutralino masses, the neutralino annihilation cross section alone is for many scenarios in the MSSM too small to be consistent with the measured relic density. Therefore, the neutralino cross section needs to be enhanced by some mechanism which could be colored (co)annihilation.
In this paper, we focus on the case where the mass of the lightest stop is very close to the one of the neutralino so that stop-antistop annihilation and stop pairannihilation become the dominant contribution to the effective cross section, and thus the relic density. This mass hierarchy is not an unnatural assumption since the tree-level mass of the lightest Higgs boson in the MSSM is bounded from above by m Z 0 | cos 2β| which requires large quantum corrections to be consistent with the observation of a SM-like 125 GeV Higgs boson [7,8]. The dominant contribution to the Higgs mass comes from the stop sector where a large trilinear coupling A t is needed in order for these corrections to be large enough, further indicating a large mass splitting between mt 1 and mt 1 [9]. The mass splitting is enhanced further through the fact that the off-diagonal entries in the sfermion mixing matrix are proportional to the associated masses of the SM partners, indicating a rather lightt 1 . The very small experimental uncertainty of the relic density in Eq. (1) requires the inclusion of radiative corrections to the annihilation cross section so that the theoretical precision matches the experimental one. However, public tools for the calculation of the relic density such as DarkSUSY [10] and MicrOMEGAs [11][12][13] only take into account the tree-level cross section with effective couplings that capture certain higher order effects.
The importance of higher-order SUSY-QCD corrections to the relic density has been shown for gaugino pairannihilation into quarks [14][15][16][17], gaugino-squark coannihilation into final states with a quark [18][19][20], squarkantisquark annihilation into electroweak final states [21], squark pair-annihilation into quarks [22] and stauantistau annihilation into heavy quarks [23]. Furthermore, the reduction of theoretical uncertainties from scheme and scale variations have been examined systematically [23,24]. Electroweak corrections to neutralino annihilation have been computed in [25][26][27]. It should be noted that the previous non-exhaustive list focuses only on one-loop corrections for relic density calculations. However, higher-order corrections in other contexts can also play an important role. The supersymmetric one-loop corrections in the strong coupling to the elastic neutralino-nucleon cross section relevant for direct detection were for example examined in Ref. [28], and one-loop EW corrections to Wino dark matter annihilation for indirect detection signals were computed in Ref. [29].
Based on these findings we present in this paper corrections of O(α s ) including Sommerfeld enhancement effects to the processest with the effectively massless quarks q ∈ {u, d, c, s}. These two processes are separate at tree level but have to be merged into one at NLO accuracy in order to obtain an infrared safe cross section. The paper is organized as follows: in Sec. II we present the color decomposed leading order cross section and discuss the phenomenological relevance of stop-antistop annihilation on the basis of a viable reference scenario. Sec. III covers details on the calculation of the virtual and real corrections, followed by the Sommerfeld resummation. In Sec. IV, we discuss the impact of the corrections on the corresponding cross section as well as the relic density for the chosen reference scenario. We conclude in Sec. V.

II. PHENOMENOLOGY OF SQUARK-ANTISQUARK ANNIHILATION
To prepare for the subsequent discussion of the higher order corrections and to clarify the notation, we start with the analytic computation of the tree-level cross section and discuss the phenomenology of the processes in Eq. (6) in the context of the neutralino relic density.

A. Leading order cross section
The Feynman diagrams for the leading order process are displayed in Fig. 1 along with the naming convention for momenta and other relevant indices. An important aspect of the processes we investigate is that both initial and final state particles are charged under SU(3) C . In order to be able to distinguish between attractive and repulsive color potentials in the context of the Coulomb corrections, it is necessary to decompose the tensor product representations under which the two incoming and   outgoing particles transform into their respective irreducible representations. The (s)quark-anti(s)quark system can be decomposed into a color octet and a color singlet whereas the decomposition of the two-gluon system reads For the decomposition of the tree-level scattering amplitudes into equivalent irreducible representations R that appear simultaneously in the initial as well as final state, the orthogonal and normalized multiplet basis elements c [R] spanning the invariant subspaces R from Ref. [30] can be used: as well as with C F = (N 2 c −1) /2Nc and N c = 3. Another important aspect in a non-Abelian theory is the treatment of internal and external polarization states. In order to include only the physical external gluon states in the transition probability, we consider two different computational approaches where we use the Feynman gauge for internal gluon lines within both possibilities. The first one is to explicitly sum only the transverse po-larizations with the help of the completeness relation T µ * which holds as an algebraic relation independently of the gauge fixing condition used for the internal propagators and where n is an arbitrary direction in momentum space that fulfills n · k = 0 and T (k) · n = 0. For some n with n 2 = 0 this is also referred to as the light-cone gauge. As there appear only two external gluons in the tree-level process, it is instructive to choose n as the momentum of the respective other gluon. The second possibility is to use −g µν as polarization sum and subtract the longitudinal polarizations by using ghosts. To arrive at the corresponding expression, we derive the two Slavnov-Taylor identities from the invariance of a general n-point function in SUSY-QCD under Becchi-Rouet-Stora (BRS) transformations [31,32]. Consequently, Eq. (13) allows to replace the longitudinal polarizations corresponding to all the terms proportional to k 1 and k 2 in Eq. (12) with ghost amplitudes. This gives for the squared matrix element summed over final-state polarizations The fermion spin sum for the quark-antiquark final state is performed in the usual way. After averaging (summing) over initial-(final-) state colors and performing the remaining phase-space integration, we obtain for the color-decomposed tree-level cross sections describing the annihilation into two gluons the expressions (σv) Tree gg,[8 S ] = 5 2 (σv) Tree gg, [1] (σv) Tree gg, with ρ = 4m 2 q/s and β = √ 1 − ρ where v = 2β corresponds to the relative velocity of the incoming squarkantisquark pair in the c.m. system and s = (p 1 + p 2 ) 2 to the squared c.m. energy. Only one color channel contributes to the annihilation into a massless quarkantiquark pair giving the cross section .
As we have to combine both processes at NLO, we define already at tree-level (σv) Tree = (σv) Tree gg + N f (σv) Tree qq (16) where N f = 4 corresponds to the number of effectively massless quark flavors.

B. Reference scenario and numerical discussion
To illustrate the importance of stop annihilation into gluons, we introduce the reference scenario given in Tab. I which has been found by performing a random scan in the MSSM with 19 free parameters considering the most important experimental constraints from searches for supersymmetry. For this scan and throughout our analysis SoftSUSY 4.1.9 [33][34][35][36] is used for the generation of the mass spectrum and mixing parameters with the option of including three-loop corrections to the mass of the CP-even Higgs boson h 0 provided by Himalaya 1.0 [37,38] turned on. Only those points that obey the Higgs mass limit 123 GeV < m h 0 < 127 GeV, feature the neutralino as lightest supersymmetric particle (LSP) and a stop as next-to-lightest supersymmetric particle (NLSP) are taken into account. We use SModelS 2.2.0 [39][40][41][42][43] and SUSY-AI [44] to exclude points that have been ruled out by LHC searches for supersymmetry. The consistency of the Higgs sector with measurements from LEP, Tevatron and the LHC is additionally checked with HiggsBounds 5.5.0 [45] and HiggsSignals 2.3.0 [46]. The module in MicrOMEGAs-5.2.13 [13] is used to check against constraints from dark matter direct detection experiments. However, unless stated otherwise we use throughout this paper MicrOMEGAs 2.4.1 [11,12] with the standard CalcHEP implementation of the MSSM for the computation of the relic density and the contributions of different (co)annihilation channels.
The latter are shown in Tab. II for the chosen reference scenario.
The largest contribution comes with 47 % from stop-antistop annihilation into gluons followed in decreasing order by stop pair-annihilation into top quarks and neutralino-stop coannihilation into a gluon and a top quark which have been previously analyzed in [22] and [20], respectively. In total, DM@NLO provides full one-loop SUSY-QCD corrections to 77 % of the effective cross section in Eq. (4).
The scenario features a bino-like neutralino which is not surprising as large wino and higgsino components would lead to other gauginos being the NLSP and the mass difference between the neutralino and the lightest stop is approximately 11 GeV. The gluino and slepton sector are chosen to be much heavier than the stop sector to ensure that they do not influence the phenomenology discussed here. In Fig. 2, the relative contributions of the three most important channels to the relic density are displayed in the M 1 -Mt R mass plane in different shades of green. We choose these two parameters as the lightest neutralino is mostly bino-like and its mass is consequently predominately given by M 1 . The Mt R parameter enters the tree-level expression of thet 1 mass and therefore these two parameters correspond to a scan in the mχ0 1 -mt 1 mass plane which in turn allows to investigate the dependence of the relic density on the LSP-NSLP mass difference. For larger mass splittings between the lightest neutralino and the stop coannihilation becomes the dominant channel whereas for small mass splittings  1 and mt 2 are in this case the DR tree-level masses, the associated pole masses of relevant particles, the bino contribution Z 11 toχ 0 1 and the neutralino relic density. All dimensionful quantities are given in GeV.  annihilation of stops is the dominant contribution. In addition, the region where the neutralino accounts for the whole dark matter content in the universe and lies within the 2σ range of the experimental value is marked in orange. This region follows an almost straight line parallel to the boundary where the neutralino is no longer the LSP.
With the knowledge that stop annihilation into gluons is important for large regions around the reference sce-nario, we turn now to the numerical comparison between our leading order cross sections for the two processes in Eq. (6) and the ones from MicrOMEGAs 2.4.1 which are all shown in Fig. 3. As a reminder that the values of of the cross section impacts the relic density only in a limited energy range, the Boltzmann distribution which is involved in the computation of the thermally averaged cross section at freeze-out temperature is shown in gray in arbitrary units. One observes that our result is about 6 % larger for both processes which has two reasons. Firstly, we set the renormalization scale which enters at treelevel only through the strong coupling to µ R = Q SUSY whereas MicrOMEGAs 2.4.1 sets the scale to twice the dark matter mass µ MO = 2mχ0 1 which is larger than µ R for the investigated scenario and therefore corresponds to a smaller strong coupling. Our choice for µ R is motivated by the fact that the besides the masses of the virtual particles in the loop, the process contains only two important scales: the mass of the lightest stop and the collisional energy s. Since most annihilations take place between s = 4m 2 t1 and the peak of the velocity distribution at s ∼ (3 TeV) 2 , Q SUSY is a suitable choice for the renormalization scale to avoid large logarithms. Secondly, MicrOMEGAs 2.4.1 calculates the running of α s in the MS-scheme using the three-loop formula in Ref. [47] with six active flavors and the SM particle content only whereas DM@NLO uses the four-loop formula from Ref. [48] in the DR-scheme [49] with six active flavors and contributions from the complete MSSM mass spectrum [50]. Considering only these two differences in the computation, the ratio should be identical for both processes, but this is not the case as MicrOMEGAs also takes into account electroweak contributions with a photon or a Z 0 propagator for the process with a quark-antiquark pair in the final state. The corresponding electroweak diagrams are not included in our calculation since the process with massless quarks is numerically insignificant for the relic density as well as the tree-level cross section compared to the annihilation into gluons as visible in Tab. II and Fig. 3 and was only added for consistency to achieve an infrared finite result.
Through comparison of the different color contributions to the combined leading order cross section depicted in Fig. 3 with the partial wave expansion of a general velocity-weighted annihilation cross section σv, it becomes apparent that the singlet and symmetric octet contributions to the cross section with two external gluons are dominated by the S-wave component s 0 Channel Contributioñ DM@NLO total [20,22] 77 % since they remain almost constant in v, whereas the antisymmetric octet part of the same process and the octet contribution to the quark-antiquark process take an inferior role and are suppressed at threshold corresponding to the S-wave and P -wave component s 1 .

III. COMPUTATIONAL DETAILS OF THE RADIATIVE CORRECTIONS
In this section, we discuss the technical details of our SUSY-QCD corrections at O(α s ) as well as the Sommerfeld enhancement. The NLO cross section with the NLO correction

A. Virtual corrections and renormalization
The virtual amplitudes consist of propagator (selfenergy), vertex and box corrections. Naively one might assume that the box corrections for the process with two final-state gluons are independent and UV finite on their own. However, they turn out to be UV divergent and fall under the renormalization of the four-squark-gluon vertex. All corresponding Feynman diagrams are shown in Figs. (4) to (12). We subtract the longitudinal gluon polarizations again through ghosts, i.e. the interference of the tree-level matrix element with the virtual amplitudes for the process with two gluons in the final state summed over the final-state polarizations can be written as where some of the ghost corrections making up the ghost amplitudes S NLO i (i = 1, 2) are shown in Figs. (9) and (10). These diagrams are regulated dimensionally in D = 4 − 2ε dimensions within the supersymmetry preserving four-dimensional helicity scheme [55][56][57] so that UV and IR divergences appear as poles of the form ε −1 and ε −2 . The standard Passarino-Veltman reduction [58,59] is used to express the one-loop amplitudes in terms of the well-known scalar integrals A 0 , B 0 , C 0 , D 0 [60][61][62]. The γ 5 -matrix which enters through the squarkquark-gluino coupling is treated in the naive scheme, i.e.  we assume that γ 5 still anti-commutes with all γ-matrices in D dimensions. The Levi-Civita symbols that occur then through traces of γ 5 with four or more γ-matrices during the evaluation of diagrams with top quarks as virtual particles are directly set to zero since they vanish anyway when being contracted with the external momenta. The UV divergences that appear in the virtual corrections are removed through the renormalization of fields, masses and the strong coupling. Within our calculation, a hybrid on-shell/DR renormalization scheme is employed where A t , A b , mt 1 , mb 1 , mb 2 along with the heavy quark masses m t , m b are treated as independent input parameters so that the mixing angles θt 1 , θt 2 and the mass of the heavier stop mt 2 depend on their definition. The trilinear couplings of the third generation, the bottom quark mass and the strong coupling are renormalized in the DR scheme while the on-shell scheme has been chosen for the top mass and the input squark masses. This particular scheme resembles the RS2 scheme introduced in Ref. [63] and was found to be robust over large regions of the parameter space for (co)annihilations involving stops in a series of previous analyses [19,20]. Since the renormalization of the gluon and the squark sector as well as the treatment of the bottom mass and the strong coupling have already been discussed in detail in the context of other processes [17,19,20], we will only cover aspects which are new to this calculation in the following such as the renormalization of ghosts and massless quarks.

Ghost wave-function renormalization
As ghost and anti-ghost share the same self-energy they can be renormalized with the same wave function renormalization constant Z c . The renormalized fields are then defined as where we need δZ c only up to O(α s ) which leads to the expansion Since the gluon is renormalized in the on-shell scheme, the same scheme is chosen for the ghost. That is, the ghost renormalization constant is obtained by requiring that the ghost Green's function has a unit residue even up to the one-loop level whereΠ denotes the derivative of the ghost self-energy whose only contribution is depicted in Fig. 13. The constant δZ c contains UV and IR divergent parts which read explicitly

Renormalization of the massless quarks
For the renormalization of massless quarks, we introduce the quark wave-function renormalization constants Z L/R q for each chirality state The renormalization constants are determined in the onshell scheme which requires the renormalized quark twopoint Green's function to have a unit residue. This condition results in the expression where the function Π L/R q (q 2 ) appears in the decomposition of the quark self-energy Π q (p) = / p P L Π L q (p 2 ) + P R Π R q (p 2 ) + Π S,L q (p 2 )P L + Π S,R q (p 2 )P R (30) whose two contributing Feynman diagrams are shown in Fig. 14. The resulting constants contain the UV and IR divergent parts where the superscripts indicating the left/right-handed chirality states are dropped here for simplicity.

B. Real corrections
The infrared divergences in the virtual corrections are compensated by including the real emission processes with q ∈ {u, d, c, s} being an effectively massless quark and where the initial squarks carry the same labels as in Fig. 1. The corresponding Feynman diagrams are shown in Figs. (15a) and (15b) where the momenta of the gluons in the first process have to be read from top to bottom starting with k 1 . As in the tree-level calculation, we use −g µν for the gluon polarization sum and subtract the longitudinal polarizations with ghosts as asymptotic states. In order to arrive at the corresponding expression, we proceed as sketched in Sec. II A by deriving the following two sets of Ward identities from BRS-invariance and where M µνρ as an expression for the squared matrix element summed over the physical final state polarizations. Eq. (37) follows from an explicit calculation with the help of Feynman rules. The final expression in Eq. (38) obeys the same structure as the one from the 2 → 2 calculation. The ghost processes are only squared with themselves and then subtracted from the matrix element squared of the actual process.
We now turn to the discussion of the treatment of infrared divergences. To make the integration over the three-particle phase space numerically accessible and to combine the real and virtual corrections to get an infrared safe cross section, we rely on the dipole subtraction methodà la Catani-Seymour [64] which has recently been extended to massive initial states in the context of dark matter calculations [65]. This method is based on the introduction of an auxiliary differential cross section dσ A which cancels the soft and collinear divergences of the differential real emission cross section pointwise but can be integrated analytically at the same time over the oneparticle phase space responsible for the soft or collinear divergence. That is, the NLO correction takes the form According to the dipole factorization formula, the auxiliary squared matrix element related to dσ A for the process with three gluons in the final state consists of 27 dipoles where the subscripts of the momenta in Eq. (33) and Eq. (34) are used to label the particles. For the precise definition of the dipoles and the underlying splitting kernels we refer to Ref. [65]. For the process containing light quarks we obtain the 15 dipoles For the explicit construction of the insertion operator which cancels the infrared divergences on the virtual side, we refer again to Ref. [65] due to the large number of terms coming from the non-factorizable color and spin structures.

C. Sommerfeld enhancement
We have discussed the fixed-order NLO corrections in the previous two subsections. However, for the nonrelativistic regime, as it is typical during freeze-out, there are also important contributions to the relic density from the exchange of n potential gluons between the incoming stop and antistop giving a correction factor proportional to (α s /v) n . This is the well-known Sommerfeld enhancement [66] of higher-order terms which can spoil the perturbativity of the cross section when the relative velocity is of the order of the strong coupling, and therefore these terms need to be resummed to all orders in perturbation theory. The fact that the tree-level cross section is dominated by S-wave annihilation as discussed in Sec. II B and visible in Fig. 3, allows to compute the Sommerfeld enhanced cross section (σv) Som = S 0, [8] +N f (σv) Tree qq, [8] + S 0, [1] (σv) Tree gg, [1] (42) by multiplying the leading contribution with the Sommerfeld factor whose computation follows the standard framework of non-relativstic QCD (NRQCD) described in Refs. [67,68]. The Green's function G [R] ( r = 0, √ s + iΓt 1 ) is defined as solution of the Schrödinger equation evaluated at the origin where is the Hamiltonian of the quasi-stoponium. The corresponding Coulomb QCD potential receives important contributions from gluon and fermion loops and reads at NLO in momentum spacẽ with the color factors and the constants where we work with n f = 5. The analytic solution for the Green's function at the origin at NLO accuracy is where the LO and NLO contributions are Here, the constants are defined through the non-relativistic velocity of the incoming particles and ψ (n) = ψ (n) (1 − κ) is the n-th derivative of ψ(z) = γ E + d/dz ln Γ(z) with the argument (1 − κ). For the computation of the Sommerfeld factor, we also need the free Greens's function We address now the choice for the Coulomb scale µ C at which the strong coupling in the QCD potential is evaluated. Following Ref. [69], we set where 2mt 1 v s is motivated by the typical momentum transfer mediated by the potential gluons. The Bohr scale µ B corresponds to twice the inverse Bohr radius r B and is obtained by iteratively solving the equation For the scenario in Tab. I, the Bohr scale takes the value µ B = 204 GeV and the associated value for the strong coupling in the MS-scheme with 6 active quark flavors is α s (µ B ) = 0.1058. As a single gluon exchange is already included in our fixed-order NLO calculation (see Fig. 7 and Fig. 8), we have to match it to the Sommerfeld enhanced cross section in order to avoid double counting. This is achieved by taking only the terms of O α 2 s in Eq. (43) into account giving the full cross section (σv) Full .
As described in Ref. [22], it is also possible to subtract the velocity-enhanced part from the fixed order calculation in order to obtain the "pure" NLO cross section which gives with the relativistic relative velocity

IV. NUMERICAL RESULTS
In this section, we discuss the impact of the corrections on the stop-antistop annihilation cross section and the corresponding impact on the theoretical uncertainty deduced from scale variations. Then, we study the impact of the full correction on the relic density for stop-antistop annihilation alone as well in conjunction with the other two important processes shown in Tab. II.

A. Annihilation cross section and its theoretical uncertainty
In Fig. 16a, we show the stop-antistop annihilation cross cross section as a function of the CM momentum p cm for the parameter point defined in Tab. I. More precisely, we show the cross section at tree-level as provided by DM@NLO (black dashed line) and by MicrOMEGAs  FIG. 14: One-loop contributions to the quark self-energy.
2.4.1 (dotted orange line), including the NLO corrections (green solid line) and the full cross section with the Sommerfeld enhancement effect (red solid line). In addition, we show the pure Sommerfeld enhanced cross section (blue dashed line) and the "pure" NLO cross section without the velocity-enhanced part (purple solid line). For small relative velocities the Coulomb corrections from the exchange of multiple gluons between the incoming particles dominate the full corrected annihilation cross section. As discussed in Sec. III C, the effect of the Coulomb corrections depends on the quadratic Casimir of the representation under which the incoming particles transform. The singlet feels an attractive force whereas the squark and antisquark transforming under an eight dimensional representation are repelled from each other. In this case, the Coulomb corrections increase the annihilation probability so that the full corrected cross section becomes larger than 100 % of the tree-level cross section for CM momenta below 88 GeV even though the LO cross section is dominated by the symmetric octet contribution which is due to the color suppression given by 1 /2Nc in the Sommerfeld factor for the eight dimensional representation. For vanishing relative velocities, the enhanced cross section even diverges and approaches the well-known Coulomb singularity which could be cured by taking the formation of bound states into account properly. However, as the Boltzmann distribution almost vanishes for momenta around p cm = 0, such effects are heavily suppressed. In contrast, the "pure" NLO correction without any enhancement corresponds to an improvement of less than ±3 % of the LO cross section such that the full corrected cross section is in very good approximation given by the pure Sommerfeld enhancement, i.e. (σv) Full ≈ (σv) Som .
The other two processes which we include in our analysis and are important in the region around the reference scenario, namelyt 1t1 → tt andχ 0 1t1 → tg, have been investigated in the context of DM@NLO in Refs. [20,22]. In contrast to the two original publications, we do not use the phase space slicing method for the real corrections in this paper but the dipole subtraction method. The implementation of the dipole approach for the two processes and the comparison between both methods were the subjects of Ref. which is used by DM@NLO, the difference between the two tree-level cross sections in Fig. 16c is due to the differences in the strong coupling as discussed in the context of the LO cross section oft 1t * 1 → gg. In the case of stop pair-annihilation, the NLO corrections cause a positive shift of about 10 % for large p cm around 600 GeV compared to the tree-level cross section whereas the correction becomes large and negative for CM momenta less than 287 GeV. For CM momenta below 50 GeV the total cross section becomes negative which is unphysical but we make in the following the assumption that this momentum region is irrelevant for the computation of the relic density due to an almost vanishing Boltzmann distribution for such low velocities. Furthermore, this unphysical behavior has already been extensively discussed in Ref. [22]. In the case of neutralino-stop coannihilation the correction is stable around 19 % for all relevant CM momenta.
We continue with the analysis of the theoretical uncertainties of the stop-antistop annihilation cross section from variations of the Coulomb and renormalization scale where we identify the central scales with the ones used in the previous discussion, i.e. µ central R = Q SUSY and µ central C = max 2mt 1 v s , µ B . In Fig. 17, we vary µ R and µ C by factors of two and show the associated values of the annihilation cross section at tree-level (blue), at NLO (green) including the Coulomb corrections (red) as well as the pure Sommerfeld enhanced cross section (purple) normalized to the corresponding cross section obtained at the central scale(s). In conjunction, the LO and NLO cross section as function of the renormalization scale for three different CM momenta are shown in Fig. 18. Within the chosen renormalization scheme, the scale dependence enters the tree-level cross section only through the strong coupling and we estimate the theoretical uncertainty to about ±5.5 %. For large CM momenta (p cm ≈ 900 GeV) the NLO correction lies within the LO uncertainty and the theoretical uncertainty is reduced to below 1 %. For intermediate energies (p cm ≈ 300 GeV) the NLO correction is no longer contained in the LO uncertainty but the uncertainty is still reduced to about ±1.5 % by including the higher-order corrections. For very small relative velocities (p cm ≈ 10 GeV) the cross section becomes non-perturbative and the NLO uncertainty is larger than the LO one. However, by including the Coulomb corrections the upper uncertainty bound for small energies is halved whereas the lower uncertainty bound increases and we have only a reduction for v → 0. As the full corrected cross section is in very good approximation given by the Sommerfeld enhancement only, we expect the same for the associated uncertainty which (a) Graphs with three gluons in the final state that are associated with the amplitude M3.
(b) Graphs with light quarks in the final state.
(c) Graphs with ghosts associated with the amplitude S1. The amplitude S2 is obtained by reversing the ghost flow.
(d) Graphs with ghosts associated with the amplitude S3. The amplitude S4 is obtained by reversing the ghost flow.
(e) Graphs with ghosts associated with the amplitude S5. The amplitude S6 is obtained by reversing the ghost flow. turns out to be the case. We note at this point that the kink in the uncertainty band of (σv) Full and (σv) Som comes from the transition from the Bohr scale to the scale of the typical momentum exchange 2mt 1 v s .

B. Impact on the relic density
At last, we investigate the impact of our radiative corrections on the neutralino relic density Ωχ0 1 h 2 by including all three processes from Tab. II which are important in a region around the chosen reference scenario and are 100 300 500 700 900 available in DM@NLO as well as for the process which is subject of this paper only. This means that the integration of the Boltzmann equation in Eq. (2) is still performed by MicrOMEGAs 2.4.1 but the cross sections are replaced by the ones implemented in DM@NLO for the specified cases and still obtained from CalcHEP for the remaining ones. Similar to Sec. II, we study the impact on the relic density in the plane spanned by M 1 and Mt R which is shown for both cases in Fig. 19.
As before, the region which is compatible up to two sigma with the Planck limit is shown in orange for the values obtained with MicrOMEGAs 2.4.1, in blue for the tree-level values from DM@NLO and in gray for the radiative corrections. In addition, the same results are presented in Fig. 20 projected into the plane of the physical neutralino and stop mass where one should highlight that this variation only comes from the scan over the parameters M 1 and Mt R whereas all other parameters in Tab. I remain fixed. The small difference between the tree-level results is again mainly due to the differences in the strong coupling.
In both cases, the favored parameter region consistent with the Planck limit is shifted towards larger stop masses for a fixed neutralino mass to compensate the increased effective annihilation cross section where this shift exceeds the experimental uncertainty. However, if we only include the radiative corrections for stop-antistop annihilation the cosmologically favored stop mass is increased by about 6.1 GeV compared to the MicrOMEGAs result whereas the additional inclusion of the higherorder corrections to the processesχ 0 1t1 → tg andt 1t1 → tt reduces this shift to about 4.3 GeV. This is due to the large negative NLO corrections for small p cm that occur for stop pair-annihilation.

V. CONCLUSION
The annihilation of colored particles which are close in mass to the dark matter candidate is an important mechanism to allow for higher dark matter masses while still being able to explain the measured relic density. In Cross section ratios  the MSSM, a theoretically well motivated candidate for such annihilation processes is the lightest stop. Based on previous analyses which show that the inclusion of higher-order corrections to the relic density exceeds the experimental uncertainty of the dark matter content in the universe, we have presented in this paper NLO SUSY-QCD corrections to stop-antistop annihilation into gluons and light quarks including QCD Coulomb corrections of O α 2 s . The two processest 1t * 1 → gg andt 1t * 1 → qq with q being an effectively massless quark are combined in our analysis since we found within our calculation that these two processes can not be treated separately at NLO accuracy in order to guarantee a well-defined and infrared safe cross section. In order to study the impact of such corrections on the annihilation cross section itself and the relic density, we have performed a random scan in the phenomenological MSSM with 19 free parameters to select a reference scenario that is consistent with the current most import experimental constraints and contains a stop with almost the same mass as the neutralino. The numerical analysis showed that the resummed cross section matched to the fixed-order NLO calculation is in very good approximation given by the Sommerfeld enhanced cross section only which can in turn be used to significantly speed up relic density scans while capturing the majority of the NLO corrections. We are confident that this result extends to simplified dark matter models containing a colored scalar similar to the MSSM as those proposed for LHC searches in Ref. [70]. In addition, we observed that the inclusion of the NLO corrections reduces the dependence of the cross section on the renormalization scale in the perturbative regime from ±5.5 % to below ±2 %. Finally, we found with respect to the impact on the relic density that the corrections to stop-antistop annihilation only can shift the cosmologically favored parameter region by a few GeV and they are therefore larger than the current experimental uncertainty. However, through the additional inclusion of the NLO SUSY-QCD corrections toχ 0 1t1 → tg andt 1t1 → tt this shift is reduced by about 30 % due to a large negative correction for the stop pair-annihilation. As in our previous studies, we conclude that the identification of parameter regions consistent with the measured relic density at the current level of precision requires the inclusion of NLO and Coulomb corrections including those covered in this work.