  The dark matter abundance plays a crucial role in the determination of the
valid parameter space of models both in the case of a discovery of dark matter
and in the context of exclusion limits. Reliable theoretical predictions of the
dark matter relic density require technically demanding precision calculations,
which were so far limited in their automation due to challenges in the
treatment of infrared divergences appearing in higher order calculations. In
particular, massive initial states need to be considered in early Universe
computations, so that the known dipole subtraction methods could not be
directly exploited. We therefore provide a full generalization of the dipole
subtraction method by Catani and Seymour to supersymmetric (SUSY) QCD with
massive initial states. All dipole splitting functions and their integrated
counterparts are given explicitly for four different dimensional schemes. To
showcase their application, we apply our results to dark matter
(co)annihilation processes in the context of the minimal supersymmetric
Standard Model. We also demonstrate the accuracy of the dipole method by
comparing our numerical results with those obtained with the phase space
slicing method. Our analytical results will facilitate future automation of
dark matter abundance calculations at next-to-leading order for both SUSY and
non-SUSY models.
