  Recently we have proved factorization of infrared divergences in NRQCD S-wave
heavy quarkonium production at high energy colliders at all orders in coupling
constant. One of the problem which still exists in the higher order pQCD
calculation of color singlet P-wave heavy quarkonium production/anihillation is
the appearance of non-canceling infrared divergences due to real soft gluons
exchange, although no such infrared divergences are present in the color
singlet S-wave heavy quarkonium. In this paper we find that since the
non-perturbative matrix element of the color singlet P-wave heavy quarkonium
production contains derivative operators, the gauge links are necessary to make
it gauge invariant and be consistent with the factorization of such
non-canceling infrared divergences at all orders in coupling constant.
