SUSY-QCD corrected and Sommerfeld enhanced stau annihilation into heavy quarks with scheme and scale uncertainties

We investigate stau-antistau annihilation into heavy quarks in the phenomenological Minimal Supersymmetric Standard Model within the DM@NLO project. We present the calculation of the corresponding cross section including corrections up to $\mathcal{O}(\alpha_s)$ and QED Sommerfeld enhancement. The numerical impact of these corrections is discussed for the cross section and the dark matter relic density, where we focus on top-quark final states and consider either neutralino or gravitino dark matter. Similarly to previous work, we find that the presented corrections should be included when calculating the relic density or extracting parameters from cosmological observations. Considering scheme and scale variations, we estimate the theoretical uncertainty that affects the prediction of the annihilation cross section and thus the prediction of the relic density.


I. INTRODUCTION
More than 80 years after its first observation [1], the existence of dark matter in our Universe is now well established by various coinciding observations (see review [2] and references therein). In the absence of a clear consensus about the exact nature of dark matter, numerous theoretical models have been developed to explain its nature. Most of such models are based on the hypothesis that dark matter consists of weakly interacting massive particles (WIMPs) which achieve the observed relic abundance through thermal freeze-out [3,4].
In the present work, we focus on the Minimal Supersymmetric Standard Model (MSSM), which provides suitable WIMP candidates [5] such as the lightest neutralino or the gravitino. Assuming that R-parity is conserved, these particles are stable, and they interact only weakly as required by the WIMP paradigm.
Over the last decades, the relic abundance of cold dark matter (CDM) within the cosmological ΛCDM model has been determined to a very good precision, where h denotes the present Hubble expansion rate in units of 100 km s −1 Mpc −1 . This interval has been obtained from the cosmic microwave background measurements by the Planck satellite [6] combined with polarization data from the WMAP mission [7]. Within a given particle physics model, such as the MSSM, the relic density of dark matter can be theoretically predicted, allowing to identify cosmologically favoured regions of parameter space. More precisely, for a dark matter candidate χ with mass m χ , the predicted relic density Ω χ h 2 is obtained through Here, ρ crit stands for the critical energy density of the Universe and n χ for today's number density of the dark matter candidate. The value of n χ corresponds to the solution of the Boltzmann equation [8][9][10] dn χ dt = −3Hn χ − σ ann v n 2 χ − n eq χ 2 . (1.3) This differential equation describing the time evolution of the number density contains the thermally averaged cross section σ ann v of the annihilating neutralinos. Dark matter annihilation cross sections typically lead to a relic density exceeding the limits given in Eq. (1.1) causing the overclosure of the Universe. As a consequence, the annihilation cross section needs to be enhanced by some mechanism, such as resonant annihilation or efficient co-annihilations, resulting in lower values of the relic density. In the present paper, we will focus on the latter case, in particular we will assume that the lightest neutralino and the lightest stau are almost degenerate in mass.
Accounting for co-annihilations, the thermally averaged annihilation cross section appearing in Eq. where the relative velocities are given by v ij = (p i · p j ) 2 − m 2 i m 2 j /(E i E j ). The ratios of equilibrium densities appearing in Eq. (1.4) are suppressed via the so-called Boltzmann factors n eq i n eq χ ∝ exp − m i − m χ T . (1.5) This implies that only co-annihilation of dark matter with particles that are almost degenerate in mass can have a sizeable impact on the relic density. In the present paper, we assume that the lightest neutralino and the lighter stau are very close in mass. The very narrow observational interval given in Eq. (1.1) clearly calls for very precise theoretical predictions. Within public codes calculating the dark matter relic density for new physics models, such as micrOMEGAs [11][12][13][14][15] or DarkSUSY [16][17][18], the underlying processes entering Eq. (1.4) are evaluated usually at tree-level, including effective couplings capturing certain higher order effects, e.g., for the Yukawa couplings. It is the goal of the DM@NLO project to provide a more accurate calculation of the annihilation cross section and thus of the dark matter relic density.
The present paper will add to the above list of processes by presenting the corrections of order α s to stauantistau annihilation. This process may be relevant in scenarios with neutralino or gravitino dark matter. In the following, we start by discussing the phenomenological impact of stau-antistau annihilation in Sec. II. In Sec. III, we will first discuss the technical details of the NLO calculation. Moreover, we will present the QED Sommerfeld enhancement included in our calculation. In Sec. IV, we will illustrate the effect of the corrections in typical scenarios within the phenomenological MSSM for both neutralino and gravitino dark matter. We will also discuss the theoretical uncertainty coming from the variations of the renormalization scale and the renormalization scheme. Conclusions are given in Sec. V.

