  The rare decay channels of Higgs boson to heavy quarkonium offer vital
opportunities to explore the coupling of Higgs to heavy quarks. We study the
semi-exclusive decay channels of Higgs boson to heavy quarkonia, i.e., $H^0\to
|(Q\bar{Q^{\prime}})[n]\rangle+\bar{Q}Q^{\prime}$ ($Q^{(\prime)}=c~\text{or}~b$
quark) within the NRQCD framework. In addition to the lower-level Fock states
$|(Q\bar{Q'})[1S]\rangle$ continent, contributions of high excited states
$|(Q\bar{Q'})[2S]\rangle$, $|(Q\bar{Q'})[3S]\rangle$,
$|(Q\bar{Q'})[4S]\rangle$, $|(Q\bar{Q'})[1P]\rangle$,
$|(Q\bar{Q'})[2P]\rangle$, $|(Q\bar{Q'})[3P]\rangle$ and
$|(Q\bar{Q'})[4P]\rangle$ are also studied. According to our study, the
contributions of high excited Fock states should be considered seriously.
Differential distributions of total decay width with respect to invariant-mass
and angles, as well as uncertainties caused by non-perturbative hadronic
non-perturbative matrix elements are discussed. If all excited heavy quarkonium
states decay to the ground spin-singlet state through electromagnetic or
hadronic interactions, we obtain the decay widths for $|(Q\bar{Q'})\rangle$
quarkonium production through $H^0$ semi-exclusive decays:
$25.10^{+11.6\%}_{-51.6\%}$ keV for $|(b\bar{c})[n]\rangle$ meson,
$3.23^{+0\%}_{-62.2\%}$ keV for $|(c\bar{c})[n]\rangle$ and
$2.36^{+0\%}_{-57.1\%}$ keV for $|(b\bar{b})[n]\rangle$, where uncertainties
are caused by adopting different non-perturbative potential models. At future
high energy LHC ($\sqrt{s}=27$ TeV), numerical results show that sizable
amounts of events for those high excited states can be produced, which implies
that one could also consider exploring the coupling properties of Higgs to
heavy quarks in these high excited states channels, especially for the
charmonium and bottomonium.
