FKS subtraction for quarkonium production at NLO

We extend the local infrared-divergence subtraction formalism, originally proposed by Frixione, Kunszt and Signer (FKS), to calculate short-distance (differential) cross section for any inclusive process involving a quarkonium particle in non-relativistic QCD (NRQCD) factorisation at next-to-leading order (NLO) accuracy in the strong coupling constant α_s. The new formulas are generally applicable to the production of an S- or P-wave quarkonium state in association with any number of elementary particles. The main new ingredients derived in this paper are the local and integrated soft counterterms for the colour-singlet and colour-octet P-wave bound states. It, therefore, paves the way to the automation of the NLO calculations for heavy quarkonium inclusive and associated production processes.