II. PHENOMENOLOGY RELATED TO STAU-ANTISTAU ANNIHILATION
Let us start by discussing the phenomenology of stauantistau annihilation in the context of the dark matter relic density. This process may become relevant in two cases.
First, as mentioned in the introduction, if the stau is very close in mass to the neutralino, neutralino-stau coannihilation as well as stau pair-annihilation will contribute in a sizeable manner to the total annihilation cross section σ ann . In this case, the prediction of the neutralino relic density is obtained directly from solving the Boltzmann equation as explained in the introduction.
The second situation is the case where the dark matter candidate is the gravitino, denoted asG, which is the spin-3/2 superpartner of the graviton, if local supergravity is assumed. In this situation, the next-to-lightest supersymmetric particle may be either the lightest gaugino or the lightest sfermion, for example the lighter stau. The phenomenology related to the additional gravitino is governed by a single additional parameter, which is the gravitino mass mG. Recent detailed discussions of gravitino dark matter within the phenomenological Minimal Supersymmetric Standard Model (pMSSM) can be found, among others, in Refs. [39][40][41].
In this situation, the gravitino relic density receives contributions from thermal and non-thermal production, The thermal contribution depends on the reheating temperature, T R , and the gluino mass, mg, and can be approximated as [42] Ω th G h 2 0.27 The corresponding full expression has been derived in Refs. [43,44]. In the present work we focus on cases where the non-thermal contribution dominates, such that it is reasonable to rely on the simplified expression given in Eq. (2.2).
The non-thermal contribution arises from the decay of the lighter stau. If R-parity is conserved, each stau can decay only into a gravitino, and the corresponding contribution to the gravitino relic density is obtained by reweighing the would-be relic density of the stau, obtained from integrating the Boltzmann equation, according to [45] Ω non−th In this context, the stau lifetime is constrained in order to preserve the abundances of light elements in the early Universe, which are well explained by primordial nucleosynthesis. In terms of stau and gravitino masses, this constraint can be approximated as [46,47] tτ 1 6100 s 1 TeV mτ 1 5 mG 100 GeV 2 6000 s , (2.4) implying that the gravitino mass is about one order of magnitude smaller than the stau mass.
In the following, we illustrate the phenomenology, and later on also the impact of higher-order corrections to the stau annihilation cross section, for two example scenarios within the pMSSM, where the 19 independent softbreaking parameters are defined at the supersymmetry (SUSY) scale Q = √ mt 1 mt 2 . In order to identify representative parameter configurations, we make use of the pMSSM analysis presented by ATLAS in Ref. [48]. The parameter points resulting from this study satisfy the imposed constraints from LHC searches, the observed Higgs mass of m h 0 ∼ 125 GeV [49][50][51], and rare decays such as, e.g., b → sγ and B s → µ + µ − . Let us also stress that the corresponding parameter region is robust against more recent searches performed by the ATLAS and CMS collaborations [52].
Inspired from the parameter points given in Ref. [48], we have chosen a typical scenario for our study related to neutralino dark matter. All pMSSM parameters related to this scenario, labelled I, are given in Tab. I, together with the corresponding neutralino and stau masses, which are relevant in our study.
As the ATLAS analysis only covers the case of neutralino dark matter, we have defined a second scenario, labelled II, for our study of gravitino dark matter. This second scenario aims at illustrating the situation in a rather simple way, the soft parameters being chosen such that only the staus are relatively light with masses just below 2 TeV, while all other states are rather heavy with masses of about 5 TeV. The corresponding soft-breaking parameters, together with the relevant masses, are displayed in Tab. I. In addition, in order to satisfy the lifetime constraint mentioned above, the gravitino mass will be chosen to be around 400 GeV. This implies a reheating temperature of T R ≈ O(10 7 ) GeV in order to meet the observed relic density in Eq. (1.1).
For each parameter point, the physical mass has been computed from the soft-breaking parameters using SPheno 3.3.3 [53,54]. In both scenarios, the lightest neutralino is a pure bino, while the lighter stau is strongly mixed, the mixing angle corresponding to cos 2 θτ ≈ 0.42 and sin 2 θτ ≈ 0.58 for Scenario I and cos 2 θτ ≈ sin 2 θτ ≈ 0.50 for Scenario II.
Coming to the calculation of the relic density, the physical mass spectrum obtained from SPheno is handed over to micrOMEGAs 2.4.1 [12] using the SUSY Les Houches Accord 2 format [55,56]. micrOMEGAs then performs the numerical integration of the Boltzmann equation based on the annihilation cross section computed by CalcHEP [57]. We will use micrOMEGAs to compute the neutralino or stau relic density, respectively. For Scenario II, the gravitino relic density will then be obtained from Eqs. (2.2) and (2.3) once the gravitino mass and the reheating temperature have been fixed.
In Tab. II, we summarize the dominant annihilation and co-annihilation channels contributing to the total annihilation cross section σ ann entering the Boltzmann equation, Eq. (1.3). For both parameter configurations, stau-antistau annihilation into top-antitop pairs is the We also indicate the resulting physical masses of the lightest neutralino and the lighter stau. The values of the remaining physical masses are not displayed here, as they are not relevant for our study. The gravitino mass for the study of Scenario II will be specified in Sec. IV C. All dimensionful quantities are given in GeV.
TABLE II. Relative contributions in percent of the dominant annihilation channels contributing to the annihilation cross section σann in the two reference scenarios I and II defined in Tab. I. Here, and ν denote arbitrary lepton and neutrino states, = e, µ, τ and ν = νe, νµ, ντ . Further contributions below 5% are omitted.
dominant contribution, followed by channels which, at the one-loop level, are insensitive to QCD corrections, such as processes including neutralinos, Higgs and gauge bosons, photons, leptons, and neutrinos. 1 Generally, top quarks are more important than other quarks in the final state due to the important top-quark Yukawa coupling, which is additionally tan β enhanced in exchanges of scalar Higgs bosons, h 0 and H 0 . Note that, since only the lighter stau is relevant in the given context, the exchange of a pseudoscalar Higgs A 0 is absent. We therefore focus on the annihilation into top quarks, i.e. we consider the processτ 1τ * 1 → tt. Providing QCD corrections to this channel means that we can correct about 32% and 26% of the total annihilation cross section for Scenario I and II, respectively.
We conclude the phenomenological discussion by illustrating the situation in the vicinity of our chosen reference scenario I. To this end, we show in Fig. 1 the regions corresponding to a relic density compatible with the range given in Eq. (1.1) obtained from the variation of the bino mass parameter M 1 and the stau mass parameter Mτ L around the values given in Table I. We also indicate the mass difference between the neutralino and the stau, and the relative contribution of the stau-antistau annihilation into top-antitop pairs which dominates this parameter region and will therefore be in the focus of our study.
As can be seen in Fig. 1, the parameter region where the relic density agrees with the limits of Eq. (1.1) closely follows the line where the neutralino and the stau are equal in mass. This illustrates the importance of coannihilations in order to obtain the observed relic density. In our Scenario I, indicated in Fig. 1 by the red dot, the small mass difference enhances the importance of the stau annihilations through the exponential factor given in Eq.
(1.5) and enables the neutralino relic density to be within the given limits.
In the remainder of this paper, we will present a higher-τ order calculation of the stau-antistau annihilation into top quarks and show the impact on the phenomenology discussed here.

