  In the paper, we calculate the fragmentation functions for $c \to \eta_c$ and
$b \to \eta_b$ up to next-to-leading-order (NLO) QCD accuracy. The ultraviolet
divergences in the real corrections are removed through operator
renormalization under the modified minimal subtraction scheme. We then obtain
the fragmentation functions $D_{c \to \eta_c}(z,\mu_F)$ and $D_{b \to
\eta_b}(z,\mu_F)$ up to NLO QCD accuracy, which are presented as figures and
fitting functions. The numerical results show that the NLO corrections are
significant. The sensitives of the fragmentation functions to the
renormalization scale and the factorization scale are analyzed explicitly.
