The dark matter abundance plays a crucial role in the determination of the
valid parameter space of models both in 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 (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 (SUSY) Standard Model (MSSM). We also
demonstrate the accuracy of the dipole method by comparing our numerical
results with those obtained with the space 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.Comment: 42 pages, 8 figures, 3 table