Squark-pair annihilation into quarks at next-to-leading order

The Minimal Supersymmetric Standard Model (MSSM) is under intense scrutiny at the LHC and in dark matter searches. Interestingly, scenarios with light squarks of the third generation remain not only viable, but also well motivated by the observed Standard-Model-like Higgs boson mass and dark matter relic density. The latter often requires important contributions from squark pair annihilation. Following up on previous work, we present in this paper a precision analysis of squark pair annihilation into quarks at next-to-leading order of QCD including Sommerfeld enhancement effects. We discuss all technical details of our one-loop, real emission and resummation calculations, their implementation in the precision tool DM@NLO, as well as the numerical impact on the annihilation cross section and cosmological relic density in phenomenological MSSM scenarios respecting in particular current LHC constraints. We demonstrate that including these radiative corrections leads to substantial shifts in the preferred parameter regions by up to 20 GeV.


I. INTRODUCTION
Strong evidence for the existence of dark matter, along with the fact that neutrinos are massive, is a compelling sign of the need for physics beyond the Standard Model. Even though dark matter still evades direct detection in Earthbased experiments such as LUX [1] and XENON1T [2], there is overwhelming evidence from cosmological data such as the cosmic microwave background that dark matter exists in the Universe. Moreover, the relic density of cold dark matter (CDM) has been determined to the unprecedented precision of Ω CDM h 2 ¼ 0.1200 AE 0.0012 ð1:1Þ as measured by the Planck satellite and interpreted within the ΛCDM cosmological model [3]. The indicated uncertainty corresponds to the 68% confidence level, and h stands for the present Hubble expansion rate H 0 in units of 100 kms −1 Mpc −1 .
The Standard Model does not contain any suitable candidate for dark matter with the required properties. The leading candidate therefore remains a weakly interacting massive particle (WIMP), which leads to the correct relic density via the freeze-out mechanism. However, alternative candidates and mechanisms do exist, e.g., in the form of the freeze-in mechanism [4][5][6][7]. The Standard Model is therefore extended to include new particles which provide the required dark matter candidate. The new particles are usually protected from decaying by introducing an ad hoc Z 2 -symmetry, where all new particles are Z 2 -odd and the Standard Model particles are Z 2 -even. One such model, which was actually not introduced to address the existence of dark matter, is the Minimal Supersymmetric Standard Model (MSSM), where the conserved Z 2 -symmetry is the R-parity. In most MSSM scenarios, the lightest supersymmetric particle (LSP) is the lightest neutralinoχ 0 1 , which is stable and interacts only weakly. The lightest neutralino is an extremely well-studied candidate for cold dark matter.
The theory prediction for its relic abundance, Ω˜χ0 1 h 2 , is related to the number density n χ of the neutralino, which can be computed by solving the Boltzmann equation [8][9][10] dn χ dt ¼ −3Hn χ − hσ ann vi½n 2 χ − ðn eq χ Þ 2 ; ð1:2Þ where H denotes the (time-dependent) Hubble parameter, n eq χ the number density in thermal equilibrium, and v the Møller velocity of the annihilating particles [8]. All specifics about the interaction of dark matter with other particles in the chosen particle physics model are contained in the annihilation cross section σ ann , which accounts for all possible annihilation and co-annihilation processes. Its thermal average can be expressed as T being the temperature. From this last equation, it becomes obvious that annihilations involving particles other than the lightest neutralino are suppressed if these particles are heavy compared to the neutralino. On the other hand, coannihilations with the next-to-lightest supersymmetric particle (NLSP) will be important or even dominant if the mass difference is rather small. Typical examples in the MSSM are co-annihilations of the neutralino with a scalar top quark or a scalar tau lepton. For smaller mass differences between the LSP and NLSP, even pair annihilations of the next-to-lightest particle contribute in a sizable manner, and can even become dominant in the total annihilation cross section. In case there are more than two almost mass degenerate particles, (co-)annihilations between all particles have to be taken into account.
In the present paper, we focus on the case where the masses of one or two squarks, the lightest stop and/or the lightest sbottom, are close to the neutralino mass. The case of a light scalar top quark is very well motivated. Scenarios with light scalar tops satisfy the experimental constraints from LHC searches and can also contribute to a successful prediction of the mass of the lightest Higgs boson in the MSSM [11].
The relic density of dark matter in scenarios with a light stop which is almost mass degenerate with the lightest neutralino is very sensitive to the mass difference of the two particles. Any small shift in the predicted relic density can cause a large shift of the parameter region where the relic density is compatible with the experimental limits given by Eq. (1.1). In this analysis we focus on next-to-leading supersymmetric QCD (SUSY-QCD) corrections to the corresponding annihilation and co-annihilation cross sections in scenarios with a light scalar quark. These corrections have the potential to significantly modify the annihilation cross section and thereby also the relic density.
Including radiative corrections to the total annihilation cross section not only shifts the parameter regions corresponding to the correct relic density, but it also reduces the theoretical uncertainty of the relic density prediction. The theoretical uncertainty from scale and scheme variations on the annihilation cross section and the neutralino relic density has been evaluated for specific subclasses of processes in Ref. [23].
After the work presented in Refs. [12][13][14][15][17][18][19], with the present paper, we make a first step towards completing the missing processes sensitive to radiative corrections of order α s . More precisely, we present such corrections for squark-squark annihilation into quark-quark pairs. The discussion of squark-antisquark annihilation into quarks and gluons is left for forthcoming publications.
In Sec. II, we start by discussing the phenomenological importance of the processes under consideration in this work. We also present two reference scenarios featuring important contributions of the processes of our interest. In Sec. III, we then detail the analytical calculation of the radiative corrections. We discuss, in particular, points that are beyond the discussion presented in Refs. [12][13][14][15][17][18][19] and analyze the impact that the radiative corrections have on the corresponding cross sections. The impact of the corrections on the relic density in the two reference scenarios is presented in Sec. IV. Our conclusions are given in Sec. V.

