  We compute the full O(alpha_s) supersymmetric QCD corrections for
neutralino-stop co-annihilation into electroweak gauge and Higgs bosons in the
Minimal Supersymmetric Standard Model (MSSM). We show that these annihilation
channels are phenomenologically relevant within the so-called phenomenological
MSSM, in particular in the light of the observation of a Higgs-like particle
with a mass of about 126 GeV at the LHC. We present in detail our calculation,
including the renormalization scheme, the infrared treatment, and the
kinematical subtleties to be addressed. Numerical results for the
co-annihilation cross sections and the predicted neutralino relic density are
presented. We demonstrate that the impact of including the corrections on the
cosmologically preferred region of parameter space is larger than the current
experimental uncertainty from WMAP data.
