  In this paper, we calculate the decay widths and branching fractions for the
decays $Z \to H+ \gamma$ up to ${\cal O}(\alpha_{s} v^2)$ accuracy within the
framework of nonrelativistic QCD, where $H$ stands for the $S$-wave quarkonium
$\eta_c$, $J/\psi$, $\eta_b$ or $\Upsilon$, respectively. To compare with the
leading-order terms, those corrections show good perturbative behavior as
expected. It is found that contributions from the next-to-leading order QCD
correction ${\cal O}(\alpha_{s}v^0)$, the relativistic correction ${\cal
O}(\alpha^{0}_{s}v^2)$ and their joint correction ${\cal O}(\alpha_{s} v^2)$
are sizable and comparable to each other, especially for the charmonium case.
Thus we need to take all of them into consideration for a sound estimation. For
a high luminosity electron-positron collider running around the $Z$-pole, due
to $Z$-boson resonance effect, sizable events could be produced from those rare
decay channels.
