  The principle of maximum conformality (PMC) provides a way to eliminate the
conventional renormalization scale ambiguity in a systematic way. By applying
the PMC scale setting, all non-conformal terms in perturbative series are
summed into the running coupling, and one obtains a unique, scale-fixed
prediction at any finite order. In the paper, we make a detailed PMC analysis
for the spin-singlet heavy quarkoniums decay (into light hadrons) at the
next-to-leading order. After applying the PMC scale setting, the decay widths
for all those cases are almost independent of the initial renormalization
scales. The PMC scales for $\eta_c$ and $h_c$ decays are below $1$ GeV, in
order to achieve a confidential pQCD estimation, we adopt several low-energy
running coupling models to do the estimation. By taking the MPT model, we
obtain: $\Gamma(\eta_{c} \to LH)=25.09^{+5.52}_{-4.28}$ MeV, $\Gamma(\eta_{b}
\to LH)=14.34^{+0.92}_{-0.84}$ MeV, $\Gamma(h_{c} \to LH)=0.54^{+0.06}_{-0.04}$
MeV and $\Gamma(h_{b} \to LH)=39.89^{+0.28}_{-0.46}$ KeV, where the errors are
calculated by taking $m_{c}\in[1.40\rm GeV,1.60\rm GeV]$ and $m_{b}\in[4.50\rm
GeV,4.70\rm GeV]$. These decay widths agree with the principle of minimum
sensitivity estimations, in which the decay widths of $\eta_{c,b}$ are also
consistent with the measured ones.