II. PHENOMENOLOGY OF SQUARK ANNIHILATION
The analysis presented in this paper concentrates on the contributions from squark-pair annihilation to the total annihilation cross section σ ann of neutralino dark matter. We investigate scenarios in the phenomenological MSSM (pMSSM), where the processes play an important role. Supersymmetry and the MSSM, in particular, have been extensively tested by searches at the Large Hadron Collider (LHC) and at experiments aiming at the detection of direct signals from elastic collisions of dark matter with heavy nuclei such as XENON1T. In order to take into account the most important experimental constraints from the searches for supersymmetry, we use the results of an analysis performed by the ATLAS Collaboration in the light of recent searches at the LHC [24]. 1 The ATLAS analysis is performed in the pMSSM with 19 parameters (defined at the 1 SUSY scale) and is based on a sample of 5 × 10 8 parameter points. Applying constraints from ATLAS SUSY searches, electroweak precision observables such as Δρ and ðg − 2Þ μ , flavor observables such as b → sγ, and requiring that the neutralino is the LSP and a dark matter candidate with the relic density less than 2 0.1208 (for details on further constraints see Ref. [24]) leads to a subset of about 300 000 viable points. We have analyzed this subset in order to examine in which regions of parameter space the above processes contribute significantly.
In order for the contribution from the annihilation of third-generation squarks to the total dark matter cross section to be significant, one (or more) scalar quarks have to be almost mass degenerate with the lightest neutralino. This is not an unnatural requirement because a light scalar top quark is necessary to explain the measured mass of the Standard Model Higgs boson within the MSSM. Moreover, scenarios where a scalar top is almost mass degenerate with the lightest neutralino are quite frequent, as the mass degeneracy gives rise to different topologies in collider searches, making their testing more challenging and their exclusion less likely. Another aspect of scenarios where scalar quarks are the NLSP is that the lightest neutralino is mostly bino-like. Higgsino-like and wino-like lightest neutralinos mostly lead to scenarios with other gauginos being the NLSP. A consequence of the lightest neutralino being bino-like is that the annihilation of dark matter is typically not efficient enough for the relic density to reach the value determined by the Planck Collaboration given in Eq. (1.1). Therefore, in scenarios with bino-like neutralino dark matter, some enhancement mechanism is needed for them to be consistent with the relic density measurements. As we will discuss below, in the scenarios analyzed here, the enhancement comes from the presence of LSP-NLSP co-annihilations as well as NLSP annihilations.

A. Reference scenarios
As mentioned above, the numerical part of the present study will be based on two reference scenarios inspired by the findings presented in Ref. [24]. More precisely, we will focus on two pMSSM scenarios, whose most relevant softbreaking parameters and particle masses are presented in Table I. It is to be noted that, although the input soft mass parameters of the two scenarios are identical to those of two actual scenarios given in Ref. [24], the resulting physical masses slightly differ from those associated with the ATLAS study due to the fact that we are using a different computational setup. The actual shift in the physical masses is small so that all experimental constraints are still satisfied and the phenomenology is not altered.
Both scenarios feature bino-like neutralinos, the bino mass parameter M 1 being smaller than the wino and Higgsino mass parameters M 2 and jμj. The key parameters of the third-generation squarks of our interest are the "lefthanded" stop and sbottom mass parameter Mq L , and the "right-handed" stop and sbottom mass parameters Mt R and Mb R . In both scenarios, squarks of the first and second generations, the sleptons, and other electroweak gauginos are heavier such that they do not influence the phenomenology discussed here.
In our setup, starting from the soft-breaking terms defined at the scale Q SUSY indicated in Table I, we obtain the physical mass spectrum using the spectrum generator SPheno 3.3.3 [27,28]. The mass spectrum is then handed over to micrOMEGAS 2.4.1 [29,30] making use of the SUSY Les Houches Accord 2 [31]. In addition to the actual value of the relic density, micrOMEGAs also provides the contributions of all individual channels contributing to σ ann given in Table II for the two chosen reference scenarios. As can be seen, the processes given in Eqs. (2.1)-(2.3) contribute in a significant manner for both scenarios. More precisely, in scenario I, the scalar top-pair annihilation is the second most important process, and together with the processes previously analyzed in Refs. [17][18][19], it makes up more than 45% of the total annihilation cross section. In this scenario, the mostly "right-handed" scalar topt 1 is the NLSP, and the mass difference between the lightest neutralino and the NLSP is about 20 GeV. Moreover, the processt 1t1 → tt is enhanced by the relatively low gluino mass. Scalar bottom quarks are heavy in this scenario, such that the corresponding annihilation channels are negligible. The processt 1t1 → tt amounts to about 30% TABLE I. Reference scenarios within the phenomenological MSSM for our numerical study. Note that only the parameters which are relevant for our analysis are given here. All dimensionful quantities are given in GeV.  [17,18] accounts for about 15%, such that our next-to-leading SUSY-QCD corrections affect almost 50% of the total annihilation cross section. The importance of different relevant contributions to the total annihilation cross section in and around this scenario is shown in the first four plots of Fig. 1. The part of the parameter space where the lightest neutralino is not the LSP and hence also not the dark matter candidate is indicated in grey. Different shades of green indicate the relative importance of several NLSP annihilation and LSP-NLSP co-annihilation processes. We see that, in general, the coannihilations are most important when the mass splitting between the LSP and the NLSP (in our case the lightest neutralino and the lightest scalar top) is larger, i.e., about 150 GeV, and the dominant contribution shifts to the NLSP annihilations as the mass splitting gets smaller. The parameter region where the dark matter relic density is within 2σ of the experimental value given in Eq. (1.1) is highlighted in all plots in orange. The neutralino relic density is computed using micrOMEGAs. For scenarios which contain lighter neutralinos and lighter stops than the reference scenario I, the total annihilation cross section is dominated by the same processes as in the reference scenario I, and the region where the relic density agrees with the experimental measurement follows an almost straight line where the mass splitting is constant. As we move along the region with the correct relic density towards scenarios with heavier LSP and NLSP, we reach a point where the LSP and NLSP are similar in mass to the light gluino (m˜g ¼ 1495.5 GeV). In these scenarios, we have three particles with almost degenerate masses, and gluino annihilations and co-annihilations with the stop dominate the total annihilation cross section.
The situation is different for scenario II. Here, the lefthanded mass parameter Mq L is much smaller than the righthanded masses Mt R and Mb R , such that the relevant physical statest 1 andb 1 are mainly left-handed with almost degenerate masses. The mass difference between them and the lightest neutralino is about 30 GeV. As a consequence, processes containing botht 1 andb 1 contribute to the annihilation cross section σ ann , as can be seen in Table II. The three processes of our interest contribute to more than 50% of the total annihilation cross section. As we shall discuss later in Sec. II B, the mixed annihilatioñ t 1b1 → tb dominates as compared to stop-pair or sbottompair annihilation. In the last four plots of Fig. 1, we show the relative importance of the channels of our interest in the vicinity of scenario II. Again, the viable region of parameter space where the relic density is within 2σ of the experimental value determined by the Planck satellite closely follows the border between the neutralino and stop LSP regions. In most scenarios along this border, the mixture of the contributing processes is similar to the one presented in Table II. However, around Mq L ∼ 1600 GeV, where the scalar top mass reaches about half of the heavy Higgs mass (m H 0 ¼ 3451.2 GeV), the composition of the contributing processes changes. The stop-antistop annihilation processes enhanced by the Higgs exchange grow in importance. In contrast to the situation around scenario I, for large masses of the neutralino dark matter and the scalar top NLSP, the annihilation and coannihilation processes are not efficient enough to produce the required observed relic density, which can be partly compensated by lowering the mass difference. For even larger masses the region where the relic density is compatible with the Planck measurement features a stop LSP, such that neutralino dark matter would be excluded for M 1 ≳ 1800 GeV.

