The next-to-leading order (NLO) ($$ \mathcal{O} $$
O
($$ {\alpha}_s^3 $$
α
s
3
)) corrections for gluon fragmentation functions to a heavy quark-antiquark pair in 3$$ {P}_J^{\left[1,8\right]} $$
P
J
1
8
states are calculated within the NRQCD factorization. We use the integration-by-parts reduction and differential equations to semi-analytically calculate the 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 the matching procedure, infrared-safe short distance coefficients and $$ \mathcal{O} $$
O
($$ {\alpha}_s^2 $$
α
s
2
) perturbative NRQCD matrix elements ⟨$$ {\mathcal{O}}^3{P}_J^{\left[1,8\right]} $$
O
3
P
J
1
8
(3$$ {S}_1^{\left[8\right]} $$
S
1
8
)⟩ are obtained simultaneously. The NLO short distance coefficients are found to have significant corrections comparing with the LO ones.