For particle physics observables at colliders such as the LHC at CERN, it has been common practice for many decades to estimate the theoretical uncertainty by studying the variations of the predicted cross sections with a priori unpredictable scales. In astroparticle physics, this has so far not been possible, since most of the observables were calculated at Born level only, so that the renormalization scheme and scale dependence could not be studied in a meaningful way. In this paper, we present the first quantitative study of the theoretical uncertainty of the neutralino dark matter relic density from scheme and scale variations. We first explain in detail how the renormalization scale enters the tree-level calculations through coupling constants, masses and mixing angles. We then demonstrate a reduction of the renormalization scale dependence through one-loop SUSY-QCD corrections in many different dark matter annihilation channels and enhanced perturbative stability of a mixed on-shell/$\bar{\rm DR}$ renormalization scheme over a pure $\bar{\rm DR}$ scheme in the top-quark sector. In the stop-stop annihilation channel, the Sommerfeld enhancement and its scale dependence are shown to be of particular importance. Finally, the impact of our higher-order SUSY-QCD corrections and their scale uncertainties are studied in three typical scenarios of the phenomenological Minimal Supersymmetric Standard Model with eleven parameters (pMSSM-11). We find that the theoretical uncertainty is reduced in many cases and can become comparable to the size of the experimental one in some scenarios.

We compute the full O(alpha_s) SUSY-QCD corrections to dark matter annihilation in the Higgs-funnel, resumming potentially large mu tan beta and A_b contributions and keeping all finite O(m_b,s,1/tan^2 beta) terms. We demonstrate numerically that these corrections strongly influence the extraction of SUSY mass parameters from cosmological data and must therefore be included in common analysis tools such as DarkSUSY or micrOMEGAs.