B. Leading order
Having shown that the processes in Eqs. (2.1)-(2.3) are important in large regions around the two scenarios introduced in the previous section, we now review important features of the leading-order cross sections of these processes. The Feynman diagrams for the processes in question are shown in Fig. 2. The matrix elements of all three processes considered here have contributions from tchannel or u-channel exchanges of strongly interacting gluinos as well as from electroweak gauginos. Therefore, the cross sections can be symbolically written as where σ s is the cross section proportional to the square of the strong coupling constant α 2 s , σ se is the cross section DM@NLO current analysis 30.5% 50.2% DM@NLO total [12][13][14][15][17][18][19] 45.6% 50.2% FIG. 1. Contribution of selected processes to the total annihilation cross section σ ann in the M 1 − Mt R or M 1 − Mq L plane around reference scenarios I or II, respectively. The orange band indicates the parameter region in agreement with the Planck limit given in Eq. (1.1) at the 2σ confidence level. The green levels indicate the relative importance of the processes that can be corrected by DM@NLO (first and fifth plots) and of selected individual processes (remaining plots). The grey region corresponds to mt 1 < m˜χ0 1 . The red dots indicate scenarios I and II of Table I. originating from the interference of the strong and electroweakly interacting parts of the scattering amplitude, and σ e is the purely electroweak cross section proportional to the square of the electromagnetic coupling constant α 2 e . The decomposition of the total cross section into contributions from different channels and interferences of the three processes under consideration here is shown in Fig. 3.
The cross sections fort 1t1 → tt andb 1b1 → bb in the top left, top right and bottom right panels in Fig. 3 show the expected hierarchy, in which the gluino t-channel and u-channel exchanges dominate the cross section and are about an order of magnitude larger than the next largest contribution, which is the interference of the gluino exchange with the electroweak t-and u-channels.  Table I. In the legend, σ Tree denotes the total tree-level cross section, the subscripts g, χ and gχ correspond to gluino exchange squared, gaugino exchange squared, and gluino-gaugino interference. The superscripts indicate the squared t-channel (T), the squared u-channel (U), the sum of both (T þ U), and the t − u interference contributions (TU). For the gaugino exchange, the superscript "Int" refers to the sum of all involved diagrams. Lower part: Contributions relative to the total tree-level result σ Tree .
The contribution from the interference between the gluino exchange diagrams and the gaugino exchange diagrams is yet another order of magnitude larger than the purely electroweak contribution. As argued before, in scenarios where the processes in Eqs. (2.1)-(2.3) are important, the lightest neutralino is bino-like and the gluino mass is relatively small. These facts imply that the neutralinosquark-quark coupling and the gluino-squark-quark coupling differ mainly by the coupling constant. Therefore, the hierarchy observed in Fig. 3 is simply due to the ratio of the different coupling constants α s ð ffiffiffiffiffiffiffiffiffiffiffiffiffi mt 1 mt 2 p Þ=α e ðm Z Þ. The only process where this hierarchy is not present is the annihilationt 1b1 → tb. Here the hierarchy observed in the other two processes is modified due to a few factors. First, there is no gluino u-channel exchange. Then, the gluino mass in this scenario is larger, and this process proceeds also through a chargino u-channel exchange. The larger gluino mass, together with the missing u-channel, suppresses the gluino contribution compared to the other two processes. Moreover, in the case of the Higgsinolike chargino exchange, the Yukawa component of the chargino-squark-quark coupling is not suppressed as in the case of the bino-like neutralino. The combination of these effects results in the interference between the gluino and the chargino exchange being suppressed with respect to the pure gluino contribution only by a factor of about 2. In addition, the electroweak contribution is comparable to the gluino-chargino interference.

C. Color decomposition
Another important aspect of the processes we investigate is the fact that both initial and final state particles carry color. The color structure of the initial and final state will be extremely relevant later in the discussion of the next-toleading SUSY-QCD corrections and their resummation. Both scalar quarks in the initial state (and also the quarks in the final state) transform under the fundamental representation of the SUð3Þ group (denoted here as 3 due to the dimensionality of the representation). The two-particle system, however, transforms under a tensor product of the corresponding representations 3 ⊗ 3 which can be decomposed via a Clebsch-Gordan decomposition into SUð3Þ-invariant subspaces as ð2:5Þ In order to construct a color basis adapted to our matrix element, we can use the Clebsch-Gordan coefficients of the decomposition where the indices a 1;2 can take the values 1 to 3 (for details see Ref. [32]). The basis is constructed by considering that SUð3Þ color symmetry is an exact symmetry of the theory, so the color is conserved between the initial and final states. That means if a pair of initial state particles transforms in an irreducible representation of the SUð3Þ group, the pair of final state particles must transform in the same representation. After proper normalization, we can combine the Clebsch-Gordan coefficients into the following basis relevant for our processes: The matrix element can be expanded in this basis as where s, t, i and j are the color indices of the incoming and the outgoing particles. Given the orthonormality of the basis, the triplet and sextet parts of the amplitude can be determined as In the case of the annihilation processt 1t1 → tt or b 1b1 → bb, the triplet and the sextet matrix elements are a linear combination of the gluino and gaugino t-channel and u-channel exchanges. At tree level the explicit expression for the triplet part of the matrix element is Analogously, the sextet part of the matrix element is The same decomposition can be performed for the process t 1b1 → tb, and the explicit results given in Eqs. ð2:14Þ where due to the orthonormality of the color basis, there is no interference between the triplet and the sextet matrix elements.
The leading-order triplet and sextet cross sections for the relevant processes are shown in Fig. 4. The general behavior of the color decomposed cross sections for the processest 1t1 → tt andb 1b1 → bb is very similar. Both processes contain identical particles in the initial state and are symmetric with respect to their interchange. Given that the color basis vector C f3;3g is antisymmetric with respect to the same interchange, the partial wave of the triplet cross section is a p-wave, making its contribution to the relic density subdominant. For these two processes only the sextet color combination contributes. In the case of the last processt 1b1 → tb, the symmetry argument does not apply, and both color combinations contain an s-wave and contribute equally to the relic density.

