  In this paper, we study the $J/\psi(\Upsilon)$ production in Z boson decay in
color-singlet model(CSM). We calculate the next-to-leading-order (NLO) QCD
correction to $Z \to Quarkonium+Q\bar{Q}$, the dominant contribution in the
CSM, with the vector and axial-vector parts in $ZQ\bar{Q}$ vertex being treated
separately. The results show that the vector and axial-vector parts have the
same K factor (the ratio of NLO result to leading-order result) 1.13 with the
renormalization scale $\mu$=2$m_c$ and $m_c=1.5GeV$, and the K factor falls to
0.918 when applying the Brodsky, Lepage, and Mackenzie(BLM) renormalization
scale scheme with obtained $\mu_{BLM}=2.28GeV$ and $m_c=1.5$GeV. By including
the contributions from the next-dominant ones, the photon and gluon
fragmentation processes, the branching ratio for $Z \to J/\psi_{prompt}+X$ is
$(7.3 \sim 10.0)\times 10^{-5}$ with the uncertainty consideration for the
renormalization scale and Charm quark mass. The results are about half of the
central value of the experimental measurement 2.1$\times10^{-4}$. Furthermore,
the $J/\psi$ energy distribution in our calculation is not well consistent with
the experimental data. Therefore, even at QCD NLO, the contribution to $Z \to
J/\psi_{prompt}+X$ from the CSM can not fully account for the experimental
measurement. And there should be contributions from other mechanisms, such as
the color-octet(COM) contributions. We define $R_{c\bar{c}}=\frac{\Gamma(Z \to
J/\psi c\bar{c}X)}{\Gamma(Z \to J/\psi X)}$ and obtain $R_{cc}=0.84$ for only
CSM contribution and $R_{cc}=0.49$ for COM and CSM contributions together. Then
$R_{cc}$ measurement could be used to clarify the COM contributions.
