Based on the framework of nonrelativistic Quantum Chromodynamics (NRQCD), we carry out next-to-leading order (NLO) QCD corrections to the decay of $Z$ boson into $chi_c$ and $chi_b$, respectively. The branching ratio of $Z to chi_{c}(chi_b)+X$ is about $10^{-5}(10^{-6})$. It is found that, for $Z to chi_c(chi_b)+X$, the single gluon fragmentation diagrams of $^3S_1^{[8]}$, which first appear at the NLO level, can provide significant contributions, leading to a great enhancement on the leading-order results. Consequently the contributions from the color octet (CO) channels will account for a large proportion of the total decay widths. Moreover, the introduction of the CO processes will thoroughly change the color singlet (CS) predictions on the ratios of $Gamma_{chi_{c1}}/Gamma_{chi_{c0}}$, $Gamma_{chi_{c2}}/Gamma_{chi_{c0}}$, $Gamma_{chi_{b1}}/Gamma_{chi_{b0}}$ and $Gamma_{chi_{b2}}/Gamma_{chi_{b0}}$, which can be regarded as an outstanding probe to distinguish the CO and CS mechanism. With regard to the CS ($^3P_J^{[1]}$) channels, the heavy quark pair associated processes serve as the leading role, however, in the case of $chi_b$, $Z to bbar{b}[^3P_J^{[1]}]+g+g$ can also contribute significantly. Summing over all the feeddown contributions from $chi_{cJ}$ and $chi_{bJ}$, respectively, we find $Gamma(Z to J/psi+X)|_{chi_c-textrm{feeddown}}=(0.28 - 2.4) times 10^{-5}$ and $Gamma(Z to Upsilon(1S)+X)|_{chi_b-textrm{feeddown}}=(0.15 - 0.49) times 10^{-6}$.