III. NEXT-TO-LEADING ORDER
In this section we will discuss the details of our analytical calculation of the full SUSY-QCD corrections to squarkpair annihilation into a pair of quarks. We first concentrate on the virtual corrections and treatment of the UV divergencies in the case of squark-pair annihilation. We continue with the discussion of the treatment of IR divergencies. Finally, we address the Sommerfeld enhancement, and its treatment and impact on the full SUSY-QCD correction to the squark-pair annihilation.
Before we discuss specific details of the next-to-leadingorder calculation, we will address the systematics of SUSY-QCD corrections to processes which at leading order have both strong and electroweak contributions [see Eq. (2.4)]. If we consider any radiative corrections (SUSY-QCD or electroweak) to the processes in question, the cross section including the next-to-leading-order corrections can be symbolically written as parts. Both classes of corrections, the SUSY-QCD and the electroweak, are ultraviolet and infrared finite and gauge independent by themselves, making them formally consistent.
The first and leading term in the NLO correction is Δσ NLO s ðα 3 s Þ which receives contributions only from SUSY-QCD corrections. In particular, these are the SUSY-QCD corrections to the gluino exchange diagrams interfered with the gluino tree-level contribution. These corrections are the main result of this analysis.
The following term Δσ NLO se ðα 2 s α e Þ receives contributions from three sources-from the interference of the SUSY-QCD corrected gluino exchange with the electroweak gaugino exchange, from the interference of the SUSY-QCD corrected electroweak gaugino exchange with the gluino diagrams, and lastly from electroweak corrections to the gluino exchange interfered with the gluino tree level. The last contribution is not included in this analysis, and even though it is formally of the same order, due to the small size of the electroweak corrections which are typically a factor 10 smaller than SUSY-QCD ones, this last contribution is the smallest of the three. In this way our analysis also provides the leading corrections in the term Δσ NLO se ðα 2 s α e Þ. The third term Δσ NLO e ðα s α 2 e Þ contains the interference of the SUSY-QCD corrected electroweak gaugino exchange with the leading-order electroweak gaugino diagrams as well as electroweak corrections to both parts of the interference between the gluino and the electroweak gaugino exchange.
The last term Δσ NLO ee ðα 3 e Þ is not considered here as it contains only electroweak corrections to the electroweak parts of the cross section.
The analysis presented here does not consider electroweak corrections as they are, for the most part, subleading and contribute about 1% to 3% correction [21,22]. In some instances, however, the electroweak corrections, and specifically the Yukawa corrections, can become important [33]. Even though we do not calculate electroweak corrections in this analysis, the leading effects of the enhanced Yukawa corrections are taken into account as described in [17]. In particular, these become relevant in the case of chargino exchange in scenario II (neutralino exchanges are not as enhanced due to the lightest neutralino being a pure bino in both scenarios).
As the discussion below shows, the SUSY-QCD corrections presented here are the dominant corrections even in scenarios with large tan β and are even more dominant owing to the presence of the Sommerfeld enhancement.

A. Virtual corrections and renormalization
The class of processes considered here-the squark-pair annihilation to a pair of quarks-includes strongly interacting particles in the initial state, in the final state, and even in the intermediate state (in the case of the gluino t-and uchannel exchanges). As a consequence, the next-to-leadingorder SUSY-QCD corrections include contributions from vertex corrections, propagator corrections, and box corrections. The corresponding diagrams are displayed in Figs. 5-7, respectively. The next-to-leading-order corrections to the squark-pair annihilation contain one-loop diagrams which are ultraviolet (UV) and infrared (IR) divergent. The UV divergencies are canceled by renormalization of the parameters of the theory and the fields. In order to cancel the IR divergences, one has to properly define an infrared-safe cross section, which is done by also including 2 → 3 processes with an additional gluon being radiated (see Fig. 8).
All one-loop diagrams have been calculated in the SUSY-invariant dimensional reduction scheme (DR) [34,35] where, similar to the minimal subtraction scheme (MS), the number of space-time dimensions is set to D ¼ 4 − 2ε in order to regularize otherwise divergent loop integrals. We have used the standard Passarino-Veltman reduction [36,37] in order to reduce the tensor loop integrals in the one-loop amplitudes to only a few scalar integrals which have then been evaluated using well-known results, e.g., Refs. [38,39]. Our analytical calculations were performed and verified with the help of publicly available tools FeynArts [40], FeynCalc [41], and FORM [42].
In order to cancel UV divergencies of the one-loop amplitude, we renormalize the MSSM by introducing counterterms to the relevant parameters and fields. Because we consider SUSY-QCD corrections to a process involving scalar quarks, quarks and intermediate gluinos, the relevant parameters are the ones that receive corrections proportional to the strong coupling constant α s . Every renormalization scheme is characterized by a careful selection and definition of its input parameters. In a series of previous analyses [17,18], we have put forward a renormalization scheme which combines the advantages of both on-shell and DR renormalization schemes and consistently treats the renormalization of the quark and squark sectors. In these sectors the input parameters are chosen to be the on-shell masses m t , mt 1 , mb 1 , mb 2 , the mass of the bottom quark m b , and trilinear couplings of the third generation A t and A b . The last three parameters are defined in the DR renormalization scheme. We define our renormalization scale as In the following, we will comment on new aspects of the renormalization relevant for our analysis presented here such as the renormalization of the mass of the gluino and its wave function.

Gluino mass and wave-function renormalization
In all our analyses we adopt a convention where not only should the complete next-to-leading-order corrections to the cross section be rendered UV finite, but also all building blocks such as the n-particle irreducible Green's functions should be UV finite as well. This choice requires us to introduce wave-function renormalization constants not only to the fields that correspond to the initial and final state particles but also to fields that give rise to internal propagators. In our case, the only strongly interacting particle that appears in a propagator in our amplitude is the gluino which has not yet been treated within the DM@NLO analysis.
In the case of the gluino, both wave function and mass have to be renormalized in order for the vertex corrections and propagator corrections to be separately UV finite. To this end, we introduce counterterms to the gluino wave function δZ L;R g and the gluino mass δmg as where Π L;R ðk 2 Þ and Π SL;SR ðk 2 Þ are form factors which receive contributions from the corresponding self-energy diagrams.
Even though the gluino is not an external particle in the processes considered in this analysis, we still require that the residue of the propagator at one-loop order is set to unity. This condition fixes the wave-function renormalization constant using the form factors as Using the gluino wavefunction counterterm renders both the propagator and vertex corrections separately UV finite. Moreover, given that the gluino is not an external particle, renormalization of its wave function is not necessary for UV finiteness of the full next-to-leading-order amplitude, so the full amplitude is independent of the gluino wave-function counterterm. This constitutes another consistency check of our analytical calculation.
The mass counterterm is determined from the on-shell condition which requires that the gluino mass m˜g, which is an input parameter, is identical to the position of the pole of the gluino propagator. It is given as 2. Some remarks on the renormalization of the squark sector Even though we have discussed the details of the renormalization of squark parameters in Ref. [17], we would like to remark here on one feature of the renormalization scheme relevant for evaluating the results of this analysis. As discussed in Ref. [17], we use the relation between the nondiagonal squark mass matrices for up-type and down-type squarks, to relate the input parameters in the whole squark sector, which are defined in different renormalization schemes. In the next step, we determine the dependence of the soft supersymmetry-breaking squark mass parameters M 2Q and M 2 fŨ;Dg of the three on-shell masses mb 1 , mb 2 , mt 1 and the input parameters contained in m 2 LR (A q , m q ). In the sbottom sector the parameters M 2Q and M 2D are contained in the matrix elements m 2 b;LL and m 2 b;RR , which are given by the on-shell masses as ð3:11Þ One notices that there are two possible values for the parameters M 2Q and M 2D and consequently also for the third parameter M 2Ũ which can be found in one of the diagonal elements of the nondiagonal scalar top mass matrix and is related to the first two parameters through The parameter M 2Q is common to both elements m 2 b;LL and m 2 t;LL . Given the freedom to choose from two possible solutions for the squark soft supersymmetry-breaking mass parameter, we can end up with two possibly very different mass matrices and two different sets of mixing matrices (three out of four masses of the squarks would be the same in both cases as they are used as input). In order to ensure a naturally small correction to the mixing matrices when changing between our scheme and the DR renormalization scheme, we always select the solution which preserves the hierarchy between the mass matrix elements in the scalar top quark sector m 2 t;LL and m 2 t;RR which was present in the pure DR scheme.

