  In our previous work [2], it has been found that sizable charmonium events
via the channel $e^+e^- \to \gamma^*/Z^0 \to H(|c\bar{c}\rangle) +\gamma$ can
be produced at the suggested super $Z$ factory, where $H(|c\bar{c}\rangle)$
represents the dominant color-singlet $S$-wave and $P$-wave charmonium states
$J/\psi$, $\eta_c$, $h_c$ and $\chi_{cJ}$ ($J=0, 1, 2$), respectively. As an
important step forward, in the present paper, we present a next-to-leading
order (NLO) QCD analysis within the framework of nonrelativistic QCD. In
different to the case of $B$ factory in which the single charmonium production
is dominated by the channel via a virtual photon, at the super $Z$ factory, its
cross-section is dominated by the channel via a $Z^0$ boson. Estimations up to
NLO level are done at both the $B$ factory and the super $Z$ factory
experimental conditions. We observe that the NLO distributions have the same
shapes as those of LO distributions, but their differences are sizable. This
indicates that a NLO calculation is necessary and important to achieve a more
accurate estimation. Due to $Z^0$ boson resonance effect, at the super $Z$
factory with a high luminosity up to $10^{36}{\rm cm}^{-2}{\rm s}^{-1}$, when
summing all the color-singlet states' contribution together, one may observe
about $8.0\times10^{4}$ charmonium events via the channel $e^+e^-\to Z^0 \to
H(|c\bar{c}\rangle)+\gamma$ in one operation year. Then, such super $Z$ factory
could provide another useful platform to study the charmonium properties, even
for the higher charmonium states.
