  Within the framework of nonrelativistic QCD (NRQCD) factorization formalism,
we compute QCD next-to-next-to-leading order (NNLO) corrections to the helicity
amplitudes as well as the decay width of $Z\to H+\gamma$, where $H$ can be
$\eta_Q (Q=c,b), J/\psi$, or $\Upsilon$. In addition, we resum the
next-to-leading logarithms (NLL) of ${m_Z^2}/{m_Q^2}$ to all orders of
$\alpha_s$ for the leading-twist helicity amplitude by employing the light-cone
factorization approach. It is worth mentioning that we obtain the analytic
expressions of the truncated NLL at $\alpha_s^2$. We find that the
$\mathcal{O}(\alpha_s)$ corrections are around 10\% for $\eta_c$ and $\Upsilon$
productions, however insignificant for $J/\psi$ and $\eta_b$ productions. The
$\mathcal{O}(\alpha_s^2)$ corrections are moderate for charmonium production,
while very small for bottomonium production. Moreover, it is found that the NLL
resummation can considerably alter the NRQCD prediction, especially for
$J/\psi$ production. Combining the NRQCD and light-cone computation, we make
phenomenological predictions on the decay widths and branching fractions. In
addition, we investigate the dependence of the theoretical results on the heavy
quark mass, and find the branching fraction of $Z\to H+\gamma$ monotonically
decreases as $m_Q$ increases.