B. Real corrections
After treating the ultraviolet divergencies and removing them by renormalization, we are left with the IR divergencies in the one-loop amplitude. The IR divergencies are also regularized by performing calculations in a general dimension D ¼ 4 − 2ε and are subsequently removed by considering an IR-safe observable, which in our case is a cross section with one additional gluon in the final state. The corresponding diagrams are depicted in Fig. 8. The cross section for the radiation of an additional gluon cannot be calculated analytically. On the other hand, the cancellation of the IR divergencies between the virtual corrections and the real radiation cross section has to be performed analytically. One of the methods for how to extract the IR divergencies out of the real radiation cross section is the phase-space slicing method [37,43,44]. This method is based on the fact that the infrared divergence is connected to a specific configuration of the momentum of the gluon. The soft infrared divergence arises when the additional gluon's energy vanishes, whereas the collinear infrared divergence arises when the additional gluon is radiated collinearly to the momentum of an external massless particle. The phasespace slicing method uses kinematical cuts to divide the three-particle phase space into regions where either one or both of the aforementioned configurations of the gluon 4-momentum occur and the remainder where the cross section is IR finite. In the singular regions the full matrix element is replaced by an approximation which can be integrated analytically making the IR divergence explicit. In the nonsingular region the full matrix element can be integrated numerically without any obstacles.
In our case all external particles are massive, so the 2 → 3 cross section contains only a soft IR divergence. In the singular region where the gluon's energy is smaller than an arbitrarily small cutoff ΔE, we use the soft gluon approximation which factorizes the differential cross section as dσ dΩ where F contains the integral over the phase space of the gluon, and therefore also the divergence. Here, k is the 4-momentum of the gluon and a and b are 4-momenta of two external particles which can emit a gluon. These integrals are given in Refs. [37,45]. The cutoff ΔE introduced to separate the singular from the nonsingular region in the three-particle phase space enters into the calculation of the soft-gluon radiation as well as into the integration of the 2 → 3 cross section over the nonsingular phase space. In principle, the dependence on this cutoff should disappear in the sum of the contributions from both phase-space regions, but in practice, the independence on the cutoff is limited by the numerical stability of the integration over the nonsingular region. We have investigated the dependence on the cutoff and found that the integration in our case is stable and independent of the cutoff for a relatively large interval of cutoffs around ΔE ¼ 10 −5 ffiffi ffi s p .
The complete result after we have included all virtual corrections, the counterterms and the real radiation 2 → 3 cross section is UV and IR finite. In contrast to the leadingorder result which consists of a cross section with two particles in the final state and is implemented in micrOMEGAs, the complete result is a consistent combination of a one-loop corrected cross section with two particles in the final state and a leading-order cross section with three particles in the final state. In DM@NLO, the complete result replaces the leading-order result of micrOMEGAs.

C. Sommerfeld resummation
When calculating the relic density in our case, an important contribution comes from the annihilation of squarks moving with nonrelativistic velocities. If annihilating, nonrelativistic particles couple to much lighter force mediators which in our case are the gluons, the annihilation cross section is modified due to the wellknown Sommerfeld effect [46]. The reason for this modification is that the exchange of n gluons between the initial state squarks (see Fig. 9) contains a correction proportional to ðα s =v rel Þ n . This correction becomes significant and can spoil the perturbative expansion when the relative velocity of the squark pair v rel is comparable to the strong coupling constant α s . In such a case these contributions have to be resummed to all orders leading to the Sommerfeld effect.
Small relative velocities occur naturally in the freeze-out regime, E kin ∼ T FO ∼ m˜χ0 1 =25, and therefore the Sommerfeld resummation is expected to be relevant in the case of dark matter annihilation in general and in our case in particular. As for our processes of interest,q iq 0 j → qq 0 , the cross section is dominated by the s-wave component (see Fig. 4). We can factorize the resummed cross section as ðσvÞ resum ¼ S 0;½3 ðσvÞ Treẽ qq→qq;½3 þ S 0;½6 ðσvÞ Treẽ qq→qq;½6 ; ð3:15Þ where we have split the leading-order cross section according to its color contribution to the triplet and sextet configurations (see Sec. II C). Here S 0;f½3;½6g indicate the corresponding s-wave Sommerfeld factors, whose evaluation we discuss in the following. In the nonrelativistic limit, the resummation of the gluon exchange diagrams as shown in Fig. 9 amounts to solving the Schrödinger equation with the corresponding Coulomb potential. The Coulomb potential including gluon loops at next-to-next-to-leading order was evaluated in [47] and extended by fermion loops in Ref. [48]. In the MS scheme, the Coulomb potential reads [49] ð3:16Þ with γ E ¼ 0.5772 being the Euler-Mascheroni constant. Furthermore, we have defined where b 0 corresponds to the one-loop β-function coefficient with C A ¼ 3 and T F ¼ 1=2. We treat the top as the only massive quark, such that we set the number of massless quarks to five (n f ¼ 5). The Coulomb potential, given in Eq. (3.16), describes the interaction of any nonrelativistic colored particles transforming in general SUð3Þ-representations R 1 and R 2 . The color structure of such a scattering process can be decomposed as The color factor C ½R is given in terms of the quadratic Casimir operators of the relevant SUð3Þ representations as In the case considered here, the two squarks in the initial state both transform under the fundamental representation of SUð3Þ, and the color decomposition is 3 ⊗ 3 ¼3 ⊕ 6.
Using the quadratic Casimir operators for the fundamental and the sextet representations, we obtain [50] ð3:21Þ where the Green's function G 0 ð0; E þ iΓt 1 Þ stands for the solution of the Schrödinger equation without any Coulomb potential. The solution to Eq. (3.23) at the origin is well known [53], and we consider here all terms up to NLO, where the LO and NLO contributions are given by In Eqs. (3.26) and (3.27) and in the following, we use the short-hand notation Moreover, ψ ðnÞ ¼ ψ ðnÞ ð1 − κÞ is the nth derivative of ψðzÞ ¼ γ E þ d=dz ln ΓðzÞ with the argument (1 − κ), state particles, the parameter v s in Eq. (3.29) is the nonrelativistic velocity of one of the incoming particles and should not be confused with v rel ¼ 2v s , the relativistic, relative velocity of the two annihilating particles. In order to calculate the Sommerfeld factor in Eq. (3.24), we also need the Green's function that solves the system without any potential term in Eq. (3.23), which is given by Finally, we need to fix the scale μ C that appears in the potential and has an impact on the evaluation of α s in the Sommerfeld factor. We follow here the treatment presented in Ref. [50] and set where 4m red v s is motivated by the typical momentum exchange of the gluons in the ladder diagram and the scale μ ½R B corresponds to twice the inverse Bohr radius r B . It is defined via ð3:33Þ In order to obtain μ ½R B , we solve Eq. (3.33) iteratively. One important remark should be made here. In the case of a repulsive QCD potential, the Sommerfeld correction causes a large reduction of the annihilation cross section which for extremely small velocities should vanish. In order to prevent unphysical results such as negative cross sections, we modify the definition of the scale μ ½R C as the smallest scale where the cross section is positive for v ¼ 0. In our specific case a scale defined in this way is only about 2% larger than the scale μ ½R C defined through Eq. (3.32).
As the box diagrams in the full NLO calculation also contain the velocity-enhanced part of the one-gluon exchange, which is at the same time already included in the Sommerfeld resummation, we have to subtract this contribution in order to avoid any double counting.
To isolate the velocity-enhanced term from the box contribution, we expand the box contribution in the relative velocity (for details see the Appendix). We then construct the subtracted cross section ðσvÞ sub NLO based on the expanded matrix element of the box diagrams given in Eq. (A9). The leading velocity-enhanced term of the subtracted cross section is we see that, given C ½R ¼ C ½R box and v s ¼ v rel =2, the two expressions differ in the scale at which the strong coupling constant is being evaluated. While in the perturbative NLO calculation α s is evaluated at the renormalization scale μ R ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi mqmq0 p , in Sommerfeld resummation the characteristic scale μ C is used. By choosing to use Eq. (3.34) to avoid the double counting, we make use of the fact that the natural scale used in the description of the interaction between incoming particles at small velocities is μ C . This is consistently used to all orders in the resummed cross section ðσvÞ resum given by Eq. (3.15). The full next-toleading-order cross section, consistently including also the Sommerfeld resummation, reads The Sommerfeld enhancement of the annihilation cross section can be caused not only by the multiple exchanges of gluons but also by the exchanges of photons or even gauge or Higgs bosons. Coincidentally, in the processes considered in this analysis, the Sommerfeld enhancement from the photon exchange has the same effect on the cross section as the Sommerfeld enhancement from the exchange of gluons. In the case of photons, the Sommerfeld effect leads to a reduction of the cross section for the processes with samesign particles in the initial state such ast 1t1 → tt and b 1b1 → bb. If the electric charges of the particles in the initial state have different sign, e.g.,t 1b1 → tb, the electromagnetic interaction is attractive and the Sommerfeld effect causes an enhancement of the cross section. Given that the photon exchanges are suppressed with respect to the exchanges of the gluons by a factor α EW ðμ C Þ=α s ðμ C Þ, we neglect their effect in this analysis.
Because of the large trilinear couplings in both scenarios, it might be interesting to study the Sommerfeld enhancement coming from the exchange of Higgs bosons [54][55][56][57][58][59]. In the regime of Sommerfeld enhancement, bound state formation can also potentially occur, giving rise to new annihilation channels and thus altering the relic density prediction. This has been previously studied for stopantistop annihilation for both gluon [60] and Higgs exchange [61]. Such studies, however, are far beyond the scope of this work and are left for future analyses.

