  We compute the next-to-next-to-leading order (NNLO) QCD corrections to the
hadronic decay rates of the pseudoscalar quarkonia, at the lowest order in
velocity expansion. The validity of NRQCD factorization for inclusive
quarkonium decay process, for the first time, is verified to relative order
$\alpha_s^2$. As a byproduct, the renormalization group equation (RGE) of the
leading NRQCD 4-fermion operator ${\cal O}_1({}^1S_0)$ is also deduced to this
perturbative order. By incorporating this new piece of correction together with
available relativistic corrections, we find that there exists severe tension
between the state-of-the-art NRQCD predictions and the measured $\eta_c$
hadronic width, and in particular the branching fraction of $\eta_c\to
\gamma\gamma$. NRQCD appears to be capable of accounting for $\eta_b$ hadronic
decay to a satisfactory degree, and our most refined prediction is ${\rm
Br}(\eta_b\to \gamma\gamma) = (4.8\pm 0.7)\times 10^{-5}$.
