  A rigorous QCD analysis of the inclusive annihilation decay rates of heavy
quarkonium states is presented. The effective-field-theory framework of
nonrelativistic QCD is used to separate the short-distance scale of
annihilation, which is set by the heavy quark mass $M$, from the
longer-distance scales associated with quarkonium structure. The annihilation
decay rates are expressed in terms of nonperturbative matrix elements of
4-fermion operators in nonrelativistic QCD, with coefficients that can be
computed using perturbation theory in the coupling constant $\alpha_s(M)$. The
matrix elements are organized into a hierarchy according to their scaling with
$v$, the typical velocity of the heavy quark. An analogous factorization
formalism is developed for the production cross sections of heavy quarkonium in
processes involving momentum transfers of order $M$ or larger. The
factorization formulas are applied to the annihilation decay rates and
production cross sections of S-wave states, up to corrections of relative order
$v^3$, and of P-wave states, up to corrections of relative order $v^2$.