D. The NLO cross section results
In this section we present the first result of our analysis, which is the impact of SUSY-QCD next-to-leading-order corrections on the annihilation cross sections of scalar top or bottom pairs. Apart from the cross section (or more precisely σv) we also show in arbitrary units the Boltzmann distribution function which is involved in the calculation of the thermal average hσvi at the freeze-out temperature (grey shaded area). It should serve as a reminder that the cross section contributes to the determination of the relic density only in a limited range in the center-of-mass momentum p cm .

Scenario I
In the first scenario introduced in Sec. II, the mass splitting between the lightest neutralino and the lightest stop quark is relatively large. As a result the dark matter annihilation cross section receives important contributions not only from the stop-pair annihilation into top quarks but also from the neutralino-stop co-annihilation into a top quark and a gluon and other final states. The results of SUSY-QCD corrections for these processes are shown in Fig. 10. As described in the previous section, the next-toleading-order cross section consists of vertex corrections, propagator corrections, box corrections, counterterm contributions and the real radiation cross section, which has to be added to render the prediction infrared finite. The corrections from each contribution are not shown in Fig. 10 due to the cancellations of the ultraviolet and infrared divergencies between the contributions, which make each contribution on its own ill defined. In addition to the NLO cross section, we also include the enhanced higher-order contributions stemming from the nonrelativistic Coulomb correction.
In the case of the annihilation of a pair of scalar top quarks, both initial particles are colored, and in the limit of vanishing relative velocity of the squark pair, the Coulomb corrections dominate the full corrected cross section. The origin of these corrections is the exchange of multiple gluons between a pair of slowly moving squarks in the initial state, and the details were discussed in Sec. III C. The effect of the Coulomb corrections strongly depends on the color multiplet, in which the pair of squarks transform. Based on the color decomposition presented in Sec. II C, the annihilation cross sectiont 1t1 → tt in scenario I is dominated by the contribution where the squark pair forms an SUð3Þ-sextet (see Fig. 4). In this representation the multiple exchange of the gluons can be described by a repulsive nonrelativistic QCD potential, as discussed in detail in Sec. III C. That is why the Coulomb corrections in this case cause a reduction of the cross section. As already discussed in the previous section, the next-to-leading-order cross section contains the one-loop contribution also included in the Coulomb enhancement. This one-loop contribution can be traced back to all box diagrams in Fig. 7, where one gluon is exchanged between the incoming squarks. The contribution from this class of diagrams dominates the one-loop cross section for small velocities and is so large that it causes the cross section with one-loop corrections to be negative (see the green dashed line in Fig. 10). As discussed in Sec. III C, in order to prevent double counting, we remove the part of the box contribution which is already included in the Coulomb resummation. This allows us to quantify the pure one-loop correction to the annihilation cross section without any enhancement (red dash-dotted line in Fig. 10). We see that the one-loop correction without the enhancement is a large positive correction of about 30%-40% over a large range of p cm .
Comparing the result for the nonenhanced NLO cross section with the full result, which is the sum of the nonenhanced NLO cross section and the Sommerfeld corrections, shows that the latter are important for all relevant values of p cm . Starting at the largest value of p cm ∼ 600 GeV which is still relevant for the determination of the relic density, we observe that the Coulomb corrections already reduce the constant 30% NLO correction by a few percent. For smaller relative velocities corresponding to p cm ∼ 150 GeV, the NLO correction is fully canceled by the Coulomb corrections, and for very slow velocities the Coulomb corrections take over and the overall correction is large and negative. For almost vanishing velocities the total cross section vanishes completely. This is due to the fact that the dynamics of the squark pair in the regime when Coulomb corrections are very large (meaning for vanishing velocities) correspond to a motion of the pair in a highly repulsive QCD potential. This in turn means that large repulsive forces repel one squark from the other, reducing the probability of annihilation and thereby reducing the cross section.
In summary, we can conclude that SUSY-QCD corrections tot 1t1 → tt are sizable either through the one-loop corrections for large p cm or the enhanced Coulomb corrections for small p cm .
The co-annihilation processes important in this scenario were discussed in detail in Refs. [17,18]. In Fig. 10 we also show the effect of the SUSY-QCD corrections on the coannihilation cross sections in scenario I. We see that the next-to-leading-order corrections in the case of co-annihilations are substantial, ranging from −30% in the case of the co-annihilation into the top quark and Higgs-boson final state to þ50% in the case of the top gluon final state.
There are a few substantial differences such as the fact that the corrections are negative in the case of co-annihilations with electroweak bosons or Higgs bosons in the final state or that there is a large difference between our leading-order result and the micrOMEGAs result for the co-annihilations with electroweak and Higgs bosons, which can be traced to a different definition of underlying parameters in our renormalization scheme. The next-to-leading-order correction with respect to the micrOMEGAs result is largely reduced to at most −10%. Given that our leading-order prediction for the co-annihilations into the top quark and a gluon coincides with the micrOMEGAs prediction, the large next-to-leading-order correction also gives directly the correction with respect to the micrOMEGAs result.

Scenario II
In the second scenario, the choice of parameters such as tan β and the gaugino and squark mass parameters causes the masses of the lightest neutralino, and the lightest scalar top and bottom quarks to be almost degenerate. This leads to different processes contributing significantly to the total dark matter annihilation cross section. The smaller mass difference renders co-annihilations ineffective, and the fact that three particles are mass degenerate leads to a larger number of annihilations. Moreover, the large value of tan β enhances the gaugino exchange in the case of the stopsbottom annihilation, and this, together with a different color structure, makes this annihilation dominant in the case of scenario II.
The full next-to-leading-order results for three dominant processes are shown in Fig. 11. The processest 1t1 → tt and b 1b1 → bb have very similar features to the annihilation of a pair of scalar top quarks in scenario I. The main difference in this scenario is the processt 1b1 → tb, which has an entirely different decomposition of the leading-order cross section in terms of the t-and u-channel exchanges combined with a different color decomposition, which is essential in explaining the behavior of the NLO cross section. Similar to the already discussed case oft 1t1 → tt, the NLO correction contains a velocity-enhanced term, which is already resummed in the Sommerfeld correction. In order to avoid double counting, we define again the nonenhanced NLO correction σ NLO B where we subtract the term which is already accounted for by the Sommerfeld resummation (red dash-dotted curve in Fig. 11). As one can see in Fig. 11, the nonenhanced NLO correction is substantial in all processes in scenario II. In the case of stop-pair or sbottom-pair annihilations, this NLO correction is compensated by a large and negative Sommerfeld correction which is derived here from a repulsive QCD potential. The color decomposition of thet 1b1 → tb shows (see Fig. 4) that in contrast to the other processes, the cross section here is dominated by the part where the initial stop and sbottom quarks transform as an SUð3Þ triplet. In this color configuration a pair of slowly moving squarks experiences an attractive strong force, which leads to a large enhancement of the annihilation cross section. Comparing the full result (solid blue line in Fig. 11) with the result containing just the nonenhanced NLO corrections, we see that the Sommerfeld enhancement is important over the whole region in p cm that is relevant for the calculation of the relic density.
It is worth mentioning that due to the large value of tan β in scenario II and to the fact that the chargino u-channel exchange gives an important contribution to the cross sectioñ t 1b1 → tb, the Yukawa corrections to the chargino exchange can give a non-negligible contribution. We have included the tan β dependent Yukawa corrections even beyond next-toleading order in the full result, and we show their effect separately in Fig. 11 (yellow dashed line). Even though the Yukawa corrections are non-negligible, they are small (about 3%) compared to the remaining SUSY-QCD corrections or to the Sommerfeld enhancement.
The full correction tot 1b1 → tb is larger than 50% over the whole relevant range of p cm and can even exceed 100% (without threatening perturbativity as this correction originates from a resummation). We will show in the next section that using the fully corrected annihilation cross sections in scenario II has a large impact on the relic density.

IV. IMPACT ON THE NEUTRALINO RELIC DENSITY
We finally come to the discussion of the impact of the corrections presented in Sec. III on the neutralino relic density Ω˜χ0 1 h 2 . To this end, we have implemented the corrections into the numerical code DM@NLO, which is used as an extension to micrOMEGAs. In practice, this means that the Boltzmann equation is still integrated using the latter, while the cross section calculation of the relevant processes (see Table II) is performed by DM@NLO instead of CalcHEP.
In the following, we will illustrate the impact of the corrections by comparing the relic density obtained using the full DM@NLO NLO calculation to the values obtained using the tree-level calculation of micrOMEGAs/CalcHEP.

A. Scenario I
We start by examining the impact of NLO corrections in the vicinity of scenario I, where we compare the relic density obtained from the micrOMEGAs calculation to the one obtained using our full NLO result as presented in Sec. III. The impact is illustrated in Fig. 12, where we show the corresponding viable regions of parameter space in the M 1 -Mt R plane. As can be seen, the favored parameter region where the calculated relic density satisfies the experimental constraint, Eq. (1.1), is shifted towards smaller mass parameters in order to compensate the increased annihilation cross section. It is important to note that this shift is larger than the width of the band which corresponds to the Planck 2σ uncertainties.
The situation changes for higher masses, where the processes discussed in this work and corrected at the NLO level are less relevant. The corrections of the remaining processes relevant in this part of parameter space are left for future work.
In Fig. 13, we show the same results in the vicinity of scenario I, but this time projected onto the plane of the physical neutralino and stop masses. Note that, here, the variation of the physical masses solely stems from varying the parameters M 1 and Mt R , respectively, while all other soft parameters, including those that, in general, may influence the neutralino and stop masses, are kept fixed to the values of Table I. From Fig. 13, we see that the cosmologically favored region of parameter space is shifted by about 7 GeV in both the neutralino and the lighter stop mass. These results lead to the conclusion that the corrections presented in this work are relevant when performing an extraction of either physical masses or fundamental model parameters from cosmological data. Let us note that this conclusion is the same as for the previous analyses of other processes entering the calculation of Ω˜χ0 1 h 2 [12][13][14][15][17][18][19]. In order to get a better understanding of the impact of the different contributions to the annihilation cross section σ ann , and consequently to the relic density Ω˜χ0 1 h 2 , we have performed a one-dimensional scan along the region where micrOMEGAs predicts the correct relic density, varying simultaneously the parameters M 1 and Mt R . The result of this scan is shown in Fig. 14 as a function of both parameters, while all other parameters are fixed to the values given in Table I. First, it can be noticed that our tree-level prediction differs from the micrOMEGAs result. This is a direct consequence of the corresponding difference in the annihilation cross sections, as discussed in Sec. II B and shown, e.g., in Fig. 10. Taking into account the corrections discussed in Sec. III, it can be seen that the total correction is split into two parts associated with the relevant classes of processes, namelyt 1t1 → tt andχ 0 1t1 → qg, qV, qϕ. Figure 10 shows that the correction to the relic density in scenario I is dominated by the corrections to the co-annihilation processes even though at leading order these processes contribute about a factor 2 less thant 1t1 → tt. This is a consequence of the Sommerfeld suppression of the annihilation cross section as discussed in Sec. III D. Moreover, we see that for lower bino/squark mass parameters M 1 =Mt R , the correction to the co-annihilation processes is numerically more important than for large mass parameters. This is explained by the fact that the relative importance of the co-annihilation processes is higher in this region of parameter space (see, e.g., Fig. 1). Moving towards higher values of M 1 , the relative importance of the stop-pair annihilation increases and, consequently, the associated correction to the relic density becomes more important.
Overall, the relic density obtained using our full (i.e., NLO including resummation) calculation is about 6% smaller than the one obtained by using micrOMEGAs. Again, we emphasize that this shift is more important than the uncertainty given by the Planck measurement, which is at the 2σ confidence level of about 2%.

