  The next-to-leading order (NLO) ($\mathcal{O}(\alpha_s^3)$) corrections for
gluon fragmentation functions to a heavy quark-antiquark pair in
${^{3}\hspace{-0.6mm}P_{J}^{[1,8]}}$ states are calculated within the NRQCD
factorization. We use integration-by-parts reduction and differential equations
to semi-analytically calculate fragmentation functions in full-QCD, and find
that infrared divergences can be absorbed by the NRQCD long distance matrix
elements. Thus, the NRQCD factorization conjecture is verified at two-loop
level via a physical process, which is free of artificial ultraviolet
divergences. Through matching procedure, infrared-safe short distance
coefficients and $\mathcal{O}(\alpha_s^2)$ perturbative NRQCD matrix elements
$\langle{\mathcal
O}^{{^{3}\hspace{-0.6mm}P_{J}^{[1,8]}}}({^{3}\hspace{-0.6mm}S_{1}^{[8]}})\rangle_{\mathrm{NLO}}$
are obtained simultaneously. The NLO short distance coefficients are found to
have significant corrections comparing with the LO ones.