III. CALCULATION DETAILS
In the present work, we focus on the annihilation of a stau-antistau pair into a top-antitop pair. At the treelevel, this process proceeds through the exchange of CPeven Higgs-boson (h 0 or H 0 ), a Z 0 -boson, or a photon in the s-channel, as shown in Fig. 2. Due to the specific structure of the associated coupling with sfermions, the exchange of a pseudoscalar Higgs boson A 0 does not contribute in the present case of two identical stau mass eigenstates in the initial state.
In the following we will discuss higher-order corrections to the diagrams shown in Fig. 2. We will review the virtual and the real O(α s )-corrections, as well as the Sommerfeld enhancement due to multiple-photon exchanges between the initial state particles.

A. NLO corrections
Virtual corrections proportional to the strong coupling constant α s only arise for the final state vertex through the exchange of a gluon or a gluino between the quarks as shown in Fig. 3. In order to regulate the arising divergences, we have evaluated the loop integrals [60] appearing in the vertex diagrams in D = 4−2 dimensions using the dimensional reduction (DR) scheme, which preserves supersymmetry. The ultraviolet divergences introduced in the loop integrals are cancelled by a renormalization of the model parameters and fields.
The O(α s ) contributions are completed by the gluon radiation diagrams shown in Fig. 4, which cancel the infrared divergences introduced by the virtual corrections that include the exchange of a massless gluon. In order to combine the virtual and the real corrections, cancelling the infrared divergences, we make use of the dipole subtraction method [61,62]. This method is based on the construction of an auxiliary cross section σ A , which includes the information about the infrared divergent behaviour of the original cross section. Moreover, the auxiliary cross section is constructed such that the one-  particle gluon phase space is factorized from the three particle phase-space and can be integrated out analytically. Using the auxiliary cross section, the next-toleading order (NLO) cross section can be written as where σ LO is the leading-order cross section, dσ R and dσ V are the differential cross sections stemming from the real emission and vertex correction diagrams, respectively. The first integration in Eq. (3.1) is performed over the three-particle phase space corresponding to the real emission diagrams, and the second integration is performed over the two-particle phase space. For this last part, the auxiliary cross section is integrated analytically over the gluon phase space. For the explicit construction of the auxiliary cross section we refer the reader to our previous work [22]. As mentioned above, the ultraviolet divergences are removed by a proper redefinition of model parameters which requires a careful definition of these parameters, i.e., choosing a renormalization scheme. In our previous works [22,24,25] we have proposed and used a mixed on-shell and DR renormalization scheme of the MSSM which is well-suited for dark matter calculations. However, there is a certain ambiguity in choosing the renormalization scheme, which we want to demonstrate here by defining an alternative scheme. The new alternative renormalization scheme differs from our standard scheme just by having the top quark mass defined in the DR scheme. The alternative renormalization scheme is also particularly suitable to study variations of the renormalization scale as all relevant parameters are defined in the DR scheme. The differences resulting from using the two alternative schemes or from varying the renormalization scale will be discussed in Sec. IV B.

B. QED Sommerfeld corrections
In the limit of low relative velocities, as it is typical during freeze-out of dark matter, annihilating particles can exchange light mediators leading to the wellknown Sommerfeld enhancement. In our case this effect is caused by the exchange of multiple photons between the incoming stau anti-stau pair, see Fig. 5. For an each exchange of a photon the cross section is corrected by a factor proportional to (α/v rel ). With α ≈ v rel during freeze-out, this contribution becomes non-perturbative and thus has to be resummed to all orders of perturbation theory.
The Sommerfeld effect in electroweak theories has been discussed intensively in the literature [32][33][34][35][36][37][38] and was studied previously for QCD in the context of DM@NLO to which we refer for more details regarding the computation and implementation [26][27][28].
As stau annihilation occurs only via an s-channel exchange, the s-wave contribution dominates the squared matrix element. Therefore, we can factorize the corrected cross section in terms of the leading order contribution with the Sommerfeld factor S 0 . The latter is evaluated by solving the Schrödinger equation The solution is the Green's function G(r; E + iΓτ 1 ) = G(r, r = 0; E + iΓτ 1 ). γ E = 0.5772 indicates the Euler-Mascheroni constant, and are parameters whereas the latter originates from the one-loop β-function including all fermions f up to the scale of the typical momentum exchange. The typical scale is taken to be the Coulomb scale µ C , 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, Given the negligible scale dependence of α em , its running can be neglected for the calculation of the Sommerfeld factor. As the Coulomb potential given in Eq. (3.5) is scale independent by itself, this implies similarly a negligible contribution of the β-function. With the NLO contribution being suppressed by an additional factor of α/(4π), it has generally only a negligible effect on the correction. Hence, we performed our final calculation by including the Coulomb potential at leading order only. For further details on the numerical evaluation, we refer to our previous papers [26][27][28].

