  We present a local subtraction scheme for computing next-to-next-to-leading
order QCD corrections to the production of a massive quark-antiquark pair from
a colourless initial state. The subtraction terms are built following the
CoLoRFulNNLO method and refined in such a way that their integration gives rise
to compact, fully analytic expressions. All ingredients necessary for a
numerical implementation of our subtraction scheme are provided in detail. As
an example, we calculate the fully differential decay rate of the Standard
Model Higgs boson to massive bottom quarks at next-to-next-to-leading order
accuracy in perturbative QCD.
