  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}$.
