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 b\bar{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}$.