  We derive the imaginary part of the potential NRQCD Hamiltonian up to order
1/m^4, when the typical momentum transfer between the heavy quarks is of the
order of Lambda_{QCD} or greater, and the binding energy E much smaller than
Lambda_{QCD}. We use this result to calculate the inclusive decay widths into
light hadrons, photons and lepton pairs, up to O(mv^3 x
(Lambda_{QCD}^2/m^2,E/m)) and O(mv^5) times a short-distance coefficient, for
S- and P-wave heavy quarkonium states, respectively. We achieve a large
reduction in the number of unknown non-perturbative parameters and, therefore,
we obtain new model-independent QCD predictions. All the NRQCD matrix elements
relevant to that order are expressed in terms of the wave functions at the
origin and six universal non-perturbative parameters. The wave-function
dependence factorizes and drops out in the ratio of hadronic and
electromagnetic decay widths. The universal non-perturbative parameters are
expressed in terms of gluonic field-strength correlators, which may be fixed by
experimental data or, alternatively, by lattice simulations. Our expressions
are expected to hold for most of the charmonium and bottomonium states below
threshold. The calculations and methodology are explained in detail so that the
evaluation of higher order NRQCD matrix elements in this framework should be
straightforward. An example is provided.
