  The next-to-next-to-next-to-leading order (N^3LO) Hamiltonian of potential
nonrelativistic QCD is derived. The complete matching of the Hamiltonian and
the contribution from the ultrasoft dynamical gluons relevant for perturbative
bound-state calculations is performed including one-, two-, and three-loop
contributions. The threshold expansion is used to disentangle and match
contributions of different scales in the effective-theory calculations. As a
physical application, the heavy-quarkonium spectrum is obtained at N^3LO for
the case of vanishing QCD beta function. Our results set the stage for a full
N^3LO analysis of the heavy-quarkonium system.
