  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 $\alpha_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.