B. Scenario II
Let us now focus on scenario II, where not only is t 1t1 → tt relevant, but the related processesb 1b1 → bb andt 1b1 → tb also give sizable contributions to the annihilation cross section σ ann . Therefore, they also need to be corrected at the NLO level including the resummation, as discussed in Sec. III. Again, we start by depicting the parameter region compatible with the measured values for the relic density (see the second row of plots in Fig. 12) in the vicinity of scenario II-in this case, in the M 1 − Mq L plane-which are the relevant neutralino and squark mass parameters. Similar to scenario I, as discussed above, the viable regions with respect to the relic density are shifted towards lower masses. Here, the shift is again more important than the uncertainty and is much larger than in scenario I. It corresponds to a shift of about 17 GeV in the bino mass parameter M 1 and about 15 GeV in the left-handed squark mass parameter Mq L . In terms of physical masses, shown in Fig. 13, this corresponds to a shift of about 17 GeV for the neutralino mass and of about 15 GeV in the lighter stop mass. Once more, these findings underline the importance of the presented corrections in the light of precision cosmology. It is to be noted that for this part of the analysis, only M 1 and Mq L have been varied, and the results have been projected on the so-obtained plane of the physical neutralino and stop masses, while all other parameters have remained fixed to the values given in Table I.
In order to decompose our full NLO prediction for the relic density into the contributions from different processes, we show in Fig. 14 the NLO corrected contributions to the relic abundance from individual processes along the region where micrOMEGAs predicts the correct relic density, simultaneously varying the parameters M 1 and Mq L . We see that using our tree-level annihilation FIG. 14. Upper part: Neutralino relic density Ω˜χ0 1 h 2 along the parameter region satisfying the experimental constraints on the relic abundance. Both the bino mass parameter M 1 and the squark mass parameter Mt R (or Mq L ) are shown around scenario I (or scenario II). In both plots we show the values obtained using micrOMEGAs (Ω MO h 2 ), our tree-level calculation (Ω Tree h 2 ), and our full one-loop calculation including the resummation (Ω F h 2 ). For scenario I, we also show the value obtained by correcting only neutralino-stop coannihilation (Ω˜χ0 1t 1 h 2 ), and the value obtained by correcting only stop-pair annihilation (Ω˜t 1t1 h 2 ). Similarly, for scenario II, we show the effect of correcting only the stop and sbottom-pair annihilations (Ω˜q 1q1 h 2 ) and the stop-sbottom annihilations (Ω˜t 1b1 h 2 ). Lower part: Impact of the different contributions relative to the relic density obtained by using micrOMEGAs. cross section with differently defined input parameters shifts the relic density by a few percent (black dashed line). In all remaining contributions we use our tree-level annihilation cross section for all relevant processes. Starting from our tree level, a very small correction of about 1% in the whole region comes from including NLO corrections only to thet 1t1 → tt andb 1b1 → bb. It is important to point out that the reason for this extremely small correction is the Sommerfeld enhancement, which in this case in fact suppresses the cross section. This is due to the repulsive nature of the dominant SUð3Þ-sextet contribution. Even if the full correction to the cross section is large and negative for small p cm (see Fig. 11), interestingly the correction to the thermal averaged cross section is still positive, leading to a drop in relic density. The largest contribution comes from the annihilation cross section oft 1b1 → tb. The first reason is the large contribution of this process to the total annihilation cross section already at tree level (see Table II). The second reason is the large Sommerfeld enhancement emerging from the attractive potential of the dominant SUð3Þ-triplet contribution, which in this case makes the full correction to the cross section extremely large. We see that depending on the dominant contribution of the color decomposition, SUð3Þ-sextet or SUð3Þ-triplet, the Sommerfeld correction either suppresses or enhances the cross section such that the total SUSY-QCD correction to the relic density over the whole range is about 25%. This results in the visible shift of the preferred parameter region, which is much larger than the experimental uncertainty given in Eq. (1.1).

V. CONCLUSION
Scenarios in the MSSM with light stops are still very appealing due to their potential to address many problems that the MSSM with heavy particles might have. In such scenarios, squark-pair annihilations into quarks are often very important processes that govern the annihilation of dark matter, which is typically the lightest neutralino.
We have analyzed two such example scenarios, which pass all current experimental constraints. We focused on the SUSY-QCD corrections to squark-pair annihilations into quarks. We have reviewed the details of the one-loop calculation and of the Sommerfeld enhancement. We have shown that the one-loop corrections of the cross sections are sizable even without the Sommerfeld enhancement. The Sommerfeld corrections are shown to cause two different effects depending on the nature of the strong force between the pair of incoming scalar quarks. In the case of annihilations between the same type of squarks, the enhancement turns into a reduction of the cross section as here the squarks experience a strong repulsive force. If the scalar quarks are different, however, the cross section is strongly enhanced due to the attractive strong force. Finally, we have investigated the impact of these corrections on the predicted relic density. We have demonstrated in our typical scenarios that even with the Sommerfeld reduction, the corrections are larger than the experimental uncertainty. In case of an enhancement, the corrections cause a 25% shift in the preferred parameter region where the relic density satisfies the experimental constraints.

ACKNOWLEDGMENTS
The authors would like to thank A. Pukhov for providing us with the possibility to interface our DM@NLO code with micrOMEGAs. J.

APPENDIX: SMALL VELOCITY EXPANSION OF THE BOX CONTRIBUTION
In order to subtract the velocity-enhanced part of the NLO contribution that is already included in the Sommerfeld resummation, we expand the corresponding contribution of the box diagrams in the relativistic relative velocity [63] v rel ¼ where the tensor integral has the arguments as in Eq. (A2) with M 2 ¼ m χ and the color factors C ½R box are given as The scalar integral for the specific arguments from Eq. (A2) can be written as Li 2 ðx 13 ; x k 12 ; x l 23 Þ − π 2 6 : ðA5Þ The generalized polylogarithm in Eq. (A5) is defined as [38] Li 2 ðx 1 ; …; The variables x ij are defined using the loop masses M i and M j as well as the invariant combinations of 4-momenta p 2 ij , as given in Eq. (A2), as In our case the only velocity dependent x ij is x 13 , which for M 1 ¼ m 1 , M 2 ¼ m χ and M 3 ¼ m 2 gives having used the relative velocity as defined in Eq. (A1). Given that the color factors C ½R box are independent of the exchanged gaugino, we can simplify the enhanced matrix element as