IV. NUMERICAL RESULTS
Let us now discuss the numerical impact of the corrections presented in Secs. III A and III B, first on the annihilation cross section itself, and then on the prediction for the relic density of dark matter. For this numerical study, we will rely on the two reference scenarios defined in Tab. I and discussed in Sec. II.
In order to compute the relic density including the corrections discussed above, our full NLO calculation including Sommerfeld corrections has been implemented in the DM@NLO package. In practice, the evaluation of the Boltzmann equation by micrOMEGAs uses cross sections computed by CalcHEP which are replaced in specific cases by the values obtained from the DM@NLO calculation. In this way, the processes included in the DM@NLO are taken into account in a consistent way throughout the calculation of the relic density and provide a more precise prediction of the relic density.

A. Annihilation cross section and its theoretical uncertainty
In Fig. 6 we show the stau-antistau annihilation cross section as a function of the center-of-mass momentum p cm for masses and couplings from the Scenario I of Tab. I. Given that in the Boltzmann equation the total cross section is thermally averaged, we also show the corresponding thermal distribution. The velocity distribution indicates the momentum range which is most relevant for the computation of the relic density.
The two different plots in Fig. 6 show the next-toleading-order annihilation cross section results for both renormalization schemes mentioned in Sec. III A. Let us first discuss the results using our standard DM@NLO renormalization scheme. We compare our results to the result from micrOMEGAs (black dashed line). We see that our leading-order (LO) result (orange dotted line) does not coincide with the micrOMEGAs cross section. One of the reasons is the different definition of the top quark mass. In the DM@NLOrenormalization scheme we use the physical on-shell top quark mass whereas micrOMEGAs uses the top quark mass in the DR scheme. The other reason is the difference in the Yukawa couplings due to the fact that micrOMEGAs uses effective couplings to include some higher-order corrections.
Including the NLO corrections decreases the cross section by about 9% as compared to the LO result, while the NLO cross section is about 7.4% larger compared to the micrOMEGAs result. The relative correction is fairly constant for a large span of the center-of-mass momentum p cm . On top of the next-to-leading-order SUSY-QCD corrections, we include also the electroweak Sommerfeld enhancement.
The Sommerfeld enhancement dominates the cross section for small relative velocities. For an attractive force such as the electromagnetic force between a stau and an anti-stau particle, the Sommerfeld enhancement increases the cross section (blue dash-dotted line in Fig.  6). The final correction to the leading-order cross section (red line in Fig. 6) after including both the SUSY-QCD NLO corrections and the electroweak Sommerfeld enhancement is relatively small given that both effects compensate each other.
In addition to the shift of the numerical result, including higher-order corrections leads to a better estimate of the theoretical uncertainty associated with the prediction. The theoretical uncertainty is the estimate of the contributions of higher orders that are not included in the actual calculation. There are several methods to estimate this uncertainty.
One possibility relies on the fact that the dependence on the renormalization scale introduced through the higher-order corrections would disappear if all orders in perturbation theory could be included. The dependence on the renormalization scale is gradually reduced by including higher-order corrections. That means in turn that the remaining dependence is an estimator for the scale-dependent parts of the missing higher-order

