  The next-to-next-to-leading order (NNLO) pQCD correction to the inclusive
decays of the heavy quarkonium $\eta_Q$ ($Q$ being $c$ or $b$) has been done in
the literature within the framework of nonrelativistic QCD. One may observe
that the NNLO decay width still has large conventional renormalization scale
dependence due to its weaker pQCD convergence, e.g. about $(^{+4\%}_{-34\%})$
for $\eta_c$ and $(^{+0.0}_{-9\%})$ for $\eta_b$, by varying the scale within
the range of $[m_Q, 4m_Q]$. The principle of maximum conformality (PMC)
provides a systematic way to fix the $\alpha_s$-running behavior of the
process, which satisfies the requirements of renormalization group invariance
and eliminates the conventional renormalization scheme and scale ambiguities.
Using the PMC single-scale method, we show that the resultant PMC conformal
series is renormalization scale independent, and the precision of the $\eta_Q$
inclusive decay width can be greatly improved. Taking the relativistic
correction $\mathcal{O}(\alpha_{s}v^2)$ into consideration, the ratios of the
$\eta_{Q}$ decays to light hadrons or $\gamma\gamma$ are: $R^{\rm
NNLO}_{\eta_c}|_{\rm{PMC}}=(3.93^{+0.26}_{-0.24})\times10^3$ and $R^{\rm
NNLO}_{\eta_b}|_{\rm{PMC}}=(22.85^{+0.90}_{-0.87})\times10^3$, respectively.
Here the errors are for $\Delta\alpha_s(M_Z) = \pm0.0011$. As a step forward,
by applying the Pad$\acute{e}$ approximation approach (PAA) over the PMC
conformal series, we obtain approximate NNNLO predictions for those two ratios,
e.g. $R^{\rm NNNLO}_{\eta_c}|_{\rm{PAA+PMC}} =(5.66^{+0.65}_{-0.55})\times10^3$
and $R^{\rm NNNLO}_{\eta_b}|_{\rm{PAA+PMC}}=(26.02^{+1.24}_{-1.17})\times10^3$.
The $R^{\rm NNNLO}_{\eta_c}|_{\rm{PAA+PMC}}$ ratio agrees with the latest PDG
value $R_{\eta_c}^{\rm{exp}}=(5.3_{-1.4}^{+2.4})\times10^3$, indicating the
necessity of a strict calculation of NNNLO terms.
