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.