contributions.
We have investigated the dependence of the cross section on the variation of the renormalization scale for both renormalization schemes (for technical details see Ref. [28]). The results are shown in Fig. 7. In the left panel, we show the impact of the variation of the renormalization scale between µ = 0.5 TeV and µ = 2 TeV on the next-to-leading-order cross section calculated in the DM@NLO renormalization scheme. The leading-order cross section is completely insensitive to the scale variation and even the next-to-leading-order cross section is only mildly sensitive in this scheme. This is simple to understand as the most prominent parameter in this case, the top quark mass, is defined in the on-shell scheme which by definition removes the renormalization scale dependence related to the top quark mass from both the leading and next-to-leading-order cross sections. The dependence of the next-to-leading-order cross section on the renormalization scale comes from the scale dependent strong coupling constant which was first introduced by the SUSY-QCD one-loop corrections. But even this dependence is only mild due to the high scale of µ = 1 TeV which is natural for this process. At such high scales, the changes in α s due to the change in scale are very small.
In the case of the NLO calculation in the alternative renormalization scheme, the top quark mass was defined in the DR scheme which leads to larger sensitivity to the change in the renormalization scale. As one can see in the right panel of Fig. 7, the scale dependence in the alternative scheme is about 3% at leading order and is reduced to per mille level at NLO.
From the investigation of the dependence of the cross section on the renormalization scale, we might conclude that the theoretical uncertainty at NLO is smaller than three per mille. There are some caveats to this conclusion. First, as we have seen, the sensitivity to scale changes depends on the renormalization scheme where (1.1) as measured by the Planck satellite at the 1σ confidence level. The leading (nextto-leading) order relic density from both our renormalization schemes is denotes by blue (black) lines. The predictions in the DM@NLO (alternative) renormalization scheme are shown using the solid (dashed) lines. As in Fig. 1, the green contours indicate the relative contribution of the processτ1τ * 1 → tt to the total annihilation cross section, based on the micrOMEGAs calculation.
only pure MS or DR schemes exhibit the full sensitivity. The other important caveat is that varying the renormalization scale highlights only the size of the scaledependent part of the higher-order corrections. In order to highlight the shortcomings of the estimation of the theoretical uncertainty by varying the renormalization scale, we compare the changes in the cross sections due to different renormalization schemes. A renormalization scheme is specified by the definition of the model parameters and the corresponding definition of the model parameter counterterms. At leading order the different definition of parameters cause a large difference between calculations in different renormalization schemes. This can be seen either by comparing the left and right panels in Fig. 6 or by constructing the ratio of the leading-order cross sections as in Fig. 8. In our case the difference between the leading-order cross sections is larger than 30%. At next-to-leading order the counterterms compensate for the difference in parameter definitions. The only difference between the next-to-leading order calculations in different renormalization schemes comes from the use of different parameters in the one-loop corrections. This difference is of a higher order and can be used as an estimate of theoretical uncertainty. In our case the difference between the NLO predictions in our two schemes is only about 4-5%. The theoretical uncertainty defined in this way also reduces with every order included in the calculation and it takes into account not only terms sensitive to changes in the renormalization scale but all terms which depend on the model parameters e.g. the masses. This definition of theoretical uncertainties is also not without flaws. To truly assess the theoretical uncertainty one would ideally use many different renormalization schemes which are not always simple to define consistently for a given model (here the MSSM). And even if this was feasible, this approach as well as the previous one, cannot capture the presence of constant terms which can be determined only by an exact calculation of the higher-order corrections whose size we are trying to estimate.
We see that in order to be conservative, in our case we should choose the variation of the renormalization scheme to define the theoretical uncertainty. We then conclude that the leading-order cross section has an uncertainty of about 30% and the cross section including the next-to-leading order corrections has still an uncertainty of 4-5%. Using the DM@NLO renormalization scheme produces smaller higher-order corrections indicating quicker convergence of the perturbative series. This is one of the reasons why we adopt this scheme again in the following and we apply it to the relic density calculation, assuming neutralino or gravitino dark matter within the pMSSM.

B. Impact on the neutralino relic density
We first consider the case, where the lightest neutralino is the dark matter candidate and the second-lightest supersymmetric particle is the lighter stau. This situation corresponds to the mass spectrum of Scenario I defined in Tab. I.
In order to study the impact of the higher-order corrections on the relic density, we vary two key parameters, namely the bino mass parameter M 1 and the left-handed stau mass parameter Mτ L around the values specified in the Scenario I. The parameter region where the relic density satisfies the experimental constraint given by Eq.
(1.1) is shown in Fig. 9 as a yellow band. The band is determined using the micrOMEGAs relic density calculation of the cross section. The width of the band corresponds to one sigma experimental uncertainty. The blue solid and dashed lines in Fig. 9 denote the predictions for the correct relic density from the leading-order calculations in the DM@NLO and the alternative renormalization scheme, respectively. The band formed by these two lines denotes the theoretical uncertainty of the leadingorder relic density calculation. Similarly, the black solid and dashed lines correspond to the next-to-leading-order relic density prediction in the two schemes. We see that the NLO predictions from both renormalization schemes are very consistent with each other and the theoretical uncertainty of the relic density determination at NLO is very small in this scenario.
In Fig. 10 we first compare the final relic density prediction from micrOMEGAs and from the NLO calculation in the DM@NLO scheme in the plane of the soft mass parameters M 1 and Mτ L , then in the plane of the corresponding physical neutralino and stau masses. Note that, also for the right plot of Fig. 10, all other input parameters are fixed to the values given in Tab. I.
We observe that the shift between the favoured region based on the micrOMEGAs/CalcHEP calculation and the one based on our full calculation amounts to a shift of about 3 GeV for the neutralino mass (for fixed stau mass) or about 5 GeV for the stau mass (for fixed neutralino mass). Most importantly, the shift is much larger than the width of the respective band corresponding to the experimental uncertainty given in Eq. (1.1). This shows the importance of including higher-order corrections, and in particular the importance of including the Sommerfeld enhancement in the present situation. Gravitino relic density for as NLSP Scenario II (DM@NLO scheme) FIG. 11. Parameter regions in the mG-TR plane which are compatible with the Planck limits given in Eq. (1.1) for the case of gravitino dark matter, where the stau relic density has been computed using micrOMEGAs (yellow) and our full NLO and Sommerfeld corrected cross section in the DM@NLO scheme (orange). All other parameters are fixed to those given for Scenario II in Tab. I. The blue contours correspond to the gravitino relic density based on the micrOMEGAs calculation. The black lines indicate the relative non-thermal contribution in percent according to Eq. (2.3) to the total gravitino relic density, again based on the micrOMEGAs calculation.

