  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.