C. Impact on the gravitino relic density
Here, we consider the case of the gravitino being the lightest supersymmetric particle. The second-lightest particle is the lighter stauτ 1 with the mass given in Tab. I for Scenario II. Let us recall that in this illustrative scenario, all other superpartners are rather heavy with masses of about 5 TeV to simplify the analysis.
In a similar way as above for neutralino dark matter, we illustrate in Fig. 11 the impact of our NLO and Sommerfeld corrections presented in Sec. III A on the favoured region of parameter space. As in the previous case, the shift between the micrOMEGAs calculation and our full calculation is more important than the Planck uncertainty given in Eq. (1.1).
Although the non-thermal contribution accounts for only about 80% of the gravitino relic density, and the process affected by the presented corrections accounts for only about 32% of the total stau annihilation crosssection, the observed shift is more important than the impact found for Scenario I. This is caused by a relatively large impact of the NLO corrections in this scenario. Here, contrary to Scenario I, the squarks and gluino are rather heavy, such that the corresponding gluino loop contribution (see Fig. 3) is suppressed. This contribution has an opposite sign with respect to the Standard Model top-gluon loop contribution. The compensation between the two is therefore reduced and the relative NLO con-tribution is more important amounting to about 70% in this scenario.
In this illustrative scenario, for a fixed value of the reheating temperature, the corrections account for a shift of about 50 GeV in the gravitino mass, which corresponds to a shift of about 10%. For a fixed gravitino mass of 450 about GeV, the reheating temperature needs to be multiplied by about a factor of two in order to still satisfy the Planck constraint.
Let us emphasize that in a situation where the stops and the gluino are closer in mass to the annihilating stau, the impact of the presented corrections is therefore expected to be reduced and similar to what has been observed in the analysis of our Scenario I.

V. CONCLUSION
We have discussed the impact of NLO SUSY-QCD corrections and the QED Sommerfeld enhancement on the cross section of stau-anti-stau annihilation into top quarks as well as their impact on the relic density in scenarios where this cross section is important. We have explored a scenario where the lightest neutralino is the lightest supersymmetric particle (LSP) and a dark matter candidate and the mass difference between the neutralino and the lighter stau is small which increases the importance of the stau annihilation. As the stau annihilations are also important in scenarios with gravitino dark matter, we have analyzed the impact of NLO corrections on the gravitino relic density in a typical scenario with gravitino dark matter.
We have analyzed different ways of defining the theoretical uncertainty. We have shown that the usual way of using the variation of the renormalization scale largely underestimates the theoretical uncertainty in this case. It is better estimated by looking at different renormalization schemes. We have shown that at leading order, the uncertainty of the relevant cross section is about 30% which translates into an uncertainty of about 5% on the relic density. This uncertainty largely reduces at nextto-leading order and we have demonstrated that at NLO the theoretical uncertainty is comparable with the experimental one.
We have demonstrated that in the studied cases the next-to-leading order corrections are important. They shift the region in the parameter space which corresponds to the experimentally determined relic density by more than the experimental uncertainty. Moreover, the theoretical uncertainty of the next-to-leading order prediction for the relic density remains below 1% making the NLO prediction for the relic